(*
  Formalization of the base expression language for the daisy framework
*)
needs "Infra/tactics.hl";;
needs "Infra/RealConstruction.hl";;
(*
  Expressions will use binary operators.
  Define them first
*)
let binop_IND, binop_REC = define_type
  "binop = Plus | Sub | Mult | Div ";;
(*
  Define an evaluation function for binary operators.
  Errors are added on the expression evaluation level later
*)
let eval_binop = new_recursive_definition binop_REC
  `(eval_binop Plus v1 v2 = v1 + v2) /\
  (eval_binop Sub v1 v2 = v1 - v2) /\
  (eval_binop Mult v1 v2 = v1 * v2) /\
  (eval_binop Div v1 v2 = v1 / v2)`;;
(*
  Define expressions parametric over some value type V.
  Will ease reasoning about different instantiations later.
*)
let exp_IND, exp_REC = define_type
  "exp = Var num
  |Param num
  | Const V
  | Binop binop exp exp";;
(*
  Define the machine epsilon for floating point operations.
  Current value: 2^(-53)
*)
let machineEpsilon =  define `machineEpsilon:real = realFromNum 1 1 53`;;
(*
  Define a perturbation function to ease writing of basic definitions
*)
let perturb = define `(perturb:real->real->real) (r:real) (e:real) = r * ((&1) + e)`;;
(*
  Define expression evaluation relation parametric by an "error" delta.
  This value will be used later to express float computations using a perturbation
  of the real valued computation by (1 + d)
*)
let eval_exp_RULES, eval_exp_IND, eval_exp_CASES = new_inductive_definition
  `(!eps env v.
      eval_exp eps env (Var v) (env v)) /\
  (!eps env v delta.
     abs delta <= eps ==>
     eval_exp eps env (Param v) (perturb (env v) delta)) /\
  (!eps env n delta.
     abs delta <= eps ==>
     eval_exp eps env (Const n) (perturb n delta)) /\
  (!eps env b e1 e2 v1 v2 delta.
     eval_exp eps env e1 v1 /\
     eval_exp eps env e2 v2 /\
     abs delta <= eps ==>
     eval_exp eps env (Binop b e1 e2) (perturb (eval_binop b v1 v2) delta))`;;

(* Inversion Lemmata, to do proofs similar to Coq proofs *)
let var_inv =
  prove(
    `!eps env v val.
       eval_exp eps env (Var v) val ==> val = env v`,
     ONCE_REWRITE_TAC
       [eval_exp_CASES] THEN
       INTRO_TAC "!eps env v val; cases_exp" THEN
       REMOVE_THEN "cases_exp" (DESTRUCT_TAC "varc | paramc | constc | binopc")
       THENL [
         REMOVE_THEN "varc" (DESTRUCT_TAC "@v2. var_eq val_eq") THEN
           RULE_ASSUM_TAC (REWRITE_RULE [injectivity "exp"]) THEN ASM_REWRITE_TAC [];
         REMOVE_THEN "paramc" (DESTRUCT_TAC "@n. @d. var_eq _ _") THEN lcontra "var_eq" "exp";
         REMOVE_THEN "constc" (DESTRUCT_TAC "@n. @d. var_eq _ _") THEN
           lcontra "var_eq" "exp";
        REMOVE_THEN "binopc"
          (DESTRUCT_TAC "@b. @e1. @e2. @v1. @v2. @d. var_eq _ _ _ _") THEN
          lcontra "var_eq" "exp"] THEN
       FIRST_ASSUM CONTR_TAC);;

let var_equiv =
  prove (
    `!eps env v val. eval_exp eps env (Var v) val <=> val = env v`,
    intros "!eps env v val"
      THEN EQ_TAC
      THENL[
        intros "eval_exp_var"
          THEN MATCH_MP_TAC var_inv
          THEN EXISTS_TAC `eps:real`
          THEN ASM_REWRITE_TAC[];
        intros "val_eq"
          THEN ASM_REWRITE_TAC[eval_exp_RULES]]);;

