Theory Refine_Reachability_Analysis

theory Refine_Reachability_Analysis
  imports
    Abstract_Reachability_Analysis
    Refine_Rigorous_Numerics
begin

lemma isDERIV_simps[simp]:
  "isDERIV i 1 xs" "isDERIV i 0 xs"
  "isDERIV i (a + b) xs   isDERIV i a xs  isDERIV i b xs"
  "isDERIV i (a - b) xs  isDERIV i a xs  isDERIV i b xs"
  "isDERIV i (a * b) xs  isDERIV i a xs  isDERIV i b xs"
  "isDERIV i (a / b) xs  isDERIV i a xs  isDERIV i b xs  interpret_floatarith b xs  0"
  "isDERIV i (-a) xs = isDERIV i a xs"
  by (auto simp: one_floatarith_def zero_floatarith_def plus_floatarith_def minus_floatarith_def
      times_floatarith_def divide_floatarith_def uminus_floatarith_def)

lemma list_of_eucl_of_env:
  assumes [simp]: "length xs = DIM('a)"
  shows "(list_of_eucl (eucl_of_env xs vs::'a)) = (map (λi. vs ! (xs ! i)) [0..<DIM('a::executable_euclidean_space)])"
  by (auto intro!: nth_equalityI simp: eucl_of_env_def eucl_of_list_inner)

context approximate_sets_ode
begin

lemma interpret_ode_fa[simp]:
  assumes [simp]: "length xs = DIM('a::executable_euclidean_space)" "length vs  DIM('a)" "length ode_e = DIM('a)"
    and mV:  "max_Var_floatariths ode_e  DIM('a)"
  shows "(eucl_of_list (interpret_floatariths (ode_fa xs) vs)::'a) = ode (eucl_of_list (interpret_floatariths xs vs))"
  unfolding ode_fa_def
  apply (auto intro!: euclidean_eqI[where 'a='a] simp: eucl_of_list_inner ode_def)
  apply (subst interpret_floatarith_subst_floatarith[where D="length vs"])
   apply (auto intro!: max_Var_floatarith_le_max_Var_floatariths_nth[le]
      mV[le])
  apply (rule interpret_floatarith_max_Var_cong)
  apply (drule max_Var_floatariths_lessI) apply simp
  apply (drule less_le_trans[OF _ mV])
  apply auto
  apply (subst nth_map)
   apply simp  using assms(2) order.strict_trans2 apply blast
  apply (subst nth_upt) apply simp
   apply (rule less_le_trans, assumption, simp)
  apply auto
  done

lemma length_ode_fa[simp]: "length (ode_fa xs) = length ode_e"
  by (auto simp: ode_fa_def)

lemma max_Var_ode_fa[le]:
  assumes "max_Var_floatariths ode_e  length xs"
  shows "max_Var_floatariths (ode_fa xs)  max_Var_floatariths xs"
  by (auto simp: ode_fa_def intro!: assms max_Var_floatariths_subst_floatarith_le)

lemma max_Var_floatariths_ode_d_expr[le]:
  "max_Var_floatariths ode_e  d  d > 0 
    max_Var_floatariths (ode_d_expr d m)  3 * d"
  by (auto simp: ode_d_expr_def
      intro!: max_Var_floatarith_FDERIV_n_floatariths[le]
        max_Var_floatarith_FDERIV_floatariths[le])

lemma interpret_ode_d_fa:
  assumes FDERIV: "(eucl_of_list (interpret_floatariths xs vs)::'a::executable_euclidean_space)  Csafe"
  assumes [simp]: "length ds = DIM('a)" "length xs = DIM('a)"
  notes [simp] = safe_length[OF FDERIV]
  shows "(eucl_of_list (interpret_floatariths (ode_d_fa n xs ds) vs)::'a) =
      ode_d n (eucl_of_list (interpret_floatariths xs vs)) (eucl_of_list (interpret_floatariths ds vs))
        (eucl_of_list (interpret_floatariths ds vs))"
  apply (transfer fixing: xs ds vs n)
  using FDERIV apply (auto simp del: isnFDERIV.simps simp: interpret_floatariths_append)
  apply (auto simp add: list_of_eucl_of_env ode_def
      ode_d_raw_def eucl_of_list_inner
      intro!: euclidean_eqI[where 'a='a])
  apply (auto simp: ode_d_fa_def )
  apply (subst interpret_floatarith_subst_floatarith[OF max_Var_floatarith_le_max_Var_floatariths_nth], simp)
  apply (rule interpret_floatarith_max_Var_cong)
  subgoal premises prems for b i
  proof -
    from prems have i: "i < max_Var_floatariths (ode_d_expr DIM('a) n)"
      by (auto dest!: max_Var_floatariths_lessI)
    also have "  3 * DIM('a)"
      by (auto intro!: max_Var_floatariths_ode_d_expr safe_max_Var[OF FDERIV])
    finally have "i < 3 * DIM('a)" .
    then show ?thesis
      using prems i
      by (auto simp: nth_append)
  qed
  done

lemma safe_at_within: "x  Csafe  at x = at x within Csafe"
  by (subst (2) at_within_open)  (auto simp: open_safe)

lemma ivlflowsD:
  assumes "ivlflows stops stopcont trap rsctn" "ivl  (plane_of ` stops) × UNIV "
  shows "ivl  (snd (stopcont ivl))"
    "(fst (stopcont ivl))  sbelow_halfspace rsctn × UNIV"
    "fst (stopcont ivl)  snd (stopcont ivl)"
    "(snd (stopcont ivl))  sbelow_halfspace rsctn × UNIV"
    "flowsto (ivl) {0..} ((snd (stopcont ivl))) ((fst (stopcont ivl))  trap)"
  using assms(1)[unfolded ivlflows_def, rule_format, OF assms(2)]
  by auto

lemma ivlflows_flowsto:
  assumes "ivlflows stops stopcont trap rsctn" "ivl  (plane_of ` stops) × UNIV"
  assumes "stopcont ivl = (x, y)"
  shows "flowsto (ivl) {0..} y (x  trap)"
  using ivlflowsD[OF assms(1,2)] assms(3)
  by auto

lemma ivlflows_emptyI: "ivlflows {} (λx. (x, x)) J K"
  by (auto simp: ivlflows_def set_of_ivl_def)

lemma plane_of_neg[simp]: "plane_of (Sctn (- normal sctn) (- pstn sctn)) = plane_of sctn"
  by (auto simp: plane_of_def)

lemmas safe_max_Var_le[intro] = safe_max_Var[le]

lemmas [simp] = safe_length

lemma continuous_on_ode_d2: "continuous_on (Csafe) ode_d2"
proof -
  have isn: "isnFDERIV DIM('a) ode_e [0..<DIM('a)] [DIM('a)..<2 * DIM('a)] (list_of_eucl x @ list_of_eucl i)
                 (Suc (Suc 0))"
    if "x  Csafe"
    for x i::'a
    by (rule safe_isnFDERIV) fact
  have "continuous_on (Csafe::'a set) (λx. ode_d_raw (Suc 0) x i j)"
     if "i  Basis" "j  Basis" for i j
    apply (rule has_derivative_continuous_on)
    apply (auto simp: ode_d_raw_def at_within_open[OF _ open_safe])
     apply (rule interpret_floatarith_FDERIV_floatariths_append)
    apply (auto simp: ode_d_expr_def 
        intro!: isDERIV_FDERIV_floatariths safe_isFDERIV_append that isFDERIV_map_Var
        safe_max_Var_le
        max_Var_floatarith_FDERIV_floatariths[le])
     apply assumption
    apply arith
    done
  then show ?thesis
    apply (auto intro!: continuous_on_blinfun_componentwise)
    subgoal for i j
    apply (rule continuous_on_eq[where f="(λx. ode_d_raw (Suc 0) x i j)"])
       apply force
      apply (subst ode_d2.rep_eq)
      apply simp
      apply (subst ode_d.rep_eq)
      apply (split if_splits)
      apply (rule conjI) apply simp
      using isn apply force
      done
    done
qed

lemmas continuous_on_ode_d2_comp[continuous_intros] = continuous_on_compose2[OF continuous_on_ode_d2]


lemma map_ode_fa_nth[simp]:
  "d  length ode_e  map (ode_fa_nth CX) [0..<d] = map ((!) (ode_fa CX)) [0..<d]"
  by (auto simp: ode_fa_nth cong: map_cong)

lemma map_ode_d_fa_nth[simp]:
  "d  length ode_e  map (ode_d_fa_nth i CX X) [0..<d] = map ((!) (ode_d_fa i CX X)) [0..<d]"
  by (auto simp: ode_d_fa_nth cong: map_cong)

lemma einterpret_euler_incr_fas:
  assumes "length ode_e = DIM('a)" "length X0 = DIM('a)" "length CX = DIM('a)"
    "DIM('a)  length vs" "max_Var_floatariths ode_e  DIM('a)"
  shows "(einterpret (euler_incr_fas X0 h CX) vs::'a::executable_euclidean_space) =
    einterpret X0 vs + (interpret_floatarith h vs) *R ode (einterpret CX vs)"
  by (simp add: euler_incr_fas_def euler_incr_fas_nth_def assms ode_fa_nth cong: map_cong)

lemma einterpret_euler_err_fas:
  assumes safe: "(einterpret CX vs::'a)  Csafe"
  assumes [simp]: "length X0 = DIM('a)" "length CX = DIM('a)" "DIM('a)  length vs"
  shows "(einterpret (euler_err_fas X0 h CX) vs::'a::executable_euclidean_space) =
    (((interpret_floatarith h vs))2/2) *R ode_d 0 (einterpret CX vs) (ode (einterpret CX vs)) (ode (einterpret CX vs))"
  using safe_length[OF safe] safe_max_Var[OF safe]
  apply (simp add: euler_err_fas_def euler_err_fas_nth_def[abs_def] euler_incr_fas_def)
  apply (subst interpret_ode_d_fa)
  by (auto simp: safe)

lemma einterpret_euler_fas1:
  assumes safe[simp]: "(einterpret CX vs::'a)  Csafe"
  assumes [simp]: "length X0 = DIM('a)" "length CX = DIM('a)" "DIM('a)  length vs"
  shows "(einterpret (take DIM('a) (euler_fas X0 h CX)) vs::'a::executable_euclidean_space) =
    einterpret X0 vs + (interpret_floatarith h vs) *R ode (einterpret X0 vs) +
    (((interpret_floatarith h vs))2/2) *R ode_d 0 (einterpret CX vs) (ode (einterpret CX vs)) (ode (einterpret CX vs))"
  using safe_length[OF safe] safe_max_Var[OF safe]
  by (simp add: euler_fas_def euler_incr_fas_def euler_incr_fas_nth_def[abs_def]
      einterpret_euler_err_fas euler_err_fas_nth_def[abs_def] interpret_ode_d_fa)

