Theory Abstract_Rigorous_Numerics

theory Abstract_Rigorous_Numerics
imports
  "HOL-Library.Parallel"
  Transfer_Euclidean_Space_Vector
  "../Refinement/Enclosure_Operations"
  "../Refinement/Refine_Vector_List"
  "../Refinement/Refine_Hyperplane"
  "../Refinement/Refine_Interval"
  "../Refinement/Refine_Invar"
  "../Refinement/Refine_Unions"
  "../Refinement/Refine_Info"
begin

section ‹misc›

lemma length_listset: "xi ∈ listset xs ⟹ length xi = length xs"
  by (induction xs arbitrary: xi) (auto simp: set_Cons_def)

lemma Nil_mem_listset[simp]: "[] ∈ listset list ⟷ list = []"
  by (induction list) (auto simp: set_Cons_def)

lemma sing_mem_listset_iff[simp]: "[b] ∈ listset ys ⟷ (∃z. ys = [z] ∧ b ∈ z)"
  ― ‹TODO: generalize to Cons?›
  by (cases ys) (auto simp: set_Cons_def)


unbundle (in autoref_syn) no funcset_syntax

definition cfuncset where "cfuncset l u X = funcset {l .. u} X ∩ Collect (continuous_on {l .. u})"
lemma cfuncset_iff: "f ∈ cfuncset l u X ⟷ (∀i∈{l .. u}. f i ∈ X) ∧ continuous_on {l .. u} f"
  unfolding cfuncset_def by auto

lemma cfuncset_continuous_onD: "f ∈ cfuncset 0 h X ⟹ continuous_on {0..h} f"
  by (simp add: cfuncset_iff)


section ‹Implementations›

subsection ‹locale for sets›

definition "product_listset xs ys = (λ(x, y). x @ y) ` ((xs::real list set) × (ys::real list set))"

abbreviation "rl_rel ≡ ⟨rnv_rel⟩list_rel"

abbreviation "slp_rel ≡ ⟨Id::floatarith rel⟩list_rel"
abbreviation "fas_rel ≡ ⟨Id::floatarith rel⟩list_rel"

type_synonym 'b reduce_argument = "'b list ⇒ nat ⇒ real list ⇒ bool" ― ‹is this too special?›

record 'b numeric_options =
  precision :: nat
  adaptive_atol :: real
  adaptive_rtol :: real
  method_id :: nat
  start_stepsize :: real
  iterations :: nat
  halve_stepsizes :: nat
  widening_mod :: nat
  rk2_param :: real
  default_reduce :: "'b reduce_argument"
  printing_fun :: "bool ⇒ 'b list ⇒ unit"
  tracing_fun :: "string ⇒ 'b list option ⇒ unit"

record 'b reach_options =
  max_tdev_thres :: "real"
  pre_split_reduce :: "'b reduce_argument"
  pre_inter_granularity :: "real"
  post_inter_granularity :: "real"
  pre_collect_granularity :: real
  max_intersection_step :: real

definition "reach_options_rel TYPE('b) = (UNIV::('b reach_options × unit) set)"
lemma sv_reach_options_rel[relator_props]: "single_valued (reach_options_rel TYPE('a))"
  by (auto simp: reach_options_rel_def single_valued_def)

definition "reduce_argument_rel TYPE('b) = (UNIV::('b reduce_argument × unit) set)"
lemma sv_reduce_argument_rel[relator_props]: "single_valued (reduce_argument_rel TYPE('a))"
  by (auto simp: reduce_argument_rel_def single_valued_def)


definition [refine_vcg_def, simp]: "max_tdev_thres_spec (ro::unit) = SPEC (λx::real. True)"
definition [refine_vcg_def, simp]: "max_intersection_step_spec (ro::unit) = SPEC (λx::real. True)"
definition [refine_vcg_def, simp]: "pre_inter_granularity_spec (ro::unit) = SPEC (λx::real. True)"
definition [refine_vcg_def, simp]: "post_inter_granularity_spec (ro::unit) = SPEC (λx::real. True)"
definition [refine_vcg_def, simp]: "pre_collect_granularity_spec (ro::unit) = SPEC (λx::real. True)"