let param_inv =
  prove
 (`!eps env n v.
  eval_exp eps env (Param n) v ==> ?d. v = perturb (env n) d /\ abs d <= eps`,
  ONCE_REWRITE_TAC
    [eval_exp_CASES] THEN
  INTRO_TAC "!eps env n v; cases_exp" THEN
  REMOVE_THEN "cases_exp" (DESTRUCT_TAC "varc | paramc | constc | binopc") THENL
   [REMOVE_THEN "varc" (DESTRUCT_TAC "@v. const_eq _") THEN
      lcontra "const_eq" "exp";
    REMOVE_THEN "paramc" (DESTRUCT_TAC "@n2. @d. c_eq v_eq abs_less") THEN
    EXISTS_TAC `d:real` THEN
    RULE_ASSUM_TAC (REWRITE_RULE [injectivity "exp"]) THEN ASM_REWRITE_TAC [];
    REMOVE_THEN "constc" (DESTRUCT_TAC "@n. @d. var_eq _ _") THEN lcontra "var_eq" "exp";
    REMOVE_THEN "binopc"
      (DESTRUCT_TAC "@b. @e1. @e2. @v1. @v2. @d. const_eq _ _ _ _") THEN
      lcontra "const_eq" "exp"]);;

let const_inv =
  prove
 (`!eps env n v.
  eval_exp eps env (Const n) v ==> ?d. v = perturb n d /\ abs d <= eps`,
  ONCE_REWRITE_TAC
    [eval_exp_CASES] THEN
  INTRO_TAC "!eps env n v; cases_exp" THEN
  REMOVE_THEN "cases_exp" (DESTRUCT_TAC "varc | paramc | constc | binopc") THENL
   [REMOVE_THEN "varc" (DESTRUCT_TAC "@v. const_eq _") THEN
   lcontra "const_eq" "exp";
    REMOVE_THEN "paramc" (DESTRUCT_TAC "@n. @d. var_eq _ _") THEN lcontra "var_eq" "exp";
    REMOVE_THEN "constc" (DESTRUCT_TAC "@n2. @d. c_eq v_eq abs_less") THEN
    EXISTS_TAC `d:real` THEN
    RULE_ASSUM_TAC (REWRITE_RULE [injectivity "exp"]) THEN ASM_REWRITE_TAC [];
    REMOVE_THEN "binopc"
      (DESTRUCT_TAC "@b. @e1. @e2. @v1. @v2. @d. const_eq _ _ _ _") THEN
   lcontra "const_eq" "exp"]);;

let binop_inv =
  prove
    ( `!eps env b e1 e2 v.
        eval_exp eps env (Binop b e1 e2) v ==>
        ?d v1 v2.
        v = perturb (eval_binop b v1 v2) d /\
            eval_exp eps env e1 v1 /\
            eval_exp eps env e2 v2 /\
            abs d <= eps`,
      INTRO_TAC "!eps env b e1 e2 v; cases_exp"
        THEN RULE_ASSUM_TAC (ONCE_REWRITE_RULE[eval_exp_CASES])
        THEN destruct "cases_exp" "varc | paramc | constc | binopc "
        THENL [
          destruct "varc" "@v. b_eq _" THEN lcontra "b_eq" "exp";
          destruct "paramc" "@n. @d. b_eq _ _" THEN lcontra "b_eq" "exp";
          destruct "constc" "@n. @d. b_eq _ _" THEN lcontra "b_eq" "exp";
          destruct "binopc" "@b2. @e1'. @e2'. @v1. @v2. @d. b_eq  v_eq eval_e1 eval_e2 abs_delta"
            THEN EXISTS_TAC `d:real` THEN EXISTS_TAC `v1:real` THEN EXISTS_TAC `v2:real`
            THEN RULE_ASSUM_TAC (REWRITE_RULE[injectivity "exp"]) THEN ASM_REWRITE_TAC[]]);;