lemma einterpret_euler_fas2:
  assumes [simp]: "(einterpret CX vs::'a)  Csafe"
  assumes [simp]: "length X0 = DIM('a)" "length CX = DIM('a)" "DIM('a)  length vs"
  shows "(einterpret (drop DIM('a) (euler_fas  X0 h CX)) vs::'a::executable_euclidean_space) =
    (((interpret_floatarith h vs))2/2) *R ode_d 0 (einterpret CX vs) (ode (einterpret CX vs)) (ode (einterpret CX vs))"
  by (simp add: euler_fas_def euler_incr_fas_def einterpret_euler_err_fas)

lemma ode_d_Suc_0_eq_ode_d2: "x  Csafe  ode_d (Suc 0) x = ode_d2 x"
  unfolding ode_d2.rep_eq by auto

lemma rk2_increment_rk2_fas_err:
  fixes h s1 s2 rkp x0 cx vs
  defines "h'  interpret_floatarith h vs"
  defines "s2'  interpret_floatarith s2 vs"
  defines "rkp'  interpret_floatarith rkp vs"
  defines "x0'  einterpret x0 vs"
  defines "cx'  einterpret cx vs"
  assumes cx_flow: "cx' = flow0 x0' (h' * s1')"
  assumes [simp]: "length x0 = DIM('a)" "length cx = DIM('a)" "DIM('a)  length vs"
  assumes safes: "x0'  Csafe" "cx'  Csafe" "(x0' + (s2' * h' * rkp') *R ode x0') Csafe"
  shows "(einterpret (rk2_fas_err rkp x0 h cx s2) vs::'a::executable_euclidean_space) =
    heun_remainder1 (flow0 x0') ode_na ode_d_na ode_d2_na 0 h' s1' -
    heun_remainder2 rkp' (flow0 x0') ode_na     ode_d2_na 0 h' s2'"
  using safes
  using safe_length[OF safes(1)] safe_max_Var[OF safes(1)]
  apply (auto simp: heun_remainder1_def heun_remainder2_def discrete_evolution_def
      ode_na_def ode_d_na_def ode_d2_na_def rk2_increment x0'_def rkp'_def h'_def s2'_def
      cx'_def euler_incr_fas_def rk2_fas_err_def rk2_fas_err_nth_def[abs_def]
        euler_incr_fas_nth_def[abs_def]
        interpret_ode_d_fa)
  apply (simp add: ode_d1_eq[symmetric] ode_d_Suc_0_eq_ode_d2 inverse_eq_divide)
  apply (simp add: algebra_simps field_simps divide_simps)
  unfolding cx'_def[symmetric] cx_flow x0'_def h'_def
  apply (simp add: algebra_simps)
  done

lemma map_rk2_fas_err_nth[simp]:
  "d = length ode_e  length b = length ode_e  map (rk2_fas_err_nth a b c e f) [0..<d] = map ((!) (rk2_fas_err a b c e f)) [0..<d]"
  unfolding rk2_fas_err_nth_def rk2_fas_err_def
  by (rule map_cong) auto