lemma reach_optns_autoref[autoref_rules]:
  includes autoref_syntax
  shows
    "(λx. RETURN (max_tdev_thres x), max_tdev_thres_spec) ∈ reach_options_rel TYPE('b) → ⟨rnv_rel⟩nres_rel"
    "(λx. RETURN (pre_inter_granularity x), pre_inter_granularity_spec) ∈ reach_options_rel TYPE('b)  → ⟨rnv_rel⟩nres_rel"
    "(λx. RETURN (post_inter_granularity x), post_inter_granularity_spec) ∈ reach_options_rel TYPE('b) → ⟨rnv_rel⟩nres_rel"
    "(λx. RETURN (pre_collect_granularity x), pre_collect_granularity_spec) ∈ reach_options_rel TYPE('b)  → ⟨rnv_rel⟩nres_rel"
    "(λx. RETURN (max_intersection_step x), max_intersection_step_spec) ∈ reach_options_rel TYPE('b)  → ⟨rnv_rel⟩nres_rel"
  by (auto simp: nres_rel_def)


record 'b approximate_set_ops =
  appr_of_ivl::"real list ⇒ real list ⇒ 'b list"
  product_appr::"'b list ⇒ 'b list ⇒ 'b list"
  msum_appr::"'b list ⇒ 'b list ⇒ 'b list"
  inf_of_appr::"'b numeric_options ⇒ 'b list ⇒ real list"
  sup_of_appr::"'b numeric_options ⇒ 'b list ⇒ real list"
  split_appr::"nat ⇒ 'b list ⇒ 'b list × 'b list"
  appr_inf_inner::"'b numeric_options ⇒ 'b list ⇒ real list ⇒ real"
  appr_sup_inner::"'b numeric_options ⇒ 'b list ⇒ real list ⇒ real"
  inter_appr_plane::"'b numeric_options ⇒ 'b list ⇒ real list sctn ⇒ 'b list dres"
  reduce_appr::"'b numeric_options ⇒ 'b reduce_argument ⇒ 'b list ⇒ 'b list"
  width_appr::"'b numeric_options ⇒ 'b list ⇒ real"
  approx_slp_dres::"'b numeric_options ⇒ nat ⇒ slp ⇒ 'b list ⇒ 'b list option dres"
  approx_euclarithform::"'b numeric_options ⇒ form ⇒ 'b list ⇒ bool dres"
  approx_isFDERIV::"'b numeric_options ⇒ nat ⇒ nat list ⇒ floatarith list ⇒ 'b list ⇒ bool dres"

primrec concat_appr where
  "concat_appr ops [] = []"
| "concat_appr ops (x#xs) = product_appr ops x (concat_appr ops xs)"


unbundle autoref_syntax

locale approximate_sets =
  fixes ops :: "'b approximate_set_ops"
    and set_of_appr::"'b list ⇒ real list set"
    and appr_rell :: "('b list × real list set) set"
    and optns :: "'b numeric_options"
  assumes appr_rell_internal: "appr_rell ≡ br set_of_appr top"
  assumes transfer_operations_rl:
    "SIDE_PRECOND (list_all2 (≤) xrs yrs) ⟹
      (xri, xrs) ∈ ⟨rnv_rel⟩list_rel ⟹
      (yri, yrs) ∈ ⟨rnv_rel⟩list_rel ⟹
      (appr_of_ivl ops xri yri, lv_ivl $ xrs $ yrs) ∈ appr_rell"
    "(product_appr ops, product_listset) ∈ appr_rell → appr_rell → appr_rell"
    "(msum_appr ops, (λxs ys. {List.map2 (+) x y |x y. x ∈ xs ∧ y ∈ ys})) ∈ appr_rell → appr_rell → appr_rell"
    "(xi, x) ∈ appr_rell ⟹ length xi = d ⟹
      (RETURN (inf_of_appr ops optns xi), Inf_specs d x) ∈ ⟨rl_rel⟩nres_rel"
    "(xi, x) ∈ appr_rell ⟹ length xi = d ⟹
      (RETURN (sup_of_appr ops optns xi), Sup_specs d x) ∈ ⟨rl_rel⟩nres_rel"
    "(ni, n) ∈ nat_rel ⟹ (xi, x) ∈ appr_rell ⟹ length xi = d ⟹
      (RETURN (split_appr ops ni xi), split_spec_params d n x) ∈ ⟨appr_rell ×⇩r appr_rell⟩nres_rel"
    "(RETURN o2 appr_inf_inner ops optns, Inf_inners) ∈ appr_rell → rl_rel → ⟨rnv_rel⟩nres_rel"
    "(RETURN o2 appr_sup_inner ops optns, Sup_inners) ∈ appr_rell → rl_rel → ⟨rnv_rel⟩nres_rel"
    "(xi, x) ∈ appr_rell ⟹ length xi = d ⟹ length (normal si) = d ⟹ d > 0 ⟹ (si, s) ∈ ⟨rl_rel⟩sctn_rel ⟹
      (nres_of (inter_appr_plane ops optns xi si), inter_sctn_specs d x s) ∈ ⟨appr_rell⟩nres_rel"
    "(xi, x) ∈ appr_rell ⟹ length xi = d ⟹ (rai, ra) ∈ reduce_argument_rel TYPE('b) ⟹
      (RETURN (reduce_appr ops optns rai xi), reduce_specs d ra x) ∈ ⟨appr_rell⟩nres_rel"
    "(RETURN o width_appr ops optns, width_spec) ∈ appr_rell → ⟨rnv_rel⟩nres_rel"
    "(nres_of o3 approx_slp_dres ops optns, approx_slp_spec fas) ∈ nat_rel → slp_rel → appr_rell → ⟨⟨appr_rell⟩option_rel⟩nres_rel"