let abs_leq_0_impl_zero =
  prove (
    `!d:real. abs d <= &0 ==> d = &0`,
    INTRO_TAC "!d; abs_leq_0"
      THEN SUBGOAL_TAC "abs_geq_0" `&0 <= abs d` [(REWRITE_TAC[REAL_ABS_POS])]
      THEN SUBGOAL_TAC "abs_geq_leq_0" `abs d <= &0 /\ &0 <= abs d` [CONJ_TAC THEN (ASM_REWRITE_TAC[])]
      THEN SUBGOAL_TAC "abs_eq_0" `abs d = &0` [ALL_TAC]
          THEN MP_TAC (ASSUME `abs d <= &0 /\ &0 <= abs d`)
          THEN ASM_REWRITE_TAC [SPECL [`abs d`; `&0:real`] REAL_LE_ANTISYM]
          THEN MP_TAC (ASSUME `abs d = &0`)
          THEN ASM_REWRITE_TAC [REAL_ABS_ZERO]);;


let perturb_0_val =
  prove (
    `!(v:real) (delta:real). abs delta <= &0 ==> perturb v delta = v`,
    intros "!v delta; abs_leq_0"
      THEN SIMP_TAC[perturb]
      THEN RULE_ASSUM_TAC (MATCH_MP abs_leq_0_impl_zero)
      THEN ASM_REWRITE_TAC[]
      THEN REAL_ARITH_TAC);;

g  (`!e:(real)exp v1:real v2:real env.
      eval_exp (&0) env e v1 /\ eval_exp (&0) env e v2 ==> v1 = v2`);;
e (MATCH_MP_TAC exp_IND THEN STRIP_TAC);;
(* Var Case *)
e (INTRO_TAC "!a v1 v2 env; eval_exps");;
e (destruct "eval_exps" "eval_v1 eval_v2");;
e (SUBGOAL_TAC "v1_eq_env_a" `v1:real = (env (a:num)):real`
    [MATCH_MP_TAC var_inv THEN EXISTS_TAC `&0` THEN ASM_SIMP_TAC[]]);;
e (SUBGOAL_TAC "v2_eq_env_a" `v2:real = (env (a:num)):real`
    [MATCH_MP_TAC var_inv THEN EXISTS_TAC `&0` THEN ASM_SIMP_TAC[]]);;
e (ASM_REWRITE_TAC[]);;
(* Param Case *)
e (STRIP_TAC);;
e (INTRO_TAC "!a v1 v2 env; eval_exps");;
e (destruct "eval_exps" "eval_v1 eval_v2");;
e (SUBGOAL_TAC "v1_eq_a" `?d. v1:real = perturb ((env (a:num)):real) d /\ abs d <= &0`
    [MATCH_MP_TAC param_inv THEN ASM_SIMP_TAC[]]);;
e (SUBGOAL_TAC "v2_eq_a" `?d. v2:real = perturb ((env (a:num)):real) d /\ abs d <= &0`
    [MATCH_MP_TAC param_inv THEN ASM_SIMP_TAC[]]);;
e (destruct "v1_eq_a" "@d1. v1_eq abs");;
e (destruct "v2_eq_a" "@d2. v2_eq abs2");;
e (USE_THEN "abs" (fun th -> ASM_REWRITE_TAC[MP (SPECL [`((env:num->real) a):real`; `d1:real`] perturb_0_val) th]));;
e (USE_THEN "abs2" (fun th -> ASM_REWRITE_TAC[MP (SPECL [`((env:num->real) a):real`; `d2:real`] perturb_0_val) th]));;
(* Const Case *)
e (STRIP_TAC);;
e (INTRO_TAC "!a v1 v2 env; eval_exps");;
e (destruct "eval_exps" "eval_v1 eval_v2");;
e (SUBGOAL_TAC "v1_eq_a" `?d. v1:real = perturb (a:real) d /\ abs d <= &0`
    [MATCH_MP_TAC const_inv THEN EXISTS_TAC `env:num->real` THEN ASM_SIMP_TAC[]]);;
e (SUBGOAL_TAC "v2_eq_a" `?d. v2:real = perturb (a:real) d /\ abs d <= &0`
    [MATCH_MP_TAC const_inv THEN EXISTS_TAC `env:num->real` THEN ASM_SIMP_TAC[]]);;