lemma rk2_increment_rk2_fas1:
  fixes h s1 s2 rkp x0 cx vs
  defines "h'  interpret_floatarith h vs"
  defines "s2'  interpret_floatarith s2 vs"
  defines "rkp'  interpret_floatarith rkp vs"
  defines "x0'  einterpret x0 vs"
  defines "cx'  einterpret cx vs"
  assumes cx_flow: "cx' = flow0 x0' (h' * s1')"
  assumes [simp]: "length x0 = DIM('a)" "length cx = DIM('a)" "DIM('a)  length vs"
  assumes safes: "(x0'::'a) Csafe" "(cx'::'a) Csafe" "(x0' + (s2' * h' * rkp') *R ode x0'::'a) Csafe"
  shows "(einterpret (take DIM('a) (rk2_fas rkp x0 h cx s2)) vs::'a::executable_euclidean_space) =
    discrete_evolution (rk2_increment rkp' (λ_. ode)) h' 0 x0' + (heun_remainder1 (flow0 x0') ode_na ode_d_na ode_d2_na 0 h' s1' -
    heun_remainder2 rkp' (flow0 x0') ode_na ode_d2_na 0 h' s2')"
  using safes  using safe_length[OF safes(1)] safe_max_Var[OF safes(1)]
  apply (auto simp: discrete_evolution_def rk2_fas_def)
  apply (subst rk2_increment_rk2_fas_err[OF cx_flow[unfolded cx'_def x0'_def h'_def]])
  subgoal by simp
  subgoal by simp
  subgoal by simp
  subgoal by (simp add: x0'_def)
  subgoal by (simp add: cx'_def)
  subgoal by (simp add: x0'_def s2'_def h'_def rkp'_def)
  subgoal using [[show_consts, show_sorts, show_types]]
    by (auto simp: x0'_def s2'_def h'_def rkp'_def rk2_increment euler_incr_fas_def
        euler_incr_fas_nth_def[abs_def] inverse_eq_divide)
  done

lemma rk2_increment_rk2_fas2:
  fixes h s1 s2 rkp x0 cx vs
  defines "h'  interpret_floatarith h vs"
  defines "s2'  interpret_floatarith s2 vs"
  defines "rkp'  interpret_floatarith rkp vs"
  defines "x0'  einterpret x0 vs"
  defines "cx'  einterpret cx vs"
  assumes cx_flow: "cx' = flow0 x0' (h' * s1')"
  assumes [simp]: "length x0 = DIM('a)" "length cx = DIM('a)" "DIM('a)  length vs"
  assumes safes: "x0' Csafe" "cx' Csafe" "(x0' + (s2' * h' * rkp') *R ode x0') Csafe"
  shows "(einterpret (drop DIM('a) (rk2_fas rkp x0 h cx s2)) vs::'a::executable_euclidean_space) =
    (heun_remainder1 (flow0 x0') ode_na ode_d_na ode_d2_na 0 h' s1' -
    heun_remainder2 rkp' (flow0 x0') ode_na      ode_d2_na 0 h' s2')"
  using safes
  apply (auto simp: discrete_evolution_def rk2_fas_def)
  apply (subst rk2_increment_rk2_fas_err[OF cx_flow[unfolded cx'_def x0'_def h'_def]])
  subgoal by simp
  subgoal by simp
  subgoal by simp
  subgoal by (simp add: x0'_def)
  subgoal by (simp add: cx'_def)
  subgoal by (simp add: x0'_def s2'_def h'_def rkp'_def)
  subgoal by (auto simp: x0'_def s2'_def h'_def rkp'_def rk2_increment euler_incr_fas_def inverse_eq_divide)
  done

subsubsection ‹safe set relation›

lemma mk_safe[le, refine_vcg]: "wd TYPE('a::executable_euclidean_space)
  mk_safe X  SPEC (λR::'a set. R = X  X  Csafe)"
  unfolding mk_safe_def
  by refine_vcg

lemma mk_safe_coll[le, refine_vcg]: "wd TYPE('a::executable_euclidean_space) 
    mk_safe_coll X  SPEC (λR::'a set. R = X  X  Csafe)"
  unfolding mk_safe_coll_def autoref_tag_defs
  by (refine_vcg FORWEAK_invarI[where I="λa b. X = b  a  a  Csafe"]) auto


lemma ode_set_spec[THEN order.trans, refine_vcg]:
  assumes [refine_vcg]: "wd TYPE('a::executable_euclidean_space)"
  shows "ode_set X  SPEC (λr. ode ` X  (r::'a set))"
  using assms wdD[OF assms(1)]
  unfolding ode_set_def
  apply (refine_vcg)
  subgoal by (auto simp: env_len_def ode_slp_def)
  subgoal premises prems
    using prems(1,2,4-)
    by (auto simp: env_len_def eucl_of_list_prod ode_def)
  done

lemmas fderiv[derivative_intros] = ode_has_derivative_safeI

lemma fderiv2[derivative_intros]:
  "x  Csafe  (ode_d1 has_derivative ode_d2 x) (at x)"
  by (frule ode_d1_has_derivative_safeI)
    (simp add: ode_d_Suc_0_eq_ode_d2)

lemma derivative_within_safe[derivative_intros]:
  "(g has_derivative h) (at x)  (g has_derivative h) (at x within Csafe)"
  by (rule has_derivative_at_withinI)

lemma cont_fderiv: "continuous_on (Csafe) ode_d1"
  by (rule has_derivative_continuous_on) (auto intro!: derivative_intros)

lemmas cont_fderiv'[continuous_intros] = continuous_on_compose2[OF cont_fderiv]

lemma continuous_on_ode1:
  "continuous_on (Csafe) (ode)"
  using fderiv
  by (auto intro!: has_derivative_continuous_on derivative_intros)

lemma continuous_on_ode[continuous_intros]:
  "continuous_on s g  (x. x  s  (g x)  Csafe)  continuous_on s (λx. ode (g x))"
  using continuous_on_ode1
  by (rule continuous_on_compose2) auto

lemma fderiv'[derivative_intros]:
  assumes "(g has_derivative g' y) (at y within X)"
  assumes "(g y)  Csafe"
  shows "((λy. ode (g y)) has_derivative
    (blinfun_apply (ode_d1 (g y)) ∘∘ g') y) (at y within X)"
  using diff_chain_within[OF assms(1) has_derivative_subset[OF fderiv]] assms(2)
  by (simp add: o_def)

lemma fderiv2'[derivative_intros]:
  assumes "(g has_derivative g' y) (at y within X)"
  assumes "(g y)  Csafe"
  shows "((λy. ode_d1 (g y)) has_derivative
    (blinfun_apply (ode_d2 (g y)) ∘∘ g') y) (at y within X)"
  using diff_chain_within[OF assms(1) has_derivative_subset[OF fderiv2]] assms(2)
  by (simp add: o_def)


subsubsection ‹step of Picard iteration›

lemma ncc_interval: "ncc {a .. b::'a::executable_euclidean_space}  a  b"
  by (auto simp: ncc_def)
lemma nonempty_interval: "nonempty {a .. b::'a::executable_euclidean_space}  a  b"
  by (auto simp: nonempty_def)

lemmas [refine_vcg] = Picard_step_def[THEN eq_refl, THEN order.trans]

lemma Basis_list_zero_mem_Basis[simp]:
  "Basis_list ! 0  Basis"
  unfolding Basis_list[symmetric]
  apply (rule nth_mem)
  apply (rule length_Basis_list_pos)
  done

lemma cfuncset_empty_iff:
  fixes l u::"'d::ordered_euclidean_space"
  shows "l  u  cfuncset l u X = {}  X = {}"
  unfolding cfuncset_def Pi_def
proof (safe, goal_cases)
  case hyps: (1 x)
  from x  X
  have "(λ_. x)  {f. x. x  {l..u}  f x  X}  Collect (continuous_on {l..u})"
    by (auto intro!: continuous_on_const)
  then show ?case using hyps by auto
qed auto

lemma lv_ivl_sings: "lv_ivl [x] [y] = (λx. [x]) ` {x .. y}"
  apply (auto simp: lv_ivl_def)
  subgoal for x by (cases x) auto
  done

lemma Picard_step_ivl_refine[le, refine_vcg]:
  assumes [refine_vcg]: "wd TYPE('a::executable_euclidean_space)"
  assumes "(X::'a set)  Csafe"
  assumes "0  h"
  shows "Picard_step_ivl X0 t0 h X  Picard_step X0 t0 h X"
proof -
  have "h'  {t0..t0 + h}  t0  h'" for h' by simp
  then show ?thesis
    unfolding Picard_step_ivl_def Picard_step_def
    apply (refine_vcg, clarsimp_all simp del: atLeastAtMost_iff)
    subgoal using 0  h by simp
    subgoal by (auto simp: euler_incr_slp_def wdD)
    subgoal by (auto simp: euler_incr_fas'_def)
    subgoal for XS l u
      apply (auto simp: lv_ivl_sings nonempty_interval
          simp del: atLeastAtMost_iff
          intro!: add_integral_ivl_bound)
      subgoal for x0 h' phi x h''
        apply (drule bspec, assumption)
        apply (drule bspec[where x="h'' - t0"], force)
      proof goal_cases
        case (1)
        have *: "map ((!) (list_of_eucl b)) [0..<DIM('a) - Suc 0] @ [b  Basis_list ! (DIM('a) - Suc 0)]
          = list_of_eucl b" for b::'a
          apply (auto intro!: nth_equalityI simp: nth_append not_less)
          using Intersection.le_less_Suc_eq by blast
        have "phi x  X" if "x  {t0 .. h''}" for x
          using 1 that by (auto simp: cfuncset_iff)
        have "x0 + (h'' - t0) *R ode b  {l .. u}" if "b  X" for b
        proof -
          from 1(17)[rule_format, OF that] assms(1)
          have "einterpret (euler_incr_fas' D) (list_of_eucl x0 @ (h'' - t0) # list_of_eucl b)  eucl_of_list ` XS"
            by (auto simp: wdD)
          also have "eucl_of_list ` XS  {l .. u}" by fact
          finally show ?thesis
            by (simp add: euler_incr_fas'_def einterpret_euler_incr_fas map_nth_append1 nth_append wdD[OF wd _] *)
        qed
        then have *: "(h'' - t0) *R ode b  {l - x0..u - x0}" if "b  X" for b using that
          by (auto simp: algebra_simps)
        show ?case
          apply (rule *)
          using 1 by (auto simp: cfuncset_iff)
      qed
      subgoal
        using assms(2)
        by (auto intro!: integrable_continuous_real continuous_intros
            simp: cfuncset_iff)
      done
    done
qed

subsubsection ‹Picard iteration›

lemma inf_le_supI[simp]:
  fixes a b c d::"'d::ordered_euclidean_space"
  shows
    "a  c  inf a b  sup c d"
    "a  d  inf a b  sup c d"
    "b  c  inf a b  sup c d"
    "b  d  inf a b  sup c d"
  by (auto simp: eucl_le[where 'a='d] eucl_inf[where 'a='d] eucl_sup[where 'a='d] inf_real_def sup_real_def
    intro!: sum_mono scaleR_right_mono)

lemmas [refine_vcg_def] = do_widening_spec_def

lemma P_iter_spec[le, refine_vcg]:
  assumes "PHI  Csafe"
  assumes "0  h"
  assumes [refine_vcg]: "wd TYPE('a::executable_euclidean_space)"
  shows "P_iter X0 h i PHI 
    SPEC (λr. case r of
        None  True
      | Some (PHI'::'a set)  nonempty PHI'  compact PHI'  (PHI''  PHI'. RETURN (Some PHI'')  Picard_step X0 0 h PHI'))"
  using assms[unfolded autoref_tag_defs]
proof (induction i arbitrary: PHI)
  case 0 then show ?case
    unfolding P_iter.simps
    by refine_vcg
next
  case (Suc i)
  show ?case
    unfolding P_iter.simps
    apply (refine_vcg Suc)
    subgoal by (auto simp: cfuncset_iff Picard_step_def algebra_simps add_increasing2)
    subgoal for lu l u b CX CX' lu' l' u' b'
      apply (simp add: nonempty_interval Picard_step_def)
      apply (safe intro!: exI[where x="{l .. u}"] compact_interval)
      subgoal by (auto simp: nonempty_interval)
      apply (rule subsetD[of CX' "{l .. u}"])
      subgoal
        apply (rule order_trans, assumption)
        unfolding atLeastatMost_subset_iff
        apply (rule disjI2)
        apply (rule conjI)
        subgoal
          apply (rule order_trans[where y="inf l' l - (if b' then ¦l' - l¦ else 0)"])
          by (simp_all split: if_split_asm add: algebra_simps add_increasing2)
        subgoal
          apply (split if_split_asm)
          apply (rule order_trans[where y="sup u' u + ¦u' - u¦"])
          by (auto split: if_split_asm simp add: algebra_simps add_increasing2)
        done
      subgoal by force
      done
    done
qed


subsubsection ‹iterate step size›
lemma Ball_cfuncset_continuous_on:
  "fcfuncset a b X. continuous_on {a..b} f"
  by (simp add: cfuncset_iff)

lemmas one_step_methodD = one_step_method_def[THEN iffD1, rule_format, le]
lemmas one_step_methodI = one_step_method_def[THEN iffD2, rule_format]

lemma cert_stepsize_lemma:
  assumes prems: " 0 < h"
    "X0  {l..u}"
    "Res  Csafe"
    "Res_ivl  Csafe"
    "{l..u}  Csafe"
    "nonempty PHI'"
    "nonempty Res"
    "x0X0.
       x0  Csafe 
       h  existence_ivl0 x0 
       (h'{0..h}. flow0 x0 h'  PHI') 
       x0 + h *R ode x0  PHI'  flow0 x0 h  Res"
    "nonempty Res_ivl"
    "x0X0.
       x0  Csafe 
       (h{0..h}.
           h  existence_ivl0 x0 
           (h'{0..h}. flow0 x0 h'  PHI') 
           x0 + h *R ode x0  PHI'  flow0 x0 h  Res_ivl)"
    "compact PHI'"
    "PHI''  PHI'"
    "RETURN (Some PHI'')  Picard_step X0 0 h PHI'"
  shows "flowpipe0 X0 h h Res_ivl Res"
proof -
  from prems
  have Ps: "RETURN (Some PHI'')  Picard_step X0 0 h PHI'"
    by simp
  from Ps have PHI':
    "compact PHI''" "PHI''  Csafe"
    "xPHI''. x  Csafe"
    "x0 h'' phi. x0  X0  0  h''  h''  h  phi  cfuncset 0 h'' PHI' 
    x0 + integral {0..h''} (λt. ode (phi t))  PHI''"
    by (auto simp: Picard_step_def)
  then obtain cx where cx: "(λt::real. cx)  cfuncset 0 0 PHI'"
    using PHI''  PHI' nonempty PHI' by (auto simp: cfuncset_def continuous_on_const nonempty_def)
  show "flowpipe0 X0 h h Res_ivl Res"
    unfolding flowpipe0_def atLeastAtMost_singleton
  proof safe
    show "0  h" using 0 < h by simp
    show safe_X0: "x  Csafe" if "x  X0" for x using that {l..u}  Csafe X0  {l..u} by blast
    show "x  Csafe" if "x  Res_ivl" for x using prems that
      by auto
    show "x  Csafe" if "x  Res" for x using prems that by (auto simp:)
    fix x0 assume "x0  X0"
    from PHI'(4)[OF x0  X0 order_refl 0  h cx]
    have "x0  PHI''" by simp
    have *: "h  existence_ivl0 x0" "s  {0 .. h}  flow0 x0 s  PHI''" for s
      using x0  X0 PHI''  PHI' 0  h PHI'(3) x0  PHI''
      by (auto
          simp: cfuncset_def Pi_iff closed_segment_eq_real_ivl ivl_integral_def
          intro!: Picard_iterate_mem_existence_ivlI[OF UNIV_I _ UNIV_I compact PHI''
            x0  PHI'' PHI''  Csafe] PHI'(4)) force+
    show h_ex: "h  existence_ivl0 x0" by fact
    have cf: "(λt::real. x0)  cfuncset 0 h PHI'" for h
      using x0  PHI'' PHI''  PHI'
      by (auto simp: cfuncset_def continuous_intros)
    have mem_PHI': "x0 + h' *R ode x0  PHI'" if "0  h'" "h'  h" for h'
      using that PHI''  PHI' PHI'(4)[OF x0  X0 that cf]
      by auto
    from this prems safe_X0
    show "flow0 x0 h  Res"
      using 0  h h_ex * PHI''  PHI' x0  X0
      by (auto simp: subset_iff dest!: bspec[where x=x0])
    fix h' assume h': "h'  {0..h}"
    then have "h'  existence_ivl0 x0"
      by (meson atLeastAtMost_iff existence_ivl_zero h_ex is_interval_1
          local.is_interval_existence_ivl local.mem_existence_ivl_iv_defined(2))
    from h' this prems safe_X0
    show "flow0 x0 h'  Res_ivl"
      using 0 < h h_ex * PHI''  PHI' x0  X0 mem_PHI' x0  PHI''
      by (auto simp: subset_iff dest!: bspec[where x=x0])
  qed
qed

lemma cert_stepsize_spec[le,refine_vcg]:
  assumes "h > 0"
  assumes "one_step_method m"
  assumes [refine_vcg]: "wd TYPE('a::executable_euclidean_space)"
  shows "cert_stepsize m X0 h i n  SPEC (λ(h', RES::'a set, RES_ivl, _). nonempty RES  nonempty RES_ivl  0 < h'  h'  h  flowpipe0 X0 h' h' RES_ivl RES)"
  using assms(1)[unfolded autoref_tag_defs]
proof (induction n arbitrary: h)
  case 0 then show ?case by simp
next
  case (Suc n)
  note [refine_vcg] = Suc.IH[of "h/2", le]
  show ?case
    unfolding cert_stepsize.simps
    using Suc.prems
    by (refine_vcg Ball_cfuncset_continuous_on one_step_methodD[OF assms(2)])
       (clarsimp_all simp: cert_stepsize_lemma)
qed


subsubsection ‹Euler step›

lemma embed_real_ivl_iff[simp]:
   "(x  {0 .. h *R (One::'a::executable_euclidean_space)}. P (x  hd Basis_list))  (x  {0 .. h}. P x)"
proof (auto simp: eucl_le[where 'a='a], goal_cases)
  case hyps: (1 x)
  have "x = x *R (One::'a)  hd Basis_list"
    by auto
  also have "P "
    apply (rule hyps[rule_format])
    using hyps
    by (auto simp: eucl_le[where 'a='a])
  finally show ?case .
qed

lemma convex_on_segmentI:
  assumes mem_convex: "convex C" "a  C" "a + j *R b  C"
  assumes "0  i" "i  j"
  shows "a + i *R b  C"
proof -
  have "a + i *R b = (1 - i / j) *R a + (i / j) *R (a + j *R b)"
    using assms
    by (auto simp: algebra_simps diff_divide_distrib)
  also have "  C"
    using assms
    by (auto simp: divide_simps intro!: convexD[OF mem_convex])
  finally show ?thesis .
qed

lemma one_step_flowpipe:
  assumes [THEN one_step_methodD, refine_vcg]: "one_step_method m"
  assumes [refine_vcg]: "wd TYPE('a::executable_euclidean_space)"
  shows "one_step X0 h m  SPEC (λ(h', _, RES_ivl, RES::'a set). 0 < h'  h'  h  flowpipe0 X0 h' h' RES_ivl RES)"
  using assms
  unfolding one_step_def
  by refine_vcg

lemma ncc_imageD:
  assumes "ncc ((λx. x ! i) ` env)"
  assumes "nth_image_precond env i"
  shows "compact ((λx. x ! i) ` env::real set)" "closed ((λx. x ! i) ` env)" "bounded ((λx. x ! i) ` env)"
    "((λx. x ! i) ` env)  {}" "convex ((λx. x ! i) ` env)"
  using assms
  by (auto simp: ncc_def nth_image_precond_def compact_eq_bounded_closed)

lemma max_Var_floatariths_ode_d_fa[le]:
  assumes [simp]: "length ode_e > 0" "max_Var_floatariths ode_e  length ode_e"
    "length cxs = length ode_e" "length ys = length ode_e"
  shows "max_Var_floatariths (ode_d_fa i cxs ys)  max (max_Var_floatariths (cxs)) (max_Var_floatariths ys)"
  apply (auto simp: ode_d_fa_def max_Var_floatariths_Max )
  using assms apply auto[1]
  apply (auto intro!: max_Var_floatarith_subst_floatarith_le max_Var_floatariths_ode_d_expr
      max_Var_floatarith_le_max_Var_floatariths_nthI max_Var_ode_fa
      simp: in_set_conv_nth)
   apply (auto simp: max_Var_floatariths_Max in_set_conv_nth)
  done

lemma max_Var_floatariths_euler_err_fas[le]:
  assumes nz: "0 < length ode_e"
    and [simp]: "max_Var_floatariths ode_e  length ode_e"
    "length xs = length ode_e"
    "length cxs = length ode_e"
  shows "max_Var_floatariths (euler_err_fas xs h cxs)
     max (max_Var_floatariths xs) (max (max_Var_floatarith h) (max_Var_floatariths cxs))"
  using nz
  by (auto simp: euler_err_fas_def[abs_def] euler_err_fas_nth_def[abs_def] map_nth_eq_self simp del: length_0_conv
      intro!: max_Var_floatariths_ode_d_fa max_Var_floatariths_map_times
        max_Var_floatariths_map_const max_Var_ode_fa; arith)

lemma max_Var_floatariths_euler_incr_fas[le]:
  assumes [simp]: "max_Var_floatariths ode_e  length ode_e"
    "length xs = length ode_e"
    "length cxs = length ode_e"
  shows "max_Var_floatariths (euler_incr_fas xs h cxs)
     max (max_Var_floatariths xs) (max (max_Var_floatarith h) (max_Var_floatariths cxs))"
  using length_ode_fa
  by (auto simp: euler_incr_fas_def euler_incr_fas_nth_def[abs_def] simp del: length_ode_fa
      intro!: max_Var_floatariths_ode_d_fa max_Var_floatariths_map_plus max_Var_floatariths_map_times
      max_Var_floatariths_map_const max_Var_ode_fa)

lemma map_euler_incr_fas_nth: "length X0 = d  map (euler_incr_fas_nth X0 h CX) [0..<d] = euler_incr_fas X0 h CX"
  by (auto simp: euler_incr_fas_def)

lemma map_euler_err_fas_nth: "length X0 = d  map (euler_err_fas_nth X0 h CX) [0..<d] = euler_err_fas X0 h CX"
  by (auto simp: euler_err_fas_def)

lemma max_Var_floatariths_euler_fas[le]:
  assumes [simp]: "max_Var_floatariths ode_e  length ode_e"
    "length xs = length ode_e"
    "length cxs = length ode_e"
  assumes nz: "0 < length ode_e"
  shows "max_Var_floatariths (euler_fas xs h cxs)  Max {max_Var_floatariths xs, max_Var_floatarith h, max_Var_floatariths cxs}"
  using nz
  by (auto simp: euler_fas_def map_euler_incr_fas_nth map_euler_err_fas_nth
      intro!: max_Var_floatariths_map_plus max_Var_floatariths_euler_incr_fas
        max_Var_floatariths_euler_err_fas)

lemma take_interpret_floatariths:
  "d < length fas  take d (interpret_floatariths fas vs) = interpret_floatariths (take d fas) vs"
  by (auto intro!: nth_equalityI)

lemma length_euler_slp_le: "2 * D  length euler_slp"
  by (auto simp: euler_fas'_def euler_slp_def intro!: order_trans[OF _ length_slp_of_fas_le])

lemma ncc_nonempty[simp]: "ncc x  nonempty x"
  by (simp add: ncc_def nonempty_def)

lemma nccD:
  assumes "ncc X"
  shows "compact X" "closed X" "bounded X" "X  {}" "convex X"
  using assms
  by (auto simp: ncc_def nth_image_precond_def compact_eq_bounded_closed)

lemma D_DIM_wdD[simp]: "wd TYPE('a::executable_euclidean_space)  D = DIM('a)"
  by (auto simp: wdD)

lemma euler_step_flowpipe:
  includes floatarith_syntax
  assumes [refine_vcg]: "wd TYPE('a::executable_euclidean_space)"
  shows "euler_step X0 h  SPEC (λ(h', _, RES_ivl, RES::'a set). 0 < h'  h'  h  flowpipe0 X0 h' h' RES_ivl RES)"
  unfolding euler_step_def THE_NRES_def
  apply (intro SPEC_rule_conjI one_step_flowpipe one_step_methodI)
    apply (refine_vcg, clarsimp_all)
  subgoal using assms by (auto simp: euler_slp_def euler_fas'_def)
  subgoal by (auto simp: euler_slp_def euler_fas'_def)
  subgoal using length_euler_slp_le assms by (auto simp: env_len_def wdD[OF wd _])
  subgoal using length_euler_slp_le assms by (auto simp: env_len_def wdD[OF wd _])
proof (goal_cases)
  case hyps: (1 X0 CX hl hu env res b x0 enve h)
  then interpret derivative_on_prod "{0 .. h}" CX "λ_. ode" "λ(t, x). ode_d1 x oL snd_blinfun"
    by unfold_locales (auto intro!: continuous_intros derivative_eq_intros
        simp: split_beta' subset_iff  wdD[OF wd _])
  from h  existence_ivl0 x0 have s_ex: "s  existence_ivl0 x0" if "0  s" "s  h" for s
    by (metis (no_types, lifting) atLeastAtMost_iff ivl_subset_existence_ivl subset_iff that)
  have "flow0 x0 (h) = flow0 x0 (0 + (h))" by simp
  also have "  eucl_of_list ` take D ` env"
    using hyps
    apply (intro euler_consistent_traj_set[where x="flow0 x0" and u = "h"])
            apply (auto intro!: 0  h flow_has_vector_derivative[THEN has_vector_derivative_at_within]
        simp: nccD discrete_evolution_def euler_increment subset_iff wdD[OF wd _]
          Let_def s_ex min_def max_def lv_ivl_sings)
    subgoal premises prems for s
    proof -
      have "interpret_floatariths (euler_fas' DIM('a)) (list_of_eucl x0 @ list_of_eucl (flow0 x0 s) @ [h])  env"
        using prems
        by (auto intro!: prems(1)[rule_format])
      then have "eucl_of_list (take D (interpret_floatariths (euler_fas' DIM('a)) (list_of_eucl x0 @ list_of_eucl (flow0 x0 s) @ [h])))
         eucl_of_list ` take D ` env"
        (is "eucl_of_list (take _ (interpret_floatariths  _ ?vs))  _")
        by auto
      also
      have "take (2 * D) (interpret_floatariths (euler_fas' DIM('a)) ?vs) =
        interpret_floatariths (map fold_const_fa (euler_fas (map floatarith.Var [0..<D]) (Var (2 * D)) (map floatarith.Var [D..<2 * D]))) ?vs"
        unfolding euler_fas'_def
        by (auto simp: euler_fas_def wdD[OF wd _] simp del: map_map
            intro!: max_Var_floatariths_map_plus max_Var_floatariths_euler_incr_fas
              max_Var_floatariths_euler_err_fas wd _
              max_Var_floatariths_fold_const_fa[le])
      then have "take D (take (2 * D) (interpret_floatariths (euler_fas' DIM('a)) ?vs)) =
        take D (interpret_floatariths (euler_fas  (map floatarith.Var [0..<D]) (Var(2 * D)) (map floatarith.Var [D..<2 * D])) ?vs)"
        by simp
      then have "take D (interpret_floatariths (euler_fas' DIM('a)) ?vs) =
        take DIM('a) (interpret_floatariths (euler_fas  (map floatarith.Var [0..<D]) (Var(2 * D)) (map floatarith.Var [D..<2 * D])) ?vs)"
        by (simp add: wdD[OF wd _])
      also have "eucl_of_list  =
          x0 + h *R ode x0 + (h2 / 2) *R (ode_d 0 (flow0 x0 s) (ode (flow0 x0 s))) (ode (flow0 x0 s))"
        by (auto simp: take_interpret_floatariths einterpret_euler_fas1 map_nth_append1 prems nth_append
            wdD[OF wd _])
      finally show ?thesis
        by (simp add: prems(10) prems(13) prems(14) prems(5) ode_d1_eq[symmetric] wdD[OF wd _])
    qed
    done
  also have "  res" using assms hyps by auto
  finally show ?case by simp
qed (auto simp: assms)

lemma length_rk2_slp_le: "2 * D  length rk2_slp"
  by (auto simp: rk2_slp_def rk2_fas'_def intro!: order_trans[OF _ length_slp_of_fas_le])

lemma max_Var_floatarith_Re[simp]: "max_Var_floatarith (Re x) = 0"
  by (auto simp: Re_def split: prod.splits)

lemma max_Var_floatariths_rk2_fas_err[le]:
  assumes nz: "0 < length ode_e"
    and [simp]: "max_Var_floatariths ode_e  length ode_e" "length x0 = length ode_e" "length cx = length ode_e"
  shows "max_Var_floatariths (rk2_fas_err rkp x0 h cx s2) 
    Max {max_Var_floatarith rkp, max_Var_floatariths x0, max_Var_floatarith h, max_Var_floatariths cx,
      max_Var_floatarith s2}"
  using nz
  unfolding rk2_fas_err_def rk2_fas_err_nth_def
  by (auto simp: rk2_fas_err_def
      intro!: max_Var_floatariths_append max_Var_floatariths_map_plus max_Var_floatariths_map_times
        max_Var_floatariths_map_const max_Var_ode_fa max_Var_floatariths_euler_incr_fas
        max_Var_floatariths_ode_d_fa; arith)

lemma max_Var_floatarith_one[simp]: "max_Var_floatarith 1 = 0"
  and max_Var_floatarith_zero[simp]: "max_Var_floatarith 0 = 0"
  by (auto simp: one_floatarith_def zero_floatarith_def)

lemma max_Var_floatariths_rk2_fas[le]:
  assumes nz: "0 < length ode_e"
    and [simp]: "max_Var_floatariths ode_e  length ode_e" "length x0 = length ode_e" "length cx = length ode_e"
  shows "max_Var_floatariths (rk2_fas rkp x0 h cx s2) 
    Max {max_Var_floatarith rkp, max_Var_floatariths x0, max_Var_floatarith h, max_Var_floatariths cx,
      max_Var_floatarith s2}"
  using nz
  by (auto simp: rk2_fas_def
      intro!: max_Var_floatariths_append max_Var_floatariths_map_plus max_Var_floatariths_map_times
        max_Var_floatariths_map_const max_Var_ode_fa max_Var_floatariths_euler_incr_fas
        max_Var_floatariths_rk2_fas_err)

lemma rk2_step_flowpipe:
  includes floatarith_syntax
  assumes [refine_vcg]: "wd TYPE('a::executable_euclidean_space)"
  shows "rk2_step X0 h  SPEC (λ(h', _, RES_ivl, RES::'a set).
    0 < h'  h'  h  flowpipe0 X0 h' h' RES_ivl RES)"
  unfolding rk2_step_def THE_NRES_def
  apply (intro one_step_flowpipe assms one_step_methodI)
  apply (refine_vcg, clarsimp_all)
  subgoal using assms by (auto simp: rk2_slp_def rk2_fas'_def)
  subgoal by (auto simp: rk2_slp_def rk2_fas'_def)
  subgoal using length_rk2_slp_le by (auto simp: env_len_def wdD[OF wd _])
  subgoal using length_rk2_slp_le by (auto simp: env_len_def wdD[OF wd _])
proof (goal_cases)
  case hyps: (1 X0 CX hl hu rk2_param env res b x0 el h)
  from assms have "D = DIM('a)" by simp
  have aux: "ode (flow0 x0 s) = ode (snd (s, flow0 x0 s))" for s
    by simp
  from hyps interpret derivative_on_prod "{0 .. h}" CX "λ_ x. ode x" "λ(t, x). ode_d1 x oL snd_blinfun"
    by unfold_locales
      (auto intro!: continuous_intros derivative_eq_intros simp: split_beta' subset_iff)

  have aux2: "blinfun_apply (ode_d1 (snd tx))  snd = blinfun_apply (ode_d1 (snd tx) oL snd_blinfun)"
    for tx::"real×'a"
    by (auto intro!: blinfun_eqI)

  have aux3: "blinfun_apply (ode_d2 (snd tx)) (snd h) oL snd_blinfun =
    (flip_blinfun (flip_blinfun (ode_d2 (snd tx) oL snd_blinfun) oL snd_blinfun)) h"
    for tx h::"real×'a"
    by (auto intro!: blinfun_eqI)


  have "flow0 x0 (h) = flow0 x0 (0 + (h))" by simp
  also have "  eucl_of_list ` take D ` env"
    using hyps assms
    apply (intro rk2_consistent_traj_set[where
      x="flow0 x0" and u = "h" and T="{0..h}" and X="CX" and p="rk2_param"
      and f = "ode_na" and f' = ode_d_na and g' = ode_d_na and f'' = ode_d2_na and g'' = ode_d2_na])
    subgoal by (simp add: 0  h)
    subgoal by simp
    subgoal by simp
    subgoal by auto
    subgoal by (auto simp add: ncc_def nonempty_def)
    subgoal
      apply (rule flow_has_vector_derivative[THEN has_vector_derivative_at_within, THEN has_vector_derivative_eq_rhs])
      subgoal by (metis (no_types, lifting) ivl_subset_existence_ivl subset_iff)
      subgoal by (force simp: ode_na_def[abs_def] ode_d_na_def[abs_def] ode_d2_na_def[abs_def])
      done
    subgoal
      unfolding ode_na_def ode_d_na_def ode_d2_na_def
      apply (rule derivative_eq_intros)
        apply (rule derivative_intros)
        apply (rule derivative_intros)
      subgoal by (auto simp: ncc_def nonempty_def)
      subgoal by force
      done
    subgoal
      unfolding ode_na_def ode_d_na_def ode_d2_na_def
      apply (rule derivative_eq_intros)
        apply (rule derivative_intros)
         apply (rule derivative_intros)
         apply (rule derivative_intros)
      subgoal by (force simp: nonempty_def)
       apply (rule derivative_intros)
      subgoal by (auto intro!: aux3)
      done
    subgoal by (rule refl)
    subgoal by (rule refl)
    subgoal
      apply (rule compact_imp_bounded)
      apply (rule compact_continuous_image)
      subgoal
        by (auto intro!: continuous_intros simp: ode_na_def ode_d_na_def ode_d2_na_def)
      subgoal by (auto simp: ncc_def intro!: compact_Times)
      done
    subgoal by auto
    subgoal by simp
    subgoal by simp
    subgoal
      apply (rule convex_on_segmentI[where j=h])
      using mult_left_le_one_le[of h "rk2_param"]
      by (auto simp: ncc_def mult_left_le_one_le mult_le_one ac_simps ode_na_def
          ode_d_na_def ode_d2_na_def dest: bspec[where x=0])
    subgoal by (simp add: ncc_def)
    subgoal by (simp add: ncc_def compact_imp_closed)
    subgoal for s1 s2
      apply (clarsimp simp add: lv_ivl_sings)
      subgoal premises prems
      proof -
        have "s2 * rk2_param * h  h"
          apply (rule mult_left_le_one_le)
          using assms prems
          by (auto intro!: mult_le_one)
        then have s2: "(s2 * h * rk2_param)  {0 .. h}"
          using prems assms by (auto simp: ac_simps)
        have s1: "h * s1  {0 .. h}" using prems
          by (auto intro!: mult_right_le_one_le)
        then have
          "interpret_floatariths (rk2_fas' D)
              (list_of_eucl x0 @ list_of_eucl (flow0 x0 (h * s1)) @ [rk2_param, h, s2])  env"
          unfolding D = _ using prems
          by (intro prems(17)[rule_format]) auto
        then have "take (2 * D) (interpret_floatariths (rk2_fas' D)
              (list_of_eucl x0 @ list_of_eucl (flow0 x0 (h * s1)) @ [rk2_param, h, s2]))  take (2 * D) ` env"
          (is "?l  _")
          by auto
        also have "?l = interpret_floatariths
         (map fold_const_fa (rk2_fas (Var (2 * D)) (map floatarith.Var [0..<D]) (Var (2 * D + 1))
          (map floatarith.Var [D..<2 * D])
           (Var (2 * D + 2))))
         (list_of_eucl x0 @ list_of_eucl (flow0 x0 (h * s1)) @ [rk2_param, h, s2])"
          (is "_ = interpret_floatariths (map fold_const_fa ?fas) ?xs")
          unfolding rk2_fas'_def
          by (auto intro!: max_Var_floatariths_rk2_fas max_Var_floatariths_fold_const_fa[le] simp: wdD[OF wd _])
        finally have "take D (interpret_floatariths ?fas ?xs)  take D ` take (2 * D) ` env"
          by auto
        also have " = take D ` env" by (auto simp: image_image wdD[OF wd _])
        finally have "eucl_of_list (take D (interpret_floatariths ?fas ?xs))  eucl_of_list ` take D ` env"
          by simp
        then have "einterpret (take D ?fas) ?xs  eucl_of_list ` take D ` env"
          by (simp add: take_interpret_floatariths wdD[OF wd _])
        also have "einterpret (take D ?fas) ?xs =
          discrete_evolution (rk2_increment (rk2_param) (λt x. ode_na (t, x))) h 0 x0 +
          heun_remainder1 (flow0 x0) ode_na ode_d_na ode_d2_na 0 h s1 -
          heun_remainder2 (rk2_param) (flow0 x0) ode_na ode_d2_na 0 h s2"
          apply (simp add: wdD[OF wd _])
          apply (subst rk2_increment_rk2_fas1[where ?s1'.0 = s1])
          subgoal by (auto simp: nth_append map_nth_append1)
          subgoal by auto
          subgoal by auto
          subgoal by auto
          subgoal by (auto simp: nth_append map_nth_append1 x0  Csafe)
          subgoal
            apply (auto simp: nth_append map_nth_append1 x0  Csafe)
            by (meson connectedD_interval existence_ivl_zero flow0_defined hyps
                mult_right_le_one_le mult_sign_intros(1) mvar.connected prems)
          subgoal
          proof -
            have "x0 + ((rk2_param * s2) * h) *R ode x0  CX"
              by (rule convex_on_segmentI[where j=h])
                 (use prems in auto simp: ncc_def mult_left_le_one_le mult_le_one
                  dest: bspec[where x=0])
            also note   Csafe
            finally show ?thesis
              by (auto simp: nth_append map_nth_append1 x0  Csafe ac_simps)
          qed
          subgoal by (auto simp: nth_append map_nth_append1 ode_na_def)
          done
        finally show ?thesis by (simp add: D = _)
      qed
      done
    done
  also have "  res" using hyps(6) by (simp add: D = _)
  finally show ?case by (simp add: D = _)
qed

lemma interpret_adapt_stepsize_fa:
  "interpret_floatarith (adapt_stepsize_fa rtol m_id e h') []
    = float_of h' * (float_of(rtol) / float_of e) powr (1 / (float_of (m_id) + 1))"
  by (auto simp: inverse_eq_divide adapt_stepsize_fa_def)

lemma choose_step_flowpipe[le, refine_vcg]:
  assumes "wd TYPE('a::executable_euclidean_space)"
  shows "choose_step X0 h  SPEC (λ(h', _, RES_ivl, (RES::'a set)). 0 < h'  h'  h  flowpipe0 X0 h' h' RES_ivl RES)"
  using assms
  unfolding choose_step_def
  by (refine_vcg rk2_step_flowpipe euler_step_flowpipe)

lemma CsafeI: "t  existence_ivl0 x  x  Csafe"
  using local.mem_existence_ivl_iv_defined(2) by blast

lemma apply_vareq: "blinfun_apply (vareq x t) = ode_d1 (flow0 x t)"
  by (auto simp: vareq_def)

lemma Dflow_has_derivative:
  "t  existence_ivl0 x  (Dflow x has_derivative blinfun_scaleR_left (ode_d1 (flow0 x t) oL Dflow x t)) (at t)"
  by (auto simp: Dflow_def blinfun.bilinear_simps scaleR_blinfun_compose_left apply_vareq CsafeI
      intro!: derivative_eq_intros mvar.flow_has_derivative[THEN has_derivative_eq_rhs] ext
        blinfun_eqI)

lemma matrix_scaleR: "matrix (blinfun_apply (h *R X)) = h *R matrix X"
  by (vector matrix_def blinfun.bilinear_simps)

lemma blinfun_of_vmatrix_matrix_matrix_mult[simp]:
  "blinfun_of_vmatrix (A ** B) = blinfun_of_vmatrix A oL blinfun_of_vmatrix B"
  including blinfun.lifting
  by transfer (auto simp: o_def matrix_vector_mul_assoc)

lemma blinfun_of_vmatrix_mat_1[simp]: "blinfun_of_vmatrix (mat 1) = 1L"
  including blinfun.lifting
  by transfer (auto simp: matrix_vector_mul_lid)

lemma blinfun_of_vmatrix_matrix[simp]:
  "blinfun_of_vmatrix (matrix (blinfun_apply A)) = A"
  including blinfun.lifting
  by transfer (auto simp: bounded_linear.linear matrix_works)

lemma inner_Basis_eq_vec_nth: "b  Basis  v  b = vec_nth v (enum_class.enum ! index Basis_list b)"
  by (auto simp: inner_vec_def vec_nth_Basis if_distrib Basis_vec_def axis_eq_axis
        index_Basis_list_axis1
      cong: if_cong)

lemma intersects_sctns_spec_nres[le, refine_vcg]:
  "intersects_sctns X' sctns  intersects_sctns_spec X' sctns"
  unfolding intersects_sctns_spec_def intersects_sctns_def
  by refine_vcg auto

lemma intersects_sections_spec_clw_ref[le, refine_vcg]:
  "intersects_sctns_spec_clw R sctns  intersects_sctns_spec R sctns"
  unfolding intersects_sctns_spec_def intersects_sctns_spec_clw_def
  by (refine_vcg FORWEAK_mono_rule[where I="λS b. ¬b  S  (plane_of ` sctns) = {}"]) auto

lemma eq_nth_iff_index:
  "distinct xs  n < length xs  i = xs ! n   index xs i = n"
  using index_nth_id by fastforce

lemma
  max_Var_floatariths_ode_e_wd:
  assumes wd: "wd (TYPE('n::enum rvec))"
  assumes "CARD('n)  K"
  shows "max_Var_floatariths ode_e  K"
  using wdD[OF wd] assms by auto

lemma nonzero_component[le, refine_vcg]: "nonzero_component s X n  SPEC (λ_. bX. b  n  0)"
  unfolding nonzero_component_def
  by refine_vcg auto

lemma
  interpret_slp_env_lenD:
  assumes "cxCX. interpret_slp (slp_of_fas (fas)) (env cx)  R"
  assumes "cx  CX"
  shows "interpret_floatariths fas (env cx)  take (length fas) ` R"
proof -
  from slp_of_fas
  have "take (length fas) (interpret_slp (slp_of_fas fas) (env cx)) = interpret_floatariths fas (env cx)"
    by auto
  moreover
  from assms(1)[rule_format, OF cx  CX]
  have "interpret_slp (slp_of_fas fas) (env cx)  R" by auto
  ultimately show ?thesis by force
qed

lemma flowpipe0_imp_flowpipe:
  assumes "flowpipe0 (fst ` X0) x1 x1 aba bba"
  shows "flowpipe X0 x1 x1 (aba × UNIV) (bba × UNIV)"
  using assms
  by (auto simp: flowpipe0_def flowpipe_def)

lemma disjoints_spec[le, refine_vcg]:
  "disjoints_spec X Y  SPEC (λb. b  (X  Y = {}))"
  unfolding disjoints_spec_def autoref_tag_defs
  by refine_vcg auto

lemma inner_eq_zero_abs_BasisI:
  "¦y¦  Basis  b  Basis  ¦y¦  b  y  b = 0"
  for b y::"'a::executable_euclidean_space"
  by (metis abs_inner inner_Basis linorder_not_le order_refl zero_less_abs_iff)

lemma abs_in_Basis_absE:
  fixes x y::"'a::executable_euclidean_space"
  assumes "abs y  Basis"
  obtains "abs y = y" | "abs y = -y"
proof -
  have "abs y = (iBasis. (abs (y  i)) *R i)"
    by (simp add: euclidean_representation abs_inner[symmetric] assms)
  also have "Basis = insert (abs y) (Basis - {abs y})" using assms by auto
  also have "(iinsert ¦y¦ (Basis - {¦y¦}). ¦y  i¦ *R i) = ¦y  ¦y¦¦ *R ¦y¦"
    apply (subst sum.insert)
    using assms
    by (auto simp: abs_inner[symmetric] inner_Basis if_distribR if_distrib
        cong: if_cong)
  finally have "¦y¦ = ¦y  ¦y¦¦ *R ¦y¦" by simp
  moreover have " = y   = - y"
    using assms
    by (auto simp: abs_real_def algebra_simps  intro!: euclidean_eqI[where 'a='a]
        simp: inner_Basis inner_eq_zero_abs_BasisI split: if_splits)
  ultimately consider "¦y¦ = y" | "¦y¦ = - y" by auto
  then show ?thesis
    by (cases; rule that)
qed

lemma abs_in_BasisE:
  fixes x y::"'a::executable_euclidean_space"
  assumes "abs y  Basis"
  obtains i where "i  Basis" "y = i" | i where "i  Basis" "y = -i"
proof -
  from abs_in_Basis_absE[OF assms]
  consider "¦y¦ = y" | "¦y¦ = - y"
    by auto
  then show ?thesis
  proof cases
    case 1 with assms have "abs y  Basis" "y = abs y" by auto
    then show ?thesis ..
  next
    case 2
    with assms have "abs y  Basis" "y = - abs y" by auto
    then show ?thesis ..
  qed
qed

lemma subset_spec_plane[le, refine_vcg]:
  "subset_spec_plane X sctn  SPEC (λb. b  X  plane_of sctn)"
  unfolding subset_spec_plane_def
  by (refine_vcg) (auto simp: plane_of_def eucl_le[where 'a='a] dest!: bspec elim!: abs_in_BasisE)

lemma subset_spec_coll_refine[le, refine_vcg]: "subset_spec_coll X Y  subset_spec X Y"
  unfolding subset_spec_coll_def autoref_tag_defs subset_spec_def
  by (refine_vcg FORWEAK_mono_rule[where I="λX b. b  X  Y"]) auto

lemma
  eventually_in_planerectI:
  fixes n::"'a::executable_euclidean_space"
  assumes "abs n  Basis"
  assumes "{l .. u}  plane n c" "l  u"
  assumes "i. i  Basis  i  abs n  l  i < x  i"
  assumes "i. i  Basis  i  abs n  x  i < u  i"
  shows "F x in at x within plane n c. x  {l .. u}"
proof -
  have "F x in at x within plane n c. x  plane n c"
    unfolding eventually_at_filter
    by simp
  then have "F x in at x within plane n c. l  abs n  x  abs n  x  abs n  u  abs n"
    apply eventually_elim
    using assms(1,2,3)
    by (auto simp: elim!: abs_in_BasisE)
  moreover
  {
    fix i assume that: "i  Basis" "i  abs n"
    have "F x in at x within plane n c. l  i < x  i" "F x in at x within plane n c. x  i < u  i"
      by (auto intro!: order_tendstoD assms tendsto_eq_intros that)
    then have "F x in at x within plane n c. l  i < x  i  x  i < u  i"
      by eventually_elim auto
  } then have "F x in at x within plane n c. i  Basis - {abs n}. l  i < x  i  x  i < u  i"
    by (auto intro!: eventually_ball_finite)
  then have "F x in at x within plane n c. i  Basis - {abs n}. l  i  x  i  x  i  u  i"
    by eventually_elim (auto intro!: less_imp_le)
  ultimately
  have "F x in at x within plane n c. iBasis. l  i  x  i  x  i  u  i"
    by eventually_elim auto
  then show ?thesis
    by eventually_elim (auto simp: eucl_le[where 'a='a])
qed

lemma mem_ivl_euclI: "k  {c..d::'x::ordered_euclidean_space}" if "i. i  Basis  c  i  k  i" "i. i  Basis  k  i  d  i"
  using that
  by (auto simp: eucl_le[where 'a='x])

lemma op_eventually_within_sctn[le, refine_vcg]:
  "op_eventually_within_sctn X sctn S 
    SPEC (λb. b  (x  X. x  S  (F x in at x within plane_of sctn. x  S)))"
  unfolding op_eventually_within_sctn_def
  apply refine_vcg
  unfolding plane_of_def autoref_tag_defs
  apply (safe intro!: eventually_in_planerectI mem_ivl_euclI)
  subgoal premises prems for a b c d e f g h i j k B
  proof cases
    assume "B = ¦normal sctn¦"
    moreover
    have "c  plane (normal sctn) (pstn sctn)" "k  plane (normal sctn) (pstn sctn)"
      using prems by auto
    ultimately show "c  B  k  B" using ¦normal sctn¦  set Basis_list
      by (auto simp: elim!: abs_in_Basis_absE)
  next
    assume B: "B  ¦normal sctn¦"
    have "k  B  {g  B .. h   B}"
      using k  X X  {g..h} B  Basis by (auto simp: eucl_le[where 'a='a])
    with B prems show ?thesis by (auto dest!: bspec elim!: abs_in_Basis_absE)
  qed
  subgoal premises prems for a b c d e f g h i j k B
  proof cases
    assume "B = ¦normal sctn¦"
    moreover
    have "d  plane (normal sctn) (pstn sctn)" "k  plane (normal sctn) (pstn sctn)"
      using prems by auto
    ultimately show "d  B  k  B" using ¦normal sctn¦  set Basis_list
      by (auto simp: elim!: abs_in_Basis_absE)
  qed (use prems in auto elim!: abs_in_BasisE simp: eucl_le[where 'a='a] dest!: bspec subsetD)
  subgoal by simp
  subgoal using [[simproc del: defined_all]] by (auto elim!: abs_in_BasisE simp: eucl_le[where 'a='a] dest!: bspec subsetD cong del: image_cong_simp)
  subgoal using [[simproc del: defined_all]] by (auto elim!: abs_in_BasisE simp: eucl_le[where 'a='a] dest!: bspec subsetD cong del: image_cong_simp)
  done

lemma Let_unit: "Let (x::unit) f = f ()"
  by auto
lemma CHECK_no_text: "CHECKs (x#ys) a = CHECKs [] a"
  by auto

lemma frontier_above_halfspace:
  "normal sctn  0  frontier (above_halfspace sctn) = plane_of sctn"
  using frontier_halfspace_ge[of "normal sctn" "pstn sctn"]
  by (auto simp: halfspace_simps plane_of_def inner_commute)

lemma
  flowpipe_subset:
  assumes "flowpipe X0 hl hu CX X1"
    and subs: "Y0  X0" "hl  il" "il  iu" "iu  hu" "CX  CY" "X1  Y1"
    and safe: "fst ` CY  fst ` Y1  Csafe"
  shows "flowpipe Y0 il iu CY Y1"
proof -
  from assms(1) have fp: "0  hl" "hl  hu" "fst ` X0  Csafe" "fst ` CX  Csafe" "fst ` X1  Csafe"
    "(x0, d0)X0. h{hl..hu}. h  existence_ivl0 x0  (flow0 x0 h, Dflow x0 h oL d0)  X1  (h'{0..h}. (flow0 x0 h', Dflow x0 h' oL d0)  CX)"
    by (auto simp: flowpipe_def)
  then show ?thesis
    unfolding flowpipe_def
    apply safe
    subgoal using subs by auto
    subgoal using subs by auto
    subgoal using subs by fastforce
    subgoal using safe by auto
    subgoal using safe by auto
    subgoal using subs by fastforce
    subgoal using subs fp by fastforce
    subgoal for x0 d0 h h'
      using subs fp
      apply -
      apply (rule subsetD, assumption)
      apply (drule bspec)
       apply (rule subsetD; assumption)
      apply safe
      apply (drule bspec[where x=h], force)
      apply auto
      done
    done
qed

lemma poincare_mapsto_unionI:
  assumes "poincare_mapsto P r U t d"
  assumes "poincare_mapsto P s U u e"
  shows "poincare_mapsto P (r  s) U (t  u) (d  e)"
  using assms
  apply (auto simp: poincare_mapsto_def)
  subgoal
    apply (drule bspec, assumption)
    by auto
  subgoal by fastforce
  done

lemma sabove_not_le_halfspace:
  "x  sabove_halfspace sctn  ¬ le_halfspace sctn x"
  by (auto simp: sabove_halfspace_def le_halfspace_def gt_halfspace_def)

lemma (in c1_on_open_euclidean) flowsto_self:― ‹TODO: move!›
  "0  T  X0  Z  fst ` X0  X  flowsto X0 T Y Z"
  by (force simp: flowsto_def intro!: bexI[where x=0])

lemma (in c1_on_open_euclidean) flowpipe_imp_flowsto:― ‹TODO: move!›
  assumes "flowpipe X0 hl hu CX Y" "hl > 0"
  shows "flowsto X0 {0<..hl} CX Y"
  using assms
  by (fastforce simp: flowsto_def flowpipe_def open_segment_eq_real_ivl
      dest: bspec[where x=hl]
      intro!: bexI[where x=hl])

lemma flowsto_source_unionI:
  "flowsto X0 T A B  flowsto Z T A B  flowsto (X0  Z) T A B"
  by (fastforce simp: flowsto_def dest: bspec)

lemma poincare_mapsto_subset:
  "poincare_mapsto P X0 U CX1 X1  X0'  X0  X1  X2  CX1  CX2  fst ` X2  Csafe
   poincare_mapsto P X0' U CX2 X2"
  by (force simp: poincare_mapsto_def)

lemma PDP_abs_lemma:
  fixes n::"'a::executable_euclidean_space"
  assumes "abs n  Basis"
  shows
    "(x, d - (blinfun_scaleR_left (f (x)) oL (blinfun_scaleR_left (inverse (f x  n)) oL (blinfun_inner_left n oL d)))) =
     (x, d - (blinfun_scaleR_left (f (x)) oL (blinfun_scaleR_left (inverse (f x  (abs n))) oL (blinfun_inner_left (abs n) oL d))))"
proof -
  consider "n  Basis" | "- n  Basis"
    using abs_in_Basis_absE[OF assms] assms by metis
  then show ?thesis
  proof cases
    case 1
    then show ?thesis by simp
  next
    case 2
    define i where "i = -n"
    with 2 have "i  Basis" "n = -i"
      by auto
    then show ?thesis
      by (auto simp: inverse_eq_divide intro!: blinfun_eqI blinfun.bilinear_simps euclidean_eqI[where 'a='a])
  qed
qed

lemma abs_in_BasisI:
  "¦n¦  Basis" if n: "n  Basis  - n  Basis" for n::"'a::executable_euclidean_space"
proof -
  consider "n  Basis" | "- n  Basis"
    using n by auto
  then show ?thesis
  proof cases
    case 1
    then show ?thesis by simp
  next
    case 2
    define i where "i = -n"
    with 2 have "i  Basis" "n = -i"
      by auto
    then show ?thesis
      by (auto simp: inverse_eq_divide intro!: blinfun_eqI blinfun.bilinear_simps euclidean_eqI[where 'a='a])
  qed
qed

lemma flowsto_poincareD:
  assumes f: "flowsto X0 T CX X1"
  assumes X1: "fst ` X1  P"
  assumes P: "(P × UNIV)  CX = {}" "closed P"
  assumes pos: "t. t  T  t > 0"
  assumes x0: "x0  fst ` X0"
  assumes "fst ` X1  K"
  shows returns_to_flowstoI: "returns_to P x0"
    and poincare_map_mem_flowstoI: "poincare_map P x0  K"
proof -
  from x0 obtain d where x0d: "(x0, d)  X0" by auto
  from flowstoE[OF f x0d] obtain h
    where h:
      "h  T"
      "h  existence_ivl0 x0"
      "(flow0 x0 h, Dflow x0 h oL d)  X1"
      and CX: "h'. h'  {0<--<h}  (flow0 x0 h', Dflow x0 h' oL d)  CX"
    by auto
  have "h > 0" by (auto intro!: pos h)
  have "flow0 x0 h  P" using X1 h by auto
  have "F t in at_right 0. t  {0<..<h}"
    using order_tendstoD(2)[OF tendsto_ident_at 0 < h, of "{0<..}"]
    by (auto simp: eventually_at_filter)
  then have "F t in at_right 0. flow0 x0 t  fst ` CX"
    by eventually_elim (use CX 0 < h open_segment_eq_real_ivl in auto)
  then have evnP: "F t in at_right 0. flow0 x0 t  P"
    by eventually_elim (use P in force)
  from h > 0 h(2) flow0 x0 h  P evnP P(2) show "returns_to P x0"
    by (rule returns_toI)
  have nin_P: "0 < s  s < h  flow0 x0 s  P" for s
    using CX[of s] P by (auto simp: open_segment_eq_real_ivl)
  have "return_time P x0 = h"
    using h X1
    by (auto intro!: return_time_eqI 0 < h h assms simp: nin_P)
  then have "poincare_map P x0 = flow0 x0 h" by (auto simp: poincare_map_def)
  also have "  fst ` X1" using h by auto
  also note _  K
  finally show "poincare_map P x0  K" .
qed

lemma
  inner_abs_Basis_eq_zero_iff:
  "abs n  Basis  x  ¦n¦ = 0  x  n = 0" for n::"'a::executable_euclidean_space"
  by (auto simp: elim!: abs_in_BasisE)

lemmas [simp] = sbelow_halfspaces_insert

lemma Int_Un_eq_emptyI: "a  (b  c) = {}" if "a  b = {}" "a  c = {}"
  using that by auto

lemma cancel_times_UNIV_subset: "A × UNIV  B × UNIV  A  B"
  by auto

lemma split_spec_coll_spec[le,refine_vcg]:
  "split_spec_coll X  SPEC (λ(A, B). X  A  B)"
  unfolding split_spec_coll_def
  by (refine_vcg)

lemma Un_snd_sing_Pair_eq:
  "(e, f)  a  f  (xa - {(e, f)}. snd x) = (xa. snd x)"
  by force

lemma let_unit: "Let X y = y ()" by simp

lemma (in c1_on_open_euclidean) flowpipe_imp_flowsto_nonneg:― ‹TODO: move!›
  assumes "flowpipe X0 hl hu CX Y"
  shows "flowsto X0 {0..} CX Y"
  using assms
  by (fastforce simp: flowsto_def flowpipe_def open_segment_eq_real_ivl
      dest: bspec[where x=hl]
      intro!: bexI[where x=hl])

lemma subset_DiffI: "A  B  A  C = {}  A  B - C"
  by auto

lemma flowsto_source_Union: "flowsto (xR. X0 x) T CX X1"
  if "x. x  R  flowsto (X0 x) T CX X1"
  using that
  by (auto simp: flowsto_def)

lemma times_subset_iff: "A × B  C × E  A = {}  B = {}  A  C  B  E"
  by auto

lemma
  flowsto_UnionE:
  assumes "finite Gs"
  assumes "flowsto X T CX (Gs)"
  obtains XGs where "X G. (X, G)  XGs  flowsto X T CX G" "Gs = snd ` XGs" "X = (fst ` XGs)"
  apply atomize_elim
  using assms
proof (induction arbitrary: X)
  case empty
  then show ?case by auto
next
  case (insert x F)
  from insert.prems obtain X1 X2 where X: "X = X1  X2" and X1: "flowsto X1 T CX x" and X2: "flowsto X2 T CX (F)"
    by (auto elim!: flowsto_unionE)
  from insert.IH[OF X2] obtain XGs where XGs:
    "X G. (X, G)  XGs  flowsto X T CX G" "F = snd ` XGs" "X2 = (aXGs. fst a)"
    by auto
  then show ?case
    using X X1 X2
    by (intro exI[where x="insert (X1, x) XGs"]) auto
qed

lemma flowsto_Union_funE:
  assumes "finite Gs"
  assumes "flowsto X T CX (Gs)"
  obtains f where "G. G  Gs  flowsto (f G) T CX G" "X = (f ` Gs)"
  apply atomize_elim
  using assms
proof (induction arbitrary: X)
  case empty
  then show ?case by auto
next
  case (insert x F)
  from insert.prems obtain X1 X2 where X: "X = X1  X2" and X1: "flowsto X1 T CX x" and X2: "flowsto X2 T CX (F)"
    by (auto elim!: flowsto_unionE)
  from insert.IH[OF X2] obtain f where f:
    "G. G  F  flowsto (f G) T CX G" "X2 = (aF. f a)"
    by auto
  then show ?case
    using X X1 X2 insert.hyps
    by (intro exI[where x="f (x:=X1)"]) (auto split: if_splits)
qed

lemma flowsto_Union_Un_funE:
  assumes "flowsto X T CX (Gs  trap)"
  assumes "finite Gs" "Gs  {}"
  obtains f where "G. G  Gs  flowsto (f G) T CX (G  trap)" "X = (f ` Gs)"
proof -
  from assms have "flowsto X T CX (g  Gs. g  trap)" by auto
  from flowsto_Union_funE[OF finite_imageI[OF finite Gs] this]
  obtain f where f: "G. G  (λg. g  trap) ` Gs  flowsto (f G) T CX G"
    "X =  (f ` ((λg. g  trap) ` Gs))"
    by auto
  define f' where "f' g = f (g  trap)" for g
  have "G  Gs  flowsto (f' G) T CX (G  trap)" for G
    using f(1)[of "G  trap"]
    by (auto simp: f'_def)
  moreover
  have "X = (f' ` Gs)"
    using f by (auto simp: f'_def)
  ultimately show ?thesis ..
qed

lemma flowsto_Diff_to_Union_funE:
  assumes "flowsto (X - trap) T CX (Gs)"
  assumes "finite Gs"
  obtains f where "G. G  Gs  flowsto (f G - trap) T CX G" "Gs  {}  X = (f ` Gs)"
  apply atomize_elim
  using assms(2,1)
proof (induction arbitrary: X)
  case empty then show ?case by simp
next
  case (insert x F)
  from insert.prems obtain X1 X2 where X: "X - trap = X1  X2" and X1: "flowsto (X1) T CX x" and X2: "flowsto X2 T CX (F)"
    by (auto elim!: flowsto_unionE)
  then have "X1 = X1 - trap" "X2 = X2 - trap" by auto
  then have X1': "flowsto (X1 - trap) T CX x" and X2': "flowsto (X2 - trap) T CX (F)"
    using X1 X2 by auto
  from insert.IH[OF X2'] obtain f where f: "G. G  F  flowsto (f G - trap) T CX G" "F  {}  X2 = (aF. f a)"
    by auto
  show ?case
    apply (cases "F = {}")
    subgoal using f X X1 X2 X1' X2' insert.hyps insert.prems by auto
    subgoal
      apply (rule exI[where x="f (x:=X1  (trap  X))"])
      apply auto
      subgoal
        using X1
        by (rule flowsto_subset) auto
      subgoal using X X1 X2 insert.hyps by auto
      subgoal using f X X1 X2 insert.hyps by auto
      subgoal using f X X1 X2 insert.hyps by auto
      subgoal using f X X1 X2 X1' X2' insert.hyps insert.prems by auto
      subgoal using f X X1 X2 X1' X2' insert.hyps insert.prems insert by auto
      subgoal using f X X1 X2 insert.hyps by (auto split: if_splits)
      done
    done
qed

lemma refine_case_list[refine_vcg]:
  assumes "xs = []  f  SPEC P"
  assumes "y ys. xs = y # ys  g y ys  SPEC P"
  shows "(case xs of []  f | (x#xs)  g x xs)  SPEC P"
  using assms
  by (auto split: list.splits)

lemma flowsto_stays_sbelow:
  assumes "flowsto X0 {0<..} CXS X1"
  assumes "fst ` X0  below_halfspace sctn"
  assumes "x d. (x, d)  CXS  ode x  normal sctn < 0"
  shows "flowsto X0 {0<..} (CXS  sbelow_halfspace sctn × UNIV) X1"
  unfolding flowsto_def
proof safe
  fix x d assume "(x, d)  X0"
  with assms obtain t where
    "t>0" "t  existence_ivl0 x" "(s{0<..<t}. (flow0 x s, Dflow x s oL d)  CXS)"
    "(flow0 x t, Dflow x t oL d)  X1"
    by (auto simp: flowsto_def subset_iff open_segment_eq_real_ivl)
  moreover
  have "s{0<..<t}. flow0 x s  sbelow_halfspace sctn"
  proof (rule ccontr, clarsimp)
    fix s assume s: "flow0 x s  sbelow_halfspace sctn" "0 < s" "s < t"
    let ?f = "λt. flow0 x t  normal sctn - pstn sctn"
    let ?f' = "λt dx. dx * (ode (flow0 x t)  normal sctn)"
    have "xa{0<..<s}. ?f s - ?f 0 = ?f' xa (s - 0)"
      by (rule mvt[OF 0 < s, of ?f ?f'])
        (use ivl_subset_existence_ivl[OF t  existence_ivl0 x] 0 < s s < t in
          auto intro!: continuous_intros derivative_eq_intros flow_has_derivative
            simp: flowderiv_def blinfun.bilinear_simps)
    then obtain s' where "?f s - ?f 0 = ?f' s' (s - 0)" "0 < s'" "s' < s"
      by (auto simp: algebra_simps)
    note this(1)
    also
    have "(flow0 x s', Dflow x s' oL d ) CXS"
      using 0 < s' s{0<..<t}. _  CXS s < t s' < s by auto
    then have "?f' s' (s - 0) < 0"
      using assms (x, d)  X0 0 < s
      by (auto simp: flowsto_def halfspace_simps algebra_simps subset_iff intro!: mult_pos_neg)
    finally have 1: "?f s < ?f 0"
      by simp
    also have "?f 0  0"
      using (x, d)  X0 assms mem_existence_ivl_iv_defined[OF t  existence_ivl0 _]
      by (auto simp: halfspace_simps subset_iff)
    finally have "?f s < 0" .
    moreover from s have "0  ?f s"
      by (auto simp: halfspace_simps not_less)
    ultimately show False by simp
  qed
  ultimately
  show "t{0<..}. t  existence_ivl0 x  (flow0 x t, Dflow x t oL d)  X1 
                (s{0<--<t}. (flow0 x s, Dflow x s oL d)  CXS  sbelow_halfspace sctn × UNIV)"
    by (auto intro!: simp: open_segment_eq_real_ivl)
qed

lemma poincare_mapsto_Union: "poincare_mapsto P (xa) S CXS PS" 
  if "x. x  xa  poincare_mapsto P x S CXS PS"
  by (force simp: poincare_mapsto_def dest!: that)

lemma diff_subset: "(xxa. f x) - (xxa. g x)  (xxa. f x - g x)"
  by auto

lemma poincare_mapsto_imp_flowsto:
  assumes "poincare_mapsto P X0 S CX X1"
  assumes "closed P"
  shows "flowsto X0 {0<..} (CX × UNIV) (fst ` X1 × UNIV)"
  unfolding flowsto_def
proof safe
  fix x0 d0 assume "(x0, d0)  X0"
  with assms obtain D where D:
    "returns_to P x0"
    "fst ` X0  S"
    "return_time P differentiable at x0 within S"
    "(poincare_map P has_derivative blinfun_apply D) (at x0 within S)"
    "(flow0 x0 (return_time P x0), D oL d0)  X1"
    "t. 0 < t  t < return_time P x0  flow0 x0 t  CX"
    by (auto simp: poincare_mapsto_def poincare_map_def)
  show "h{0<..}.
    h  existence_ivl0 x0  (flow0 x0 h, Dflow x0 h oL d0)  fst ` X1 × UNIV 
    (h'{0<--<h}. (flow0 x0 h', Dflow x0 h' oL d0)  CX × UNIV)"
    by (auto intro!: bexI[where x="return_time P x0"] return_time_exivl D assms return_time_pos
        simp: open_segment_eq_real_ivl)
qed

lemma flowsto_poincare_mapsto_trans_flowsto:
  assumes "flowsto X0 T CX X1'"
  assumes "poincare_mapsto P X1 S CY X2"
  assumes "closed P"
  assumes "fst ` X1'  fst ` X1"
  assumes "X1'  CX  CY × UNIV  CZ"
  assumes "t. t  T  t  0"
  shows "flowsto X0 {0<..} CZ (fst ` X2 × UNIV)"
proof -
  have X1D: "(a, b)  X1'  c. (a, c)  X1" for a b using assms(4) by force
  from poincare_mapsto_imp_flowsto[OF assms(2,3)]
  have "flowsto X1 {0<..} (CY × UNIV) (fst ` X2 × UNIV)" .
  then have "flowsto X1' {0<..} (CY × UNIV) (fst ` X2 × UNIV)"
    by (auto simp: flowsto_def dest!: X1D)
  from flowsto_trans[OF assms(1) this]
  show ?thesis
    apply (rule flowsto_subset)
    using assms
    by (auto intro!: add_nonneg_pos)
qed

lemma eq_blinfun_inner_left[intro]:
  "(λx. x  n) = blinfun_apply (blinfun_inner_left n)"
  by force

lemma flowsto_union_DiffE:
  assumes "flowsto X0 T CX (Y  Z)"
  obtains X1 where "X1  X0" "flowsto X1 T CX Y" "flowsto (X0 - X1) T CX Z"
proof -
  let ?X1 = "{xX0. flowsto {x} T CX Y}"
  from assms have "?X1  X0" "flowsto ?X1 T CX Y" "flowsto (X0 - ?X1) T CX Z"
    by (auto simp: flowsto_def)
  thus ?thesis ..
qed

lemma eucl_less_le_trans:
  fixes a b::"'a::executable_euclidean_space"
  shows "eucl_less a b  b  c  eucl_less a c"
  by (force simp: eucl_less_def[where 'a='a] eucl_le[where 'a='a])

lemma le_eucl_less_trans:
  fixes a b::"'a::executable_euclidean_space"
  shows "a  b  eucl_less b c  eucl_less a c"
  by (force simp: eucl_less_def[where 'a='a] eucl_le[where 'a='a])

lemma flowsto_source_UnionI:
  assumes "i. i  I  flowsto i T CXS (f i)"
  shows "flowsto (I) T CXS ((f ` I))"
  apply (auto simp: flowsto_def)
  subgoal premises prems for y a b
    using assms[unfolded flowsto_def, OF y  I, rule_format, OF _  y] prems
    by auto
  done

lemma poincare_mapsto_UnionI:
  assumes pm[unfolded poincare_mapsto_def, rule_format]: "i. i  I  poincare_mapsto p (X0 i) S (CX i) (X1 i)"
  assumes R: "i x. i  I  x  X1 i  x  R"
  shows "poincare_mapsto p (xI. X0 x) S ((xI. CX x)) R"
  unfolding poincare_mapsto_def
proof (safe del: conjI, goal_cases)
  case (1 x0 d0 i)
  moreover
  have "fst ` (X0 ` I)  S"
    proof (safe, goal_cases)
      case (1 _ x0 d0 i)
      from this pm[OF 1]
      show ?case by auto
    qed
  ultimately show ?case using pm[OF 1]
    by (auto intro!: R)
qed

lemma tendsto_at_top_translateI:
  assumes "(f  l) (at_top::real filter)"
  shows "((λx. f (x + t)::'a::topological_space)  l) at_top"
proof (rule topological_tendstoI)
  fix S::"'a set" assume "open S" "l  S"
  from topological_tendstoD[OF assms this]
  obtain N where "n. n  N  f n  S" by (auto simp: eventually_at_top_linorder)
  then have "n. n  N - t  f (n + t)  S" by auto
  then show "F x in at_top. f (x + t)  S"
    unfolding eventually_at_top_linorder
    by blast
qed

lemma tendsto_at_top_translate_iff:
  "((λx. f (x + t)::'a::topological_space)  l) at_top  (f  l) (at_top::real filter)"
  using tendsto_at_top_translateI[of f l t]
    tendsto_at_top_translateI[of "λx. f (x + t)" l "- t"]
  by auto

lemma stable_on_mono:
  "stable_on A B" if "stable_on C B" "A  C"
  using that
  unfolding stable_on_def
  by fastforce

lemma
  flowsto_mapsto_avoid_trap:
  assumes "flowsto (X0 - trap × UNIV) {0<..} CX P"
  assumes trapprop: "stable_on (fst ` (CX  P)) trap"
  shows "flowsto (X0 - trap × UNIV) {0<..} CX (P - trap × UNIV)"
  unfolding flowsto_def
proof (rule, goal_cases)
  case (1 xd)
  from assms(1)[unfolded flowsto_def, rule_format, OF this] obtain h x0 d0 where
    "xd = (x0, d0)" "0 < h"
    "h  existence_ivl0 (x0)"
    "(flow0 x0 h, Dflow x0 h oL d0)  P"
    "(h'{0<--<h}. (flow0 x0 h', Dflow x0 h' oL d0)  CX)"
    by auto
  then show ?case
    using 1 trapprop
    apply (auto intro!: bexI[where x=h] dest!: stable_onD simp: open_segment_eq_real_ivl image_Un)
    subgoal for s by (cases "s = h") auto
    done
qed

end

lemma map_prod_def': "map_prod f g x = (f (fst x), g (snd x))"
  by (cases x) auto

lemmas rel_prod_br = br_rel_prod

lemmas lvivl_relI = lv_relivl_relI

end