assumes approx_euclarithform[unfolded comps, autoref_rules]:
  "(nres_of o2 approx_euclarithform ops optns, approx_form_spec) ∈ Id → appr_rell → ⟨bool_rel⟩nres_rel"
assumes approx_isFDERIV[unfolded comps, autoref_rules]:
  "(λN xs fas vs. nres_of (approx_isFDERIV ops optns N xs fas vs), isFDERIV_spec) ∈
  nat_rel → ⟨nat_rel⟩list_rel → ⟨Id⟩list_rel →  appr_rell → ⟨bool_rel⟩nres_rel"
assumes set_of_appr_nonempty[simp]: "set_of_appr X ≠ {}"
assumes length_set_of_appr: "xrs ∈ set_of_appr xs ⟹ length xrs = length xs"
assumes set_of_appr_project: "xrs ∈ set_of_appr xs ⟹ (⋀i. i ∈ set is ⟹ i < length xs) ⟹
    map ((!) xrs) is ∈ set_of_appr (map ((!) xs) is)"
assumes set_of_apprs_ex_Cons: "xrs ∈ set_of_appr xs ⟹ ∃r. r#xrs ∈ set_of_appr (b#xs)"
assumes set_of_apprs_Nil[simp]: "xrs ∈ set_of_appr [] ⟷ xrs = []"
begin

abbreviation "reach_optns_rel ≡ reach_options_rel TYPE('b)"

definition "appr_rel = appr_rell O ⟨lv_rel⟩set_rel"
lemmas [autoref_rel_intf] = REL_INTFI[of appr_rel i_appr]

definition [simp]: "op_concat_listset xs = concat ` listset xs"

lemma [autoref_op_pat_def]: "concat ` listset xs ≡ OP op_concat_listset $ xs"
  by simp

lemma list_all2_replicate [simp]: "list_all2 (≤) xs xs" for xs::"'a::order list"
  by (auto simp: list_all2_iff in_set_zip)

lemma length_appr_of_ivl[simp]:
  "length (appr_of_ivl ops xs ys) = length xs"
  if "list_all2 (≤) xs ys"
  using transfer_operations_rl(1)[of xs ys xs ys] that
    apply (simp add: appr_rell_internal br_def lv_ivl_def)
  apply (auto simp: appr_rell_internal br_def list_all2_lengthD dest!: length_set_of_appr)
  using length_set_of_appr by fastforce

definition [simp]: "op_atLeastAtMost_appr = atLeastAtMost"
lemma [autoref_op_pat]: "atLeastAtMost ≡ OP op_atLeastAtMost_appr"
  by simp

definition card_info::"_ set ⇒ nat nres" where [refine_vcg_def]: "card_info x = SPEC top" ― ‹‹op_set_wcard››

definition [simp]: "set_of_coll X = X"

definition [simp]: "ivl_to_set X = X"