e (destruct "v1_eq_a" "@d1. v1_eq abs");;
e (destruct "v2_eq_a" "@d2. v2_eq abs2");;
e (SUBGOAL_TAC "abs_d1_0" `d1 = &0` [MATCH_MP_TAC abs_leq_0_impl_zero THEN ASM_SIMP_TAC[]]);;
e (SUBGOAL_TAC "abs_d2_0" `d2 = &0` [MATCH_MP_TAC abs_leq_0_impl_zero THEN ASM_SIMP_TAC[]]);;
e (ASM_REWRITE_TAC[]);;
(* Binop Case *)
e (REPEAT STRIP_TAC);;
e (FIND_ASSUM (LABEL_TAC "IH1") `!(v1:real) (v2:real) (env:num->real).
          eval_exp (&0) (env:num->real) (a1:(real)exp) (v1:real) /\
          eval_exp (&0) (env:num->real) (a1:(real)exp) (v2:real)
   ==> (v1:real) = (v2:real)`);;
e (FIND_ASSUM (LABEL_TAC "IH2") `!(v1:real) (v2:real) (env:num->real).
          eval_exp (&0) (env:num->real) (a2:(real)exp) (v1:real) /\
          eval_exp (&0) (env:num->real) (a2:(real)exp) (v2:real)
   ==> (v1:real) = (v2:real)`);;
e (SUBGOAL_TAC "v1_eq_op_a" `?d v1' v2'. v1:real = perturb (eval_binop (a0:binop) (v1':real) (v2':real)) d /\ eval_exp (&0) env a1 v1' /\ eval_exp (&0) env a2 v2' /\ abs d <= &0`
    [MATCH_MP_TAC binop_inv THEN ASM_SIMP_TAC[]]);;
e (SUBGOAL_TAC "v2_eq_op_a" `?d v1' v2'. v2:real = perturb (eval_binop (a0:binop) (v1':real) (v2':real)) d /\ eval_exp (&0) env a1 v1' /\ eval_exp (&0) env a2 v2' /\ abs d <= &0`
    [MATCH_MP_TAC binop_inv THEN ASM_SIMP_TAC[]]);;
e (destruct "v1_eq_op_a" "@d1. @v11. @v21. v1_eq eval_a11 eval_a21 abs_0");;
e (destruct "v2_eq_op_a" "@d2. @v12. @v22. v2_eq eval_a12 eval_a22 abs_0");;
e (SUBGOAL_TAC "abs_d1_0" `d1 = &0` [MATCH_MP_TAC abs_leq_0_impl_zero THEN ASM_SIMP_TAC[]]);;
e (SUBGOAL_TAC "abs_d2_0" `d2 = &0` [MATCH_MP_TAC abs_leq_0_impl_zero THEN ASM_SIMP_TAC[]]);;
e (SUBGOAL_TAC "v11_eq_v12" `v11:real = v12:real` [REMOVE_THEN "IH1" (MATCH_MP_TAC)]);;
e (EXISTS_TAC `env:num->real` THEN ASM_SIMP_TAC[]);;
e (SUBGOAL_TAC "v21_eq_v22" `v21:real = v22:real` [REMOVE_THEN "IH2" (MATCH_MP_TAC)]);;
e (EXISTS_TAC `env:num->real` THEN ASM_SIMP_TAC[]);;
e (ASM_REWRITE_TAC[]);;
let eval_0_det = top_thm();;

(* Alias for proof similarity *)
let Rsub_eq_Ropp_Rplus = real_sub;;

(** weird: cannot name lemma Rabs_err_simpl **)
let abs_err_simpl  =
  prove (
    `!(a:real) (b:real).
     abs (a - (a * (&1 + b))) = abs (a * b)`,
    intros "!a b"
      THEN REWRITE_TAC [REAL_ADD_LDISTRIB; REAL_MUL_RID]
      THEN SUBGOAL_TAC "arith_simp" `(a:real) - ((a:real) + (a:real) * (b:real)) = -- (a * b)` [REAL_ARITH_TAC]
      THEN REMOVE_THEN "arith_simp" (fun th -> REWRITE_TAC [th])
      THEN REWRITE_TAC [REAL_ABS_NEG]);;

