Library Interval.Tactics.Plot_helper

From Coq Require Import Reals List.

Require Import Sig.
Require Import Generic_proof.
Require Import Interval_helper.
Require Import Xreal Interval Tree Reify Prog.

Definition reify_var : R.

Definition plot2 (f : R → R) (ox dx oy dy : R) (h : Z) (l : list (Z × Z)) :=
  ∀ i x, (ox + dx × INR i ≤ x ≤ ox + dx × INR (S i))%R →
  (oy ≤ f x ≤ oy + dy × IZR h)%R →
  let r := nth i l (0%Z, h) in
  (oy + dy × IZR (fst r) ≤ f x ≤ oy + dy × IZR (snd r))%R.

Module PlotTacticAux (F : FloatOps with Definition radix := Zaux.radix2 with Definition sensible_format := true) (I : IntervalOps with Module F := F).

Module F' := FloatExt F.
Module IH := IntervalTacticAux I.
Import IH.

Definition plot1 (f : R → R) (ox dx : R) (l : list I.type) :=
  ∀ i x, (ox + dx × INR i ≤ x ≤ ox + dx × INR (S i))%R →
  contains (I.convert (nth i l I.nai)) (Xreal (f x)).

Lemma plot_ext :
  ∀ f g u d l,
  (∀ x, f x = g x) →
  plot1 f u d l → plot1 g u d l.

Fixpoint bound_plot_aux prec (fi : I.type → I.type) (ui di : I.type) (xz : Z) (i : nat) (acc : I.type) : I.type :=
  match i with
  | O ⇒ acc
  | S i ⇒
    let xz := Z.succ xz in
    let xi := I.add prec ui (I.mul prec di (I.fromZ prec xz)) in
    bound_plot_aux prec fi ui di xz i (I.join (fi xi) acc)
  end.

Definition bound_plot prec hyps pf cf oxi dxi nb :=
  let bounds := R.merge_hyps prec hyps ++ map (T.eval_bnd prec) cf in
  let fi xi := nth 0 (A.BndValuator.eval prec pf (xi :: bounds)) I.nai in
  bound_plot_aux prec fi oxi dxi 0%Z (Z.to_nat nb) (fi oxi).

Fixpoint sample_plot_aux prec (gi : I.type → I.type → I.type) (check : I.type → bool) (ui di zi2 : I.type) (fi : I.type → I.type) (i : nat) (mz xz rz : Z) (acc : list I.type) : list I.type :=
  match i with
  | O ⇒ acc
  | S i ⇒
    let xz' := Z.pred xz in
    let xi1 := I.add prec ui (I.mul prec di (I.fromZ prec xz')) in
    let xi2 := I.add prec ui (I.mul prec di (I.fromZ prec xz)) in
    let xi := I.join xi1 xi2 in
    let zi1 := fi xi1 in
    let yi := fi xi in
    let c := andb
      (orb (check (I.meet yi (I.lower_extent zi1))) (check (I.meet yi (I.upper_extent zi1))))
      (orb (check (I.meet yi (I.lower_extent zi2))) (check (I.meet yi (I.upper_extent zi2)))) in
    let yizi :=
      if c then (yi, zi1)
      else let fi := gi xi in (fi xi, fi xi1) in
    let mz' :=
      if Z.eqb mz xz' then (mz - Z.div2 (3 × rz))%Z
      else if c then mz
      else (xz' - 1 - Z.div2 rz)%Z in
    let mz' :=
      if Z.ltb mz' 0%Z then 0%Z else mz' in
    let firz :=
      if Z.eqb mz mz' then (fi, rz)
      else
        let xi0 := I.add prec ui (I.mul prec di (I.fromZ prec mz')) in
        (gi (I.join xi0 xi1), (xz' - mz'))%Z in
    let acc := fst yizi :: acc in
    sample_plot_aux prec gi check ui di (snd yizi) (fst firz) i mz' xz' (snd firz) acc
  end.

Definition sample_plot prec deg check hyps pf cf oxi dxi nb :=
  let hyps := R.merge_hyps prec hyps in
  let gi :=
    let bounds := hyps ++ map (T.eval_bnd prec) cf in
    let fi xi := nth 0 (A.BndValuator.eval prec pf (xi :: bounds)) I.nai in
    let bounds := A.TaylorValuator.TM.var :: map A.TaylorValuator.TM.const bounds in
    fun yi ⇒
      let zi := fi yi in
      let fi := nth 0 (A.TaylorValuator.eval prec deg yi pf bounds) A.TaylorValuator.TM.dummy in
      fun xi ⇒
        let zi' := A.TaylorValuator.TM.eval (prec, deg) fi yi xi in
        if I.subset xi yi then I.meet zi zi' else zi' in
  let ui := I.add prec oxi (I.mul prec dxi (I.fromZ prec 0)) in
  let vi := I.add prec oxi (I.mul prec dxi (I.fromZ prec nb)) in
  sample_plot_aux prec gi check oxi dxi I.whole (gi (I.join ui vi)) (Z.to_nat nb) 0%Z nb nb nil.

Lemma sample_plot_correct :
  ∀ prec deg check vars hyps pf cf oxi dxi ox dx nb l,
  contains (I.convert oxi) (Xreal ox) →
  contains (I.convert dxi) (Xreal dx) →
  sample_plot prec deg check hyps pf cf oxi dxi (Zpos nb) = l →
  eval_hyps hyps vars (
    plot1 (fun t ⇒ eval_real' pf (t :: vars) cf) ox dx l).

Definition clamp_lower (v : Basic.float Basic.radix2) (h : Z) :=
  match v with
  | Basic.Fzero ⇒ 0%Z
  | Basic.Fnan ⇒ 0%Z
  | Basic.Float true _ _ ⇒ 0%Z
  | Basic.Float false m e ⇒
    let v := Z.shiftl (Zpos m) e in
    if Z.leb h v then h else v
  end.

Definition clamp_upper (v : Basic.float Basic.radix2) (h : Z) :=
  match v with
  | Basic.Fzero ⇒ 0%Z
  | Basic.Fnan ⇒ h
  | Basic.Float true _ _ ⇒ 0%Z
  | Basic.Float false m e ⇒
    let v:=
      match e with
      | Z0 ⇒ Zpos m
      | Zpos e' ⇒ Z.shiftl (Zpos m) e
      | Zneg e' ⇒ Z.shiftl (Zpos m + (Z.shiftl 1 (Zpos e')) - 1) e
      end in
    if Z.leb h v then h else v
  end.

Definition clamp (xi : I.type) (h : Z) :=
  (clamp_lower (F.toF (I.lower xi)) h, clamp_upper (F.toF (I.upper xi)) h).

Theorem clamp_correct :
  ∀ xi h x,
  contains (I.convert xi) (Xreal x) →
  (0 ≤ x ≤ IZR h)%R →
  let yi := clamp xi h in
  (IZR (fst yi) ≤ x ≤ IZR (snd yi))%R.

Fixpoint clamp_plot prec (vi ei : I.type) (h : Z) (l : list I.type) : list (Z × Z) :=
  match l with
  | nil ⇒ nil
  | cons yi l ⇒
    let r := clamp (I.mul prec (I.sub prec yi vi) ei) h in
    cons r (clamp_plot prec vi ei h l)
  end.

Lemma affine_transf :
  ∀ oy dy y1 y2 y : R,
  (0 < dy)%R →
  (oy + dy × y1 ≤ y ≤ oy + dy × y2)%R ↔ (y1 ≤ (y - oy) / dy ≤ y2)%R.

Lemma clamp_plot_correct :
  ∀ prec oyi dyi f ox dx oy dy h l1 l2,
  (0 < dy)%R →
  contains (I.convert oyi) (Xreal oy) →
  contains (I.convert dyi) (Xreal (/dy)) →
  clamp_plot prec oyi dyi h l1 = l2 →
  plot1 f ox dx l1 →
  plot2 f ox dx oy dy h l2.

Definition get_bounds (prec : F.precision) (l : list I.type): F.type × F.type :=
  let yi :=
    match l with
    | cons hi l ⇒ List.fold_left I.join l hi
    | nil ⇒ I.empty
    end in
  
  let yl := I.lower yi in
  let yu := I.upper yi in
  let yw := F.sub_UP prec yu yl in
  (F.sub_DN prec yu yw, F.add_UP prec yl yw).

Ltac unify_eq native :=
  match goal with
  | |- ?f ?p1 = ?p2 ⇒
    match native with
    | true ⇒
      let p1 := eval hnf in p1 in
      let p := eval native_compute in (f p1) in
      instantiate (p2 := p) ;
      native_cast_no_check (eq_refl p2)
    | false ⇒
      let p1 := eval hnf in p1 in
      let p := eval vm_compute in (f p1) in
      instantiate (p2 := p) ;
      vm_cast_no_check (eq_refl p2)
    end
  end.

Ltac plot1_aux1 prec x1 x2 w h d native tac_b :=
  let x1 := reify x1 constr:(@nil R) in
  let x2 := reify x2 constr:(@nil R) in
  let ox := eval vm_compute in (I.lower (T.eval_bnd prec x1)) in
  let dx := eval vm_compute in (F.div_UP prec (F.sub_UP prec (I.upper (T.eval_bnd prec x2)) ox) (F.fromZ_DN prec (Zpos w))) in
  let oxr := eval cbv -[IZR Rdiv] in (proj_val (I.F.convert ox)) in
  let dxr := eval cbv -[IZR Rdiv] in (proj_val (I.F.convert dx)) in
  match goal with
  | |- plot1 ?f ?ox' ?dx' ?p ⇒
    unify ox' oxr ;
    unify dx' dxr ;
    let fapp := eval cbv beta in (f reify_var) in
    let vars := constr:((reify_var :: nil)%list) in
    let vars := get_vars fapp vars in
    let vars :=
      match get_vars fapp vars with
      | (reify_var :: ?vars)%list ⇒ vars
      end in
    eapply plot_ext ; [
      let t := fresh "t" in
      intros t ; hide_lhs ;
      let fapp := eval cbv beta in (f t) in
      reify_partial fapp (t :: vars) ;
      exact (fun H ⇒ H) |] ;
    find_hyps vars ;
    let y1y2 := tac_b prec ox dx w in
    let thr := eval vm_compute in (F.div_UP prec (F.sub_UP prec (snd y1y2) (fst y1y2)) (F.fromZ_DN prec (Zpos h))) in
    apply (sample_plot_correct prec) with
      (deg := d) (nb := w) (l := p)
      (check := fun yi ⇒ F'.le' (F.sub_UP prec (I.upper yi) (I.lower yi)) thr)
      (1 := I.singleton_correct ox)
      (2 := I.singleton_correct dx) ;
    unify_eq native
  end.

Ltac plot2_aux prec x1 x2 w d native tac_t tac_b :=
  match goal with
  | |- plot2 ?f ?ox ?dx ?oy' ?dy' (Zpos ?h) ?p2 ⇒
    let p1 := fresh "__p1" in
    evar (p1 : list I.type) ;
    let Hp := fresh "__Hp" in
    assert (Hp: plot1 f ox dx p1) by plot1_aux1 prec x1 x2 w h d native tac_t ;
    revert Hp ;
    let y1y2 := tac_b prec in
    let oy := constr:(fst y1y2) in
    let dy := eval vm_compute in (F.div_UP prec (F.sub_UP prec (snd y1y2) oy) (F.fromZ_DN prec (Zpos h))) in
    let oyr := eval cbv -[IZR Rdiv] in (proj_val (I.F.convert oy)) in
    let dyr := eval cbv -[IZR Rdiv] in (proj_val (I.F.convert dy)) in
    unify oy' oyr ;
    unify dy' dyr ;
    refine (clamp_plot_correct prec _ _ _ _ _ oyr dyr _ _ _ _ (I.singleton_correct oy) (J.inv_correct prec _ _ (I.singleton_correct dy)) _) ;
    [ try apply IZR_lt ;
      apply Rdiv_lt_0_compat ;
      now apply IZR_lt
    | unify_eq false ]
  end.

Definition get_threshold prec hyps pf cf ox dx w :=
  let w' := 50%Z in
  let dx := I.mul prec (I.singleton dx) (I.div prec (I.fromZ prec (Zpos w)) (I.fromZ prec w')) in
  let yi := bound_plot prec hyps pf cf (I.singleton ox) dx w' in
  (I.lower yi, I.upper yi).

Ltac plot_get_threshold prec ox dx w :=
  match goal with
  | |- eval_hyps ?hyps _ (plot1 (fun t ⇒ eval_real' ?pf (t :: _) ?cf) _ _ _) ⇒
    eval vm_compute in (get_threshold prec hyps pf cf ox dx w)
  end.

Ltac plot_get_bounds prec :=
  match goal with
  | |- plot1 _ _ _ ?p → _ ⇒
    let p := eval vm_compute in p in
    eval vm_compute in (get_bounds prec p)
  end.

Ltac plot_y_get_threshold y1 y2 prec ox dx w := constr:((y1, y2)).

Ltac plot_y_get_bounds y1 y2 prec := constr:((y1, y2)).

Ltac do_plot f x1 x2 prec degree width height native :=
  let p := fresh "__p2" in
  evar (p : list (Z × Z)) ;
  refine (_: plot2 f _ _ _ _ (Zpos height) p) ;
  plot2_aux prec x1 x2 width degree native plot_get_threshold plot_get_bounds.

Ltac do_plot_y f x1 x2 y1 y2 prec degree width height native :=
  let p := fresh "__p2" in
  let y1 := reify y1 constr:(@nil R) in
  let y2 := reify y2 constr:(@nil R) in
  let y1 := eval vm_compute in (I.lower (T.eval_bnd prec y1)) in
  let y2 := eval vm_compute in (I.upper (T.eval_bnd prec y2)) in
  evar (p : list (Z × Z)) ;
  refine (_: plot2 f _ _ _ _ (Zpos height) p) ;
  plot2_aux prec x1 x2 width degree native ltac:(plot_y_get_threshold y1 y2) ltac:(plot_y_get_bounds y1 y2).

End PlotTacticAux.