definition "ivls_of_sets X = do {
    XS ← (sets_of_coll (X:::clw_rel appr_rel) ::: ⟨⟨appr_rel⟩list_wset_rel⟩nres_rel);
    FORWEAK XS (RETURN (op_empty_coll:::clw_rel lvivl_rel))
      (λX. do {
        (i, s) ← op_ivl_rep_of_set X;
        RETURN (mk_coll (op_atLeastAtMost_ivl i s:::⟨lv_rel⟩ivl_rel))
      })
      (λa b. RETURN (a ∪ b))
  }"
sublocale autoref_op_pat_def ivls_of_sets .

definition "ivl_rep_of_set_coll X = do { Xs ← sets_of_coll (X:::clw_rel appr_rel); op_ivl_rep_of_sets Xs}"
sublocale autoref_op_pat_def ivl_rep_of_set_coll .

definition [simp]: "sets_of_ivls X = X"
sublocale autoref_op_pat_def sets_of_ivls .

definition "op_intersects X sctn = (do {
    ii ← Inf_inner X (normal sctn);
    si ← Sup_inner X (normal sctn);
    RETURN (ii ≤ pstn sctn ∧ si ≥ pstn sctn)
  })"
sublocale autoref_op_pat_def op_intersects .

definition "sbelow_sctns_coll XS sctns = do {
    XS ← (sets_of_coll XS ::: ⟨⟨appr_rel⟩list_wset_rel⟩nres_rel);
    FORWEAK XS (RETURN True) (λX. sbelow_sctns X sctns) (λa b. RETURN (a ∧ b))
  }"
sublocale autoref_op_pat_def sbelow_sctns_coll .

definition "below_sctn_coll XS sctn = do {
    XS ← (sets_of_coll XS ::: ⟨⟨appr_rel⟩list_wset_rel⟩nres_rel);
    FORWEAK XS (RETURN True) (λX. below_sctn X sctn) (λa b. RETURN (a ∧ b))
  }"
sublocale autoref_op_pat_def below_sctn_coll .


definition "intersects_clw X sctn = (do {
    XS ← sets_of_coll (X:::clw_rel appr_rel);
    FORWEAK XS (RETURN False) (λX. op_intersects X sctn) (λa b. RETURN (a ∨ b))
  })"
sublocale autoref_op_pat_def intersects_clw .

definition "op_subset X ivl = do {
    (i', s') ← ((ivl_rep ((ivl))));
    (i, s) ← (op_ivl_rep_of_set ((X::'a::executable_euclidean_space set)));
    RETURN (((i' ≤ i):::bool_rel) ∧ ((s ≤ s'):::bool_rel))
  }"
sublocale autoref_op_pat_def op_subset .

definition [simp]: "subset_spec_coll X ivl = do {
    XS ← (sets_of_coll X:::⟨⟨appr_rel⟩list_wset_rel⟩nres_rel);
    FORWEAK XS (RETURN True)
      (λX. op_subset X ivl)
      (λx y. RETURN (x ∧ y))
  }"
sublocale autoref_op_pat_def subset_spec_coll .

definition "op_project_set X b y = inter_sctn_spec X (Sctn b y)"
sublocale autoref_op_pat_def op_project_set .

definition [simp]: "project_set_clw X b y = do {
    XS ← (sets_of_coll (X:::clw_rel appr_rel));
    FORWEAK XS (RETURN op_empty_coll) (λX. do {
      P ← op_project_set X b y;
      RETURN (mk_coll P)
    }) (λX Y. RETURN (Y ∪ X))
  }"
sublocale autoref_op_pat_def project_set_clw .

definition "subset_spec_ivls X Y = do {
    Ys ← sets_of_coll Y; FORWEAK Ys (RETURN False) (λY. op_subset X Y) (λa b. RETURN (a ∨ b))
  }"
sublocale autoref_op_pat_def subset_spec_ivls .

definition "subset_spec_ivls_clw M X Y = do {
    X ← split_along_ivls M X Y;
    X ← sets_of_coll (sets_of_ivls X);
    FORWEAK X (RETURN True) (λX. subset_spec_ivls X Y) (λa b. RETURN (a ∧ b))
  }"
sublocale autoref_op_pat_def subset_spec_ivls_clw .