(**
  TODO: Check wether we need Rabs (n * machineEpsilon) instead
**)

let const_abs_err_bounded =
  prove (
    `!(n:real) (nR:real) (nF:real) (cenv:num ->real).
      eval_exp (&0) cenv (Const n) nR /\
      eval_exp machineEpsilon cenv (Const n) nF ==>
      abs (nR - nF) <= abs n * machineEpsilon`,
    intros "!n nR nF cenv; eval_real eval_float"
      THEN RULE_ASSUM_TAC (MATCH_MP const_inv)
      THEN destruct "eval_real" "@d1. nR_eq abs_d1_0"
      THEN destruct "eval_float" "@d2. nF_eq abs_d2_leq"
      THEN ASM_SIMP_TAC []
      THEN USE_THEN "abs_d1_0"
      (fun th ->
         ASM_REWRITE_TAC [MP (SPECL [`(n:real)`; `d1:real`] perturb_0_val) th])
      THEN REWRITE_TAC [perturb;abs_err_simpl; REAL_ABS_MUL]
      THEN MATCH_MP_TAC REAL_LE_LMUL
      THEN ASM_REWRITE_TAC [REAL_ABS_POS]);;

(**
  TODO: Maybe improve bound by adding interval constraint and proving that it is leq than maxAbs of bounds
  (nlo <= cenv n <= nhi)%R ->   (Rabs (nR - nF) <= Rabs ((Rmax (Rabs nlo) (Rabs nhi)) * machineEpsilon))%R.
**)
let param_abs_err_bounded =
  prove (
    `!(n:num) (nR:real) (nF:real) (cenv:num->real).
      eval_exp (&0) cenv (Param n) nR /\
      eval_exp machineEpsilon cenv (Param n) nF ==>
      (abs (nR - nF) <= abs (cenv n) * machineEpsilon)`,
    intros "!n nR nF cenv; eval_real eval_float"
      THEN RULE_ASSUM_TAC (MATCH_MP param_inv)
      THEN destruct "eval_real" "@d1. nR_eq abs_d1_0"
      THEN destruct "eval_float" "@d2. nF_eq abs_d1_leq"
      THEN ASM_SIMP_TAC[]
      THEN USE_THEN "abs_d1_0"
      (fun th ->
         ASM_REWRITE_TAC [MP (SPECL [`(cenv n:real)`; `d1:real`] perturb_0_val) th])
      THEN REWRITE_TAC [perturb; abs_err_simpl; REAL_ABS_MUL]
      THEN MATCH_MP_TAC REAL_LE_LMUL THEN ASM_REWRITE_TAC[REAL_ABS_POS]);;
(*
  Using the parametric expressions, define boolean expressions for conditionals
*)
let bexp_INDUCT, bexp_REC = define_type
  "bexp = leq (V)exp (V)exp
        | less (V)exp (V)exp";;

(*
  Define evaluation of boolean expressions
*)
let bval_exp_RULES, bval_exp_IND, bval_exp_CASES = new_inductive_definition
  `(!eps env e1 e2 v1 v2.
      eval_exp eps env e1 v1 /\
      eval_exp eps env e2 v2 ==>
      bval eps env (leq e1 e2) (v1 <= v2))/\
  (!eps env e1 e2 v1 v2.
      eval_exp eps env e1 v1 /\
      eval_exp eps env e2 v2 ==>
      bval eps env (less e1 e2) (v1 < v2))`;;

(* Simplify arithmetic later by making > >= only abbreviations *)
let gr_simps  = define
  `gr = \(e1:(V)exp). \(e2:(V)exp). less e2 e1`;;
let greq_simps = define
  `greq = \(e1:(V)exp). \(e2:(V)exp). leq e2 e1`;;