(**
  Formalization of real valued interval arithmetic
  Used in soundness proofs for error bound validator.
**)
Require Import Coq.Reals.Reals Coq.micromega.Psatz Coq.QArith.Qreals.
Require Import Daisy.Infra.Abbrevs Daisy.Infra.RealSimps.

(**
  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.
**)
Notation valid intv :=
  (IVlo intv <= IVhi intv)%R.

Notation contained x intv :=
  (IVlo intv <= x <= IVhi intv)%R.

Lemma contained_implies_valid (a:R) (iv:interval) :
  contained a iv ->
  valid iv.
Proof.
  intros inIntv; apply (Rle_trans _ a _); destruct inIntv; auto.
Qed.

(**
  Subset definition; when an interval is a subinterval of another
**)
Notation isSupersetInterval iv1 iv2 :=
  ((IVlo iv2 <= IVlo iv1)%R /\ (IVhi iv1 <= IVhi iv2)%R).

Lemma contained_superset a iv1 iv2:
  contained a iv1 ->
  isSupersetInterval iv1 iv2 ->
  contained a iv2.
Proof.
  intros; simpl in *. lra.
Qed.

Definition pointInterval (x:R) := mkInterval x x.

Lemma contained_implies_subset (a:R) (iv:interval):
  contained a iv ->
  isSupersetInterval (pointInterval a) iv.
Proof.
  intros containedIntv.
  unfold pointInterval in *; destruct containedIntv; split; auto.
Qed.

Lemma validPointInterval (a:R) :
  contained a (pointInterval a).
Proof.
  split; simpl; apply Req_le; auto.
Qed.

(**
  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)).

Lemma min4_correct (a b c d:R) :
  (let m := min4 a b c d in
   m <= a /\ m <= b /\ m <= c /\ m <= d)%R.
Proof.
  intros m.
  repeat split; unfold m, min4.
  - apply Rmin_l.
  - assert (forall c,  Rmin b c <= b)%R by (apply Rmin_l).
    apply (Rle_trans _ (Rmin b (Rmin c d)) _); auto.
    apply Rmin_r.
  - assert (Rmin c d <= c)%R by (apply Rmin_l).
    assert (Rmin b (Rmin c d) <= c)%R.
    + apply (Rle_trans _ (Rmin c d) _); auto. apply Rmin_r.
    + apply (Rle_trans _ (Rmin b (Rmin c d)) _); auto. apply Rmin_r.
  - assert (Rmin c d <= d)%R by (apply Rmin_r).
    assert (Rmin b (Rmin c d) <= d)%R.
    + apply (Rle_trans _ (Rmin c d) _); auto. apply Rmin_r.
    + apply (Rle_trans _ (Rmin b (Rmin c d)) _); auto. apply Rmin_r.
Qed.

Lemma max4_correct (a b c d:R) :
  (let m := max4 a b c d in
   a <= m /\ b <= m /\ c <= m /\ d <= m)%R.
Proof.
  intros m.
  repeat split; unfold m, max4.
  - apply Rmax_l.
  - assert (forall c,  b <= Rmax b c)%R by (apply Rmax_l).
    apply (Rle_trans _ (Rmax b (Rmax c d)) _); auto.
    apply Rmax_r.
  - assert (c <= Rmax c d)%R by (apply Rmax_l).
    assert (c <= Rmax b (Rmax c d))%R.
    + apply (Rle_trans _ (Rmax c d) _); auto. apply Rmax_r.
    + apply (Rle_trans _ (Rmax b (Rmax c d)) _); auto. apply Rmax_r.
  - assert (d <= Rmax c d)%R by (apply Rmax_r).
    assert (d <= Rmax b (Rmax c d))%R.
    + apply (Rle_trans _ (Rmax c d) _); auto. apply Rmax_r.
    + apply (Rle_trans _ (Rmax b (Rmax c d)) _); auto. apply Rmax_r.
Qed.

(**
 Asbtract interval update function, parametric by actual operator applied.
**)
Definition absIntvUpd (op:R->R->R) (iv1:interval) (iv2:interval) :=
  (
    min4 (op (IVlo iv1) (IVlo iv2))
         (op (IVlo iv1) (IVhi iv2))
         (op (IVhi iv1) (IVlo iv2))
         (op (IVhi iv1) (IVhi iv2)),
    max4 (op (IVlo iv1) (IVlo iv2))
         (op (IVlo iv1) (IVhi iv2))
         (op (IVhi iv1) (IVlo iv2))
         (op (IVhi iv1) (IVhi iv2))
  ).

Definition widenInterval (Iv:interval) (v:R) := mkInterval (IVlo Iv - v) (IVhi Iv + v).

Definition negateInterval (iv:interval) := mkInterval (- IVhi iv) (- IVlo iv).

Lemma interval_negation_valid (iv:interval) (a:R) :
  contained a iv -> contained (- a) (negateInterval iv).
Proof.
  destruct iv as [ivlo ivhi]; simpl; intros [ivlo_le_a a_le_ivhi].
  split; apply Ropp_ge_le_contravar; apply Rle_ge; auto.
Qed.

Definition invertInterval (iv:interval) := mkInterval (/ IVhi iv) (/ IVlo iv).

Lemma interval_inversion_valid (iv:interval) (a:R) :
  (IVhi iv < 0 \/ 0 < IVlo iv -> contained a iv -> contained (/ a) (invertInterval iv))%R.
Proof.
  destruct iv as [ivlo ivhi]; simpl;
    intros [ivhi_lt_zero | zero_lt_ivlo]; intros [ivlo_le_a a_le_ivhi];
      assert (ivlo <= ivhi)%R by (eapply Rle_trans; eauto);
      split.
  - assert (- / a <= - / ivhi)%R.
    + assert (0 < - ivhi)%R by lra.
      repeat rewrite Ropp_inv_permute; try lra.
      eapply RIneq.Rinv_le_contravar; try lra.
    + apply Ropp_le_ge_contravar in H0;
        repeat rewrite Ropp_involutive in H0; lra.
  - assert (- / ivlo <= - / a)%R.
    + repeat rewrite Ropp_inv_permute; try lra.
       eapply RIneq.Rinv_le_contravar; try lra.
    + apply Ropp_le_ge_contravar in H0;
        repeat rewrite Ropp_involutive in H0; lra.
  - eapply RIneq.Rinv_le_contravar; try lra.
  - eapply RIneq.Rinv_le_contravar; try lra.
Qed.


Definition addInterval (iv1:interval) (iv2:interval) :=
  absIntvUpd Rplus iv1 iv2.

Lemma interval_addition_valid iv1 iv2 a b:
  contained a iv1 ->
  contained b iv2 ->
  contained (a + b) (addInterval iv1 iv2).
Proof.
  unfold 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 subtractInterval (iv1:interval) (iv2:interval) :=
  addInterval iv1 (negateInterval iv2).

Corollary interval_subtraction_valid iv1 iv2 a b:
  contained a iv1 ->
  contained b iv2 ->
  contained (a - b) (subtractInterval iv1 iv2).
Proof.
  unfold subtractInterval.
  intros contained_1 contained_iv2.
  rewrite Rsub_eq_Ropp_Rplus.
  apply interval_addition_valid; auto.
  apply interval_negation_valid; auto.
Qed.

Definition multInterval (iv1:interval) (iv2:interval) :=
  absIntvUpd Rmult iv1 iv2.

(**
  Prove validity of multiplication for the interval lattice.
  Prove is structurally the same as the proof done Jourdan et al. in Verasco, in FloatIntervalsForward.v
  TODO: Credit
**)
Lemma interval_multiplication_valid (iv1:interval) (iv2:interval) (a:R) (b:R) :
  contained a iv1 ->
  contained b iv2 ->
  contained (a * b) (multInterval iv1 iv2).
Proof.
  unfold multInterval, absIntvUpd, IVlo, IVhi.
  destruct iv1 as [ivlo1 ivhi1]; destruct iv2 as [ivlo2 ivhi2]; simpl.
  intros [lo_leq_a a_leq_hi] [lo_leq_b b_leq_hi].
  pose proof (min4_correct (ivlo1 * ivlo2) (ivlo1 * ivhi2) (ivhi1 * ivlo2) (ivhi1 * ivhi2))
    as [leq_lolo [leq_lohi [leq_hilo leq_hihi]]];
    pose proof (max4_correct (ivlo1 * ivlo2) (ivlo1 * ivhi2) (ivhi1 * ivlo2) (ivhi1 * ivhi2))
    as [lolo_leq [lohi_leq [hilo_leq hihi_leq]]].
  split.
  (* Lower Bound *)
  - destruct (Rle_lt_dec a 0).
    + destruct (Rle_lt_dec ivhi2 0).
      * pose proof (Rmult_le_compat_neg_l a b ivhi2 r b_leq_hi) as ahi2_leq_ab.
        pose proof (Rmult_le_compat_neg_l ivhi2 a ivhi1 r0 a_leq_hi) as hihi_leq_ahi2.
        eapply Rle_trans.
        apply leq_hihi. rewrite Rmult_comm.
        eapply Rle_trans. apply hihi_leq_ahi2.
        rewrite Rmult_comm; auto.
      * eapply Rle_trans.
        apply leq_lohi.
        apply Rlt_le in r0.
        pose proof (Rmult_le_compat_l ivhi2 ivlo1 a).
        rewrite Rmult_comm.
        eapply (Rle_trans). apply H; auto.
        rewrite Rmult_comm.
        eapply (Rmult_le_compat_neg_l a); auto.
    + destruct (Rle_lt_dec ivlo2 0).
      * eapply Rle_trans.
        apply leq_hilo.
        rewrite Rmult_comm.
        eapply Rle_trans.
        apply (Rmult_le_compat_neg_l ivlo2 a ivhi1 r0 a_leq_hi).
        rewrite Rmult_comm.
        eapply (Rmult_le_compat_l a); auto.
        exact (Rlt_le _ _ r).
      * eapply Rle_trans.
        apply leq_lolo.
        apply Rlt_le in r; apply Rlt_le in r0.
        rewrite Rmult_comm.
        apply (Rle_trans _ (ivlo2 * a)).
        eapply (Rmult_le_compat_l ivlo2 ); auto.
        rewrite Rmult_comm.
        eapply (Rmult_le_compat_l a); auto.
  - destruct (Rle_lt_dec a 0).
    + eapply Rle_trans.
      apply (Rmult_le_compat_neg_l a ivlo2); auto.
      destruct (Rle_lt_dec ivlo2 0).
      * rewrite Rmult_comm.
        eapply Rle_trans.
        eapply (Rmult_le_compat_neg_l ivlo2 ivlo1); auto.
        rewrite Rmult_comm; auto.
      * rewrite Rmult_comm.
        eapply Rle_trans.
        eapply (Rmult_le_compat_l ivlo2 a ivhi1); auto.
        apply Rlt_le; auto.
        rewrite Rmult_comm; auto.
    + eapply Rle_trans.
      apply (Rmult_le_compat_l) with (r2 := ivhi2); auto.
      apply Rlt_le; auto.
      destruct (Rle_lt_dec ivhi2 0).
      * rewrite Rmult_comm. eapply Rle_trans.
        eapply (Rmult_le_compat_neg_l) with (r1:= ivlo1); auto.
        rewrite Rmult_comm; auto.
      * rewrite Rmult_comm; eapply Rle_trans.
        apply (Rmult_le_compat_l) with (r2 := ivhi1); auto.
        apply Rlt_le; auto.
        rewrite Rmult_comm; auto.
Qed.

Definition divideInterval (iv1:interval) (iv2:interval) :=
  multInterval iv1 (mkInterval (/ (IVhi iv2)) (/ (IVlo iv2))).

Corollary interval_division_valid a b iv1 iv2:
  (IVhi iv2 < 0 \/ 0 < IVlo iv2 -> contained a iv1 -> contained b iv2 -> contained (a / b) (divideInterval iv1 iv2))%R.
Proof.
  intros nodiv0 c_a c_b.
  unfold divideInterval.
  unfold Rdiv.
  eapply interval_multiplication_valid; auto.
  apply interval_inversion_valid; auto.
Qed.

Lemma contained_leq_maxAbs a iv:
  contained a iv ->
  (Rabs a <= RmaxAbsFun iv)%R.
Proof.
  intros contained_a.
  unfold RmaxAbsFun.
  simpl in contained_a; destruct contained_a.
  eapply RmaxAbs; auto.
Qed.

Lemma contained_leq_maxAbs_val a iv:
  contained a iv ->
  (a <= RmaxAbsFun (iv))%R.
Proof.
  intros contained_a; unfold RmaxAbsFun.
  assert (Rabs a <= RmaxAbsFun iv)%R by (apply contained_leq_maxAbs; auto).
  eapply Rle_trans.
  apply Rle_abs.
  auto.
Qed.

Lemma contained_leq_maxAbs_neg_val a iv:
  contained a iv ->
  (- a <= RmaxAbsFun (iv))%R.
Proof.
  intros contained_a; unfold RmaxAbsFun.
  assert (Rabs a <= RmaxAbsFun iv)%R by (apply contained_leq_maxAbs; auto).
  eapply Rle_trans.
  apply Rle_abs.
  rewrite Rabs_Ropp.
  auto.
Qed.

Lemma distance_gives_iv a b e iv:
  contained a iv ->
  (Rabs (a - b) <= e)%R ->
  contained b (widenInterval iv e).
Proof.
  intros contained_a abs_le; simpl in *.
  destruct contained_a as [lo_a a_hi].
  rewrite Rabs_minus_sym in abs_le.
  unfold Rabs in abs_le.
  destruct Rcase_abs in abs_le; try lra.
Qed.

Lemma minAbs_positive_iv_is_lo a b:
  (0 < a)%R ->
  (a <= b)%R ->
  RminAbsFun (a,b) = a.
Proof.
  intros; unfold RminAbsFun; simpl.
  assert (0 < b)%R by lra.
  assert (Rabs a = a)%R as Rabs_pos_a by (apply Rabs_right; lra).
  assert (Rabs b = b)%R as Rabs_pos_b by (apply Rabs_right; lra).
  rewrite Rabs_pos_a, Rabs_pos_b.
  rewrite Rmin_left; lra.
Qed.

Lemma minAbs_negative_iv_is_hi a b:
  (b < 0)%R ->
  (a <= b)%R ->
  (RminAbsFun (a,b) = - b)%R.
Proof.
  intros; unfold RminAbsFun; simpl.
  assert (Rabs a = - a)%R as Rabs_neg_a by (apply Rabs_left; lra).
  assert (Rabs b = - b)%R as Rabs_neg_b by (apply Rabs_left; lra).
  rewrite Rabs_neg_a, Rabs_neg_b.
  rewrite Rmin_right; lra.
Qed.