definition [simp]: "REMDUPS x = x"
sublocale autoref_op_pat_def REMDUPS .

definition "split_along_ivls2 M X Y = do {
    Xs ← sets_of_coll X;
    Rs ←FORWEAK Xs (RETURN op_empty_coll) (λX. do {
      (I, N) ← split_intersecting Y (mk_coll X);
      split_along_ivls M (mk_coll X) I
    }) (λx y. RETURN (y ∪ x));
    RETURN (REMDUPS Rs)
  }"
sublocale autoref_op_pat_def split_along_ivls2 .

definition [simp]: "op_list_of_eucl_image X = list_of_eucl ` X"
lemma [autoref_op_pat_def]: "list_of_eucl ` X ≡ OP op_list_of_eucl_image $ X" by simp

definition [simp]: "op_eucl_of_list_image X = (eucl_of_list ` X::'a::executable_euclidean_space set)"
lemma [autoref_op_pat_def]: "eucl_of_list ` X ≡ OP op_eucl_of_list_image $ X" by simp

definition [simp]: "op_take_image n X = take n ` X"
lemma [autoref_op_pat_def]: "take n ` X ≡ OP op_take_image $ n $ X" by simp

definition [simp]: "op_drop_image n X = drop n ` X"
lemma [autoref_op_pat_def]: "drop n ` X ≡ OP op_drop_image $ n $ X" by simp

definition "approx_slp_appr fas slp X = do {
    cfp ← approx_slp_spec fas DIM('a::executable_euclidean_space) slp X;
    (case cfp of
      Some cfp ⇒ do {
        RETURN ((eucl_of_list ` cfp::'a set):::appr_rel)
      }
      | None ⇒ do {
        SUCCEED
      }
    )
  }"
sublocale autoref_op_pat_def approx_slp_appr .

definition "mig_set D (X::'a::executable_euclidean_space set) = do {
    (i, s) ← op_ivl_rep_of_set (X:::appr_rel);
    let migc = mig_componentwise i s;
    ASSUME (D = DIM('a));
    let norm_fas = ([Norm (map floatarith.Var [0..<D])]:::⟨Id⟩list_rel);
    let env = (list_of_eucl ` ({migc .. migc}:::appr_rel):::appr_rell);
    (n::real set) ← approx_slp_appr  norm_fas (slp_of_fas norm_fas) env;
    (ni::real) ← Inf_spec (n:::appr_rel);
    RETURN (rnv_of_lv ni::real)
  }"
sublocale autoref_op_pat_def mig_set .

lemma appr_rel_br: "appr_rel = br (λxs. eucl_of_list ` (set_of_appr xs)::'a set) (λxs. length xs = DIM('a::executable_euclidean_space))"
  unfolding appr_rel_def lv_rel_def set_rel_br
  unfolding appr_rell_internal br_chain o_def
  using length_set_of_appr
  by (auto dest!: brD intro: brI)

definition print_set::"bool ⇒ 'a set ⇒ unit" where "print_set _ _ = ()"

definition trace_set::"string⇒'a set option⇒unit" where "trace_set _ _ = ()"

definition print_msg::"string ⇒ unit" where "print_msg _ = ()"

abbreviation "CHECKs ≡ λs. CHECK (λ_. print_msg s)"

definition "ncc (X::'a::executable_euclidean_space set) ⟷ X ≠ {} ∧ compact X ∧ convex X"

definition "ncc_precond TYPE('a::executable_euclidean_space) ⟷ (∀(Xi, X::'a set) ∈ appr_rel. ncc X)"

end

lemma prod_relI': "⟦(a,fst ab')∈R1; (b,snd ab')∈R2⟧ ⟹ ((a,b),ab')∈⟨R1,R2⟩prod_rel"
  by  (auto simp: prod_rel_def)

lemma lv_relivl_relI:
  "((xs', ys'), {eucl_of_list xs..eucl_of_list ys::'a::executable_euclidean_space}) ∈ ⟨lv_rel⟩ivl_rel"
  if [simp]: "xs' = xs" "ys' = ys" "DIM('a) = length xs" "length ys = length xs"
  by (force simp: ivl_rel_def set_of_ivl_def
      intro!:  brI lv_relI prod_relI[of _ "eucl_of_list xs" _ _ "eucl_of_list ys"])


end