Theory Example_Utilities

theory Example_Utilities
imports
  Init_ODE_Solver
begin

definition "true_form = Less (floatarith.Num 0) (floatarith.Num 1)"

lemma open_true_form[intro, simp]: "open_form true_form"
  by (auto simp: true_form_def)

lemma max_Var_form_true_form[simp]: "max_Var_form true_form = 0"
  by (auto simp: true_form_def)

lemma interpret_form_true_form[simp]: "interpret_form true_form = (λ_. True)"
  by (auto simp: true_form_def)

lemmas [simp] = length_aforms_of_ivls

declare INF_cong_simp [cong] SUP_cong_simp [cong] image_cong_simp [cong del]

declare [[ cd_patterns "_ = interpret_floatariths ?fas _" "_ = interpret_floatarith ?fa _"]]

concrete_definition reify_example for i j k uses reify_example

hide_const (open) Print.file_output
definition "file_output s f =
  (if s = STR ''''      then f (λ_. ())
  else if s = STR ''-'' then f print
  else                                  Print.file_output s f)"

definition "aforms_of_point xs = aforms_of_ivls xs xs"

definition "unit_matrix_list D = concat (map (λi. map (λj. if i = j then 1 else 0::real) [0..<D]) [0..<D])"

definition "with_unit_matrix D X = (fst X @ unit_matrix_list D, snd X @ unit_matrix_list D)"

definition "list_interval l u = {x. list_all2 (≤) l x ∧ list_all2 (≤) x u}"

context includes lifting_syntax begin
lemma list_interval_transfer[transfer_rule]:
  "((list_all2 A) ===> (list_all2 A) ===> rel_set (list_all2 A))
    list_interval list_interval"
  if [transfer_rule]: "(A ===> A ===> (=)) (≤) (≤)" "bi_total A"
  unfolding list_interval_def
  by transfer_prover
end

lemma in_list_interval_lengthD: "x ∈ list_interval a b ⟹ length x = length a"
  by (auto simp: list_interval_def list_all2_lengthD)

context includes floatarith_notation begin

definition "varvec_fas' D C = ((map Var [0..<D]) @
      concat (map (λb.
        (map (λi. (Num (C i)) +
          Var (D + D * D) * (mvmult_fa D D (map Var [D..<D + D * D]) (map Num ((replicate D 0)[i:=1])) ! b)) [0..<D])) [0..<D]))"

definition "varvec_fas D C = ((map Var [0..<D]) @
      concat (map (λi. (map (λj. (Num (C i)) + Var (D + D * D) * Var (D + D * i + j)) [0..<D])) [0..<D]))"

lemma ― ‹for illustration›
  assumes[simp]: "D=3" "rf = real_of_float"
  shows "interpret_floatariths (varvec_fas D (λi. [a, b, c] ! i))
  [a, b, c, d11, d12, d13,
            d21, d22, d23,
            d31, d32, d33, 2] =
  [rf a, rf b, rf c,
  rf a + 2 * rf d11, rf a + 2 * rf d12, rf a + 2 * rf d13,
  rf b + 2 * rf d21, rf b + 2 * rf d22, rf b + 2 * rf d23,
  rf c + 2 * rf d31, rf c + 2 * rf d32, rf c + 2 * rf d33]"
  by (simp add: varvec_fas_def mvmult_fa_def eval_nat_numeral)

definition "vareq_projections
    n  ― ‹dimension›
    ps ― ‹pairs of coordinates to project onto›
    ds ― ‹partial derivatives w.r.t. which variables›
    cs ― ‹(color) coding for partial derivatives›
  =
  [(i + n * (x + 1)::nat, i + n * (y + 1), c). (i, c) ← zip ds cs, (x, y) ← ps]"

definition "varvec_aforms_line D X line =
  approx_floatariths
    30
    (varvec_fas D (λi. float_of (fst (X ! i))))
    (take (D + D*D) X @ line)"

definition "varvec_aforms_head D X s = varvec_aforms_line D X (aforms_of_point [s])"
definition "varvec_aforms_vec D X s = varvec_aforms_line D (map (λx. (fst x, zero_pdevs)) X) [aform_of_ivl 0 s]"

definition
  "shows_aforms_vareq
      n                ― ‹dimension›
      ps               ― ‹pairs of coordinates to project onto›
      ds               ― ‹partial derivatives w.r.t. which variables›
      csl              ― ‹color coding for partial derivatives ('arrow' heads)›
      csh              ― ‹color coding for partial derivatives (lines)›
      s                ― ‹scale vectors for partial derivatives›
      (no_str::string) ― ‹default string if no C1 info is present›
      X                ― ‹affine form with C1 info›
   =
    (case (varvec_aforms_head n X s, varvec_aforms_vec n X s) of (Some X, Some Y) ⇒
        shows_sep (λ(x, y, c). shows_segments_of_aform x y X c) shows_nl (vareq_projections n ps ds csl) o shows_nl
      o shows_sep (λ(x, y, c). shows_segments_of_aform x y Y c) shows_nl (vareq_projections n ps ds csh) o shows_nl
    | _ ⇒ shows_string no_str o shows_nl)"

abbreviation "print_string s ≡ print (String.implode s)"
abbreviation "print_show s ≡ print_string (s '''')"

value [code] "print_show (shows_aforms_vareq 3 [(x, y). x ← [0..<3], y ← [0..<3], x < y]
  [0..<3] [''0x0000ff'', ''0x00ff00'', ''0xff0000''] [''0x0000ff'', ''0x00ff00'', ''0xff0000'']
  (FloatR 1 (-1)) ''# no C1 info''
    ((((λ(a, b). aforms_of_ivls a b) (with_unit_matrix 3 ([10, 20, 30], [12, 22, 32]))))))"

method_setup guess_rhs = ‹
let
  fun compute ctxt var lhs =
    let
      val lhs' = Code_Evaluation.dynamic_value_strict ctxt lhs;
      val clhs' = Thm.cterm_of ctxt lhs';
      val inst = Thm.instantiate (TVars.empty, Vars.make [(var, clhs')]);
    in PRIMITIVE inst end;
  fun eval_schematic_rhs ctxt t = (case try (HOLogic.dest_eq o HOLogic.dest_Trueprop) t of
      SOME (lhs, Var var) => compute ctxt var lhs
    | _ => no_tac);
in
  Scan.succeed (fn ctxt =>
    SIMPLE_METHOD (HEADGOAL (SUBGOAL (fn (t, _) => eval_schematic_rhs ctxt t))))
end
›

lemma length_aforms_of_point[simp]: "length (aforms_of_point xs) = length xs"
  by (auto simp: aforms_of_point_def)

definition "aform2d_plot_segments x y a = shows_segments_of_aform x y a ''0x000000''"

lemma list_of_eucl_prod[simp]: "list_of_eucl (x, y) = list_of_eucl x @ list_of_eucl y"
  by (auto simp: list_of_eucl_def Basis_list_prod_def intro!: nth_equalityI)

lemma list_of_eucl_real[simp]: "list_of_eucl (x::real) = [x]"
  by (auto simp: list_of_eucl_def Basis_list_real_def)

lemma Joints_aforms_of_ivls_self[simp]: "xs ∈ Joints (aforms_of_ivls xs xs)"
  by (auto intro!: aforms_of_ivls)

lemma Joints_aforms_of_ivls_self_eq[simp]: "Joints (aforms_of_ivls xs xs) = {xs}"
  apply (auto )
  by (auto simp: aforms_of_ivls_def Joints_def valuate_def aform_val_def
      intro!: nth_equalityI)

lemma local_lipschitz_c1_euclideanI:
  fixes T::"real set" and X::"'a::euclidean_space set"
    and f::"real ⇒ 'a ⇒ 'a"
  assumes f': "⋀t x. t ∈ T ⟹ x ∈ X ⟹ (f t has_derivative f' t x) (at x)"
  assumes cont_f': "⋀i. i ∈ Basis ⟹ continuous_on (T × X) (λ(t, x). f' t x i)"
  assumes "open T"
  assumes "open X"
  shows "local_lipschitz T X f"
  using assms
  apply (intro c1_implies_local_lipschitz[where f'="λ(t, x). Blinfun (f' t x)"])
     apply (auto simp: bounded_linear_Blinfun_apply has_derivative_bounded_linear split_beta'
      intro!: has_derivative_Blinfun continuous_on_blinfun_componentwise)
  apply (subst continuous_on_cong[OF refl]) defer apply assumption
  apply auto
  apply (subst bounded_linear_Blinfun_apply)
   apply (rule has_derivative_bounded_linear)
  by auto

definition "list_aform_of_aform (x::real aform) = (fst x, list_of_pdevs (snd x))"
primrec split_aforms_list:: "(real aform) list list
   ⇒ nat ⇒ nat ⇒ (real aform) list list" where
  "split_aforms_list Xs i 0 = Xs"
| "split_aforms_list Xs i (Suc n) = split_aforms_list (concat (map (λx. (let (a, b) = split_aforms x i in [a, b])) Xs)) i n"

definition "shows_aforms x y c X = shows_lines (map (λb. (shows_segments_of_aform x y b c ''⏎'')) X)"

end

definition "the_odo ode_fas safe_form = the(mk_ode_ops ode_fas safe_form)"

locale ode_interpretation =
  fixes safe_form::form and safe_set::"'a::executable_euclidean_space set"
    and ode_fas::"floatarith list"
    and ode::"'a ⇒ 'a"
    and finite::"'n::enum"
  assumes dims: "DIM('a) = CARD('n)"
  assumes len: "length ode_fas = CARD('n)"
  assumes safe_set_form: "safe_set = {x. interpret_form safe_form (list_of_eucl x)}"
  assumes interpret_fas: "⋀x. x ∈ safe_set ⟹ einterpret ode_fas (list_of_eucl x) = ode x"
  assumes odo: "mk_ode_ops ode_fas safe_form ≠ None"
  assumes isFDERIV: "⋀xs. interpret_form safe_form xs ⟹
    isFDERIV (length ode_fas) [0..<length ode_fas] ode_fas xs"
begin

abbreviation "odo ≡ the_odo ode_fas safe_form"
lemmas odo_def = the_odo_def

lemma odo_simps[simp]: "ode_expression odo = ode_fas" "safe_form_expr odo = safe_form"
  using odo
  by (auto simp: odo_def ode_expression_mk_ode_ops safe_form_expr_mk_ode_ops)

lemma safe_set: "safe_set = aform.Csafe odo"
  using odo dims safe_set_form isFDERIV
  unfolding aform.Csafe_def aform.safe_def aform.safe_form_def aform.ode_e_def
  by (auto simp: mk_ode_ops_def safe_set_form len split: if_splits)

lemma ode: "⋀x. x ∈ safe_set ⟹ ode x = aform.ode odo x"
  by (auto simp: aform.ode_def aform.ode_e_def interpret_fas)

sublocale auto_ll_on_open ode safe_set
  by (rule aform.auto_ll_on_open_congI[OF safe_set[symmetric] ode[symmetric]])

lemma ode_has_derivative_ode_d1: "(ode has_derivative blinfun_apply (aform.ode_d1 odo x)) (at x)"
  if "x ∈ safe_set" for x
proof -
  from aform.fderiv[OF that[unfolded safe_set]]
  have "(aform.ode odo has_derivative blinfun_apply (aform.ode_d1 odo x)) (at x)"
    by simp
  moreover
  from topological_tendstoD[OF tendsto_ident_at open_domain(2) that]
  have "∀⇩F x' in at x. x' ∈ safe_set" .
  then have "∀⇩F x' in at x. aform.ode odo x' = ode x'"
    by eventually_elim (auto simp: ode)
  ultimately show ?thesis
    by (rule has_derivative_transform_eventually) (auto simp: ode that)
qed

sublocale c1_on_open_euclidean ode "aform.ode_d1 odo" safe_set
  apply unfold_locales
  subgoal by simp
  subgoal by (simp add: ode_has_derivative_ode_d1)
  subgoal by (rule aform.cont_fderiv') (auto intro!: continuous_intros simp: safe_set)
  done

sublocale transfer_eucl_vec for a::'a and n::'n
  by unfold_locales (simp add: dims)

lemma flow_eq: "t ∈ existence_ivl0 x ⟹ aform.flow0 odo x t = flow0 x t"
  and Dflow_eq: "t ∈ existence_ivl0 x ⟹ aform.Dflow odo x t = Dflow x t"
  and ex_ivl_eq: "t ∈ aform.existence_ivl0 odo x ⟹ aform.existence_ivl0 odo x = existence_ivl0 x"
  and poincare_mapsto_eq: "closed a ⟹ aform.poincare_mapsto odo a b c d e = poincare_mapsto a b c d e"
  and flowsto_eq: "aform.flowsto odo = flowsto"
      apply -
  subgoal by (rule flow0_cong[symmetric]) (auto simp: safe_set ode)
  subgoal by (rule Dflow_cong[symmetric]) (auto simp: safe_set ode)
  subgoal by (rule existence_ivl0_cong[symmetric]) (auto simp: safe_set ode)
  subgoal
    apply (subst aform.poincare_mapsto_cong[OF safe_set[symmetric]])
    by (auto simp: ode)
  subgoal
    apply (intro ext)
    apply (subst flowsto_congI[OF safe_set ode])
    by (auto simp: safe_set)
  done

definition "avf ≡ λx::'n rvec. cast (aform.ode odo (cast x)::'a)::'n rvec"

context includes lifting_syntax begin
lemma aform_ode_transfer[transfer_rule]: "((=) ===> rel_ve ===> rel_ve) aform.ode aform.ode"
  unfolding aform.ode_def
  by transfer_prover
lemma cast_aform_ode: "cast (aform.ode odo (cast (x::'n rvec))::'a) = aform.ode odo x"
  by transfer simp

lemma aform_safe_transfer[transfer_rule]: "((=) ===> rel_ve ===> (=)) aform.safe aform.safe"
  unfolding aform.safe_def
  by transfer_prover

lemma aform_Csafe_transfer[transfer_rule]: "((=) ===> rel_set rel_ve) aform.Csafe aform.Csafe"
  unfolding aform.Csafe_def
  by transfer_prover

lemma cast_safe_set: "(cast ` safe_set::'n rvec set) = aform.Csafe odo"
  unfolding safe_set
  by transfer simp

lemma aform_ode_d_raw_transfer[transfer_rule]:
  "((=) ===> (=) ===> rel_ve ===> rel_ve ===> rel_ve ===> rel_ve) aform.ode_d_raw aform.ode_d_raw"
  unfolding aform.ode_d_raw_def
  by transfer_prover

lemma
  aform_ode_d_raw_aux_transfer:
  "((=) ===> rel_ve ===> rel_ve ===> rel_ve)
    (λx xb xa. if xb ∈ aform.Csafe x then aform.ode_d_raw x 0 xb 0 xa else 0)
    (λx xb xa. if xb ∈ aform.Csafe x then aform.ode_d_raw x 0 xb 0 xa else 0)"
  by transfer_prover

lemma aform_ode_d1_transfer[transfer_rule]:
  "((=) ===> rel_ve ===> rel_blinfun rel_ve rel_ve) aform.ode_d1 aform.ode_d1"
  apply (auto simp: rel_blinfun_def aform.ode_d1_def intro!: rel_funI)
  unfolding aform.ode_d.rep_eq
  using aform_ode_d_raw_aux_transfer
  apply -
  apply (drule rel_funD, rule refl)
  apply (drule rel_funD, assumption)
  apply (drule rel_funD; assumption)
  done

lemma cast_bl_transfer[transfer_rule]:
  "(rel_blinfun (=) (=) ===> rel_blinfun rel_ve rel_ve) id_blinfun cast_bl"
  by (auto simp: rel_ve_cast rel_blinfun_def intro!: rel_funI dest!: rel_funD)
lemma cast_bl_transfer'[transfer_rule]:
  "(rel_blinfun rel_ve rel_ve ===> rel_blinfun (=) (=)) id_blinfun cast_bl"
  apply (auto simp: rel_ve_cast rel_blinfun_def cast_cast intro!: rel_funI dest!: rel_funD)
  by (subst cast_cast) auto

lemma rel_blinfun_eq[relator_eq]: "rel_blinfun (=) (=) = (=)"
  by (auto simp: Rel_def rel_blinfun_def blinfun_ext rel_fun_eq intro!: rel_funI ext)

lemma cast_aform_ode_D1:
  "cast_bl (aform.ode_d1 odo (cast (x::'n rvec))::'a⇒⇩L'a) =
    (aform.ode_d1 odo x::'n rvec ⇒⇩L 'n rvec)"
  by transfer simp

end

definition "vf ≡ λx. cast (ode (cast x))"
definition "vf' ≡ λx::'n rvec. cast_bl (aform.ode_d1 odo (cast x::'a))
  ::'n rvec ⇒⇩L 'n rvec"
definition "vX ≡ cast ` safe_set"
sublocale a?: transfer_c1_on_open_euclidean a n ode "aform.ode_d1 odo" safe_set vf vf' vX
  for a::'a and n::'n
  by unfold_locales
    (simp_all add: dims vf_def vf'_def vX_def)

sublocale av: transfer_c1_on_open_euclidean a n "aform.ode odo" "aform.ode_d1 odo"
  "(aform.Csafe odo)" avf vf' vX
  for a::'a and n::'n
     apply unfold_locales
  unfolding vX_def
  by (simp_all add: dims avf_def  safe_set)

lemma vflow_eq: "t ∈ v.existence_ivl0 x ⟹ aform.flow0 odo x t = v.flow0 x t"
  thm flow_eq[of t "cast x"] flow_eq[of t "cast x", untransferred]
  apply (subst flow_eq[of t "cast x", untransferred, symmetric])
   apply simp
  unfolding avf_def vX_def cast_aform_ode cast_safe_set
  ..

lemma vf'_eq: "vf' = aform.ode_d1 odo"
  unfolding vf'_def cast_aform_ode_D1 ..

lemma vDflow_eq: "t ∈ v.existence_ivl0 x ⟹ aform.Dflow odo x t = v.Dflow x t"
  apply (subst Dflow_eq[of t "cast x", untransferred, symmetric])
   apply simp
  unfolding avf_def vX_def cast_aform_ode cast_safe_set vf'_eq
  ..

lemma vex_ivl_eq: "t ∈ aform.existence_ivl0 odo x ⟹ aform.existence_ivl0 odo x = v.existence_ivl0 x"
  apply (subst ex_ivl_eq[of t "cast x", untransferred, symmetric])
  unfolding avf_def vX_def cast_aform_ode cast_safe_set vf'_eq
  by auto

context includes lifting_syntax begin
lemma id_cast_eucl1_transfer_eq: "(λx. x) = (λx. (fst x, 1⇩L o⇩L snd x o⇩L 1⇩L))"
  by auto
lemma cast_eucl1_transfer[transfer_rule]:
  "(rel_prod (=) (rel_blinfun (=) (=)) ===> rel_prod rel_ve (rel_blinfun rel_ve rel_ve)) (λx. x) cast_eucl1"
  unfolding cast_eucl1_def id_cast_eucl1_transfer_eq
  apply transfer_prover_start
       apply (transfer_step)
      apply (transfer_step)
     apply (transfer_step)
    apply (transfer_step)
   apply (transfer_step)
  apply simp
  done

end

lemma avpoincare_mapsto_eq:
  "aform.poincare_mapsto odo a (b::'n eucl1 set) c d e = av.v.poincare_mapsto a b c d e"
  if "closed a"
  unfolding avf_def vX_def cast_aform_ode cast_safe_set vf'_eq
  by auto

lemma vpoincare_mapsto_eq:
  "aform.poincare_mapsto odo a (b::'n eucl1 set) c d e = v.poincare_mapsto a b c d e"
  if "closed a"
proof -
  have "closed (cast ` a::'a set)" using that
    by transfer auto
  from poincare_mapsto_eq[of "cast ` a::'a set"
      "cast_eucl1 ` b::('a × 'a ⇒⇩L 'a) set"
      "cast ` c::'a set" "cast ` d::'a set" "cast_eucl1 ` e::('a × 'a ⇒⇩L 'a) set", OF this, untransferred]
  have "v.poincare_mapsto a b c d e = av.v.poincare_mapsto a b c d e"
    by auto
  also have "… = aform.poincare_mapsto odo a (b::'n eucl1 set) c d e"
    unfolding avf_def vX_def cast_aform_ode cast_safe_set vf'_eq
    by auto
  finally show ?thesis by simp
qed

lemma avflowsto_eq: "aform.flowsto odo = (av.v.flowsto::'n eucl1 set ⇒ _)"
proof (intro ext, goal_cases)
  case (1 a b c d)
  have "av.v.flowsto a b c d = aform.flowsto odo a b c d"
    unfolding avf_def vX_def cast_aform_ode cast_safe_set vf'_eq
    by auto
  then show ?case by simp
qed

lemma vflowsto_eq: "aform.flowsto odo = (v.flowsto::'n eucl1 set ⇒ _)"
proof (intro ext, goal_cases)
  case (1 a b c d)
  have "aform.flowsto odo (cast_eucl1 ` a::'a c1_info set) b
    (cast_eucl1 ` c)  (cast_eucl1 ` d) =
    flowsto (cast_eucl1 ` a::'a c1_info set) b (cast_eucl1 ` c)  (cast_eucl1 ` d)"
    by (subst flowsto_eq) auto
  from this[untransferred] have "v.flowsto a b c d = av.v.flowsto a b c d" by auto
  also have "… = aform.flowsto odo a b c d"
    unfolding avf_def vX_def cast_aform_ode cast_safe_set vf'_eq
    by auto
  finally show ?case by simp
qed

context includes lifting_syntax begin
lemma flow1_of_list_transfer[transfer_rule]:
  "(list_all2 (=) ===> rel_prod rel_ve (rel_blinfun rel_ve rel_ve))
   flow1_of_list flow1_of_list"
  unfolding flow1_of_list_def blinfun_of_list_def o_def flow1_of_vec1_def
  by transfer_prover

lemma c1_info_of_appr_transfer[transfer_rule]:
  "(rel_prod (list_all2 (=)) (rel_option (list_all2 (=))) ===> rel_set (rel_prod rel_ve (rel_blinfun rel_ve rel_ve)))
    aform.c1_info_of_appr
    aform.c1_info_of_appr"
  unfolding aform.c1_info_of_appr_def
  by transfer_prover

lemma c0_info_of_appr_transfer[transfer_rule]:
  "((list_all2 (=)) ===> rel_set rel_ve) aform.c0_info_of_appr aform.c0_info_of_appr"
  unfolding aform.c0_info_of_appr_def
  by transfer_prover

lemma aform_scaleR2_transfer[transfer_rule]:
  "((=) ===> (=) ===> rel_set (rel_prod A B) ===> rel_set (rel_prod A B))
    scaleR2 scaleR2"
  if [unfolded Rel_def, transfer_rule]: "((=) ===> B ===> B) (*⇩R) (*⇩R)"
  unfolding scaleR2_def
  by transfer_prover
lemma scaleR_rel_blinfun_transfer[transfer_rule]: "((=) ===> rel_blinfun rel_ve rel_ve ===> rel_blinfun rel_ve rel_ve) (*⇩R) (*⇩R)"
  apply (auto intro!: rel_funI simp: rel_blinfun_def blinfun.bilinear_simps)
  apply (drule rel_funD)
   apply assumption
  apply (rule scaleR_transfer[THEN rel_funD, THEN rel_funD])
   apply auto
  done
lemma c1_info_of_appre_transfer[transfer_rule]:
  "(rel_prod (rel_prod (=) (=)) (rel_prod (list_all2 (=)) (rel_option (list_all2 (=)))) ===>
      rel_set (rel_prod rel_ve (rel_blinfun rel_ve rel_ve)))
    aform.c1_info_of_appre
    aform.c1_info_of_appre"
  unfolding aform.c1_info_of_appre_def
  by transfer_prover

lemma c1_info_of_apprs_transfer[transfer_rule]:
  "((=) ===> rel_set (rel_prod rel_ve (rel_blinfun rel_ve rel_ve)))
    aform.c1_info_of_apprs
    aform.c1_info_of_apprs"
  unfolding aform.c1_info_of_apprs_def
  by transfer_prover
lemma c1_info_of_appr'_transfer[transfer_rule]:
  "(rel_option (list_all2 (=)) ===> rel_set (rel_prod rel_ve (rel_blinfun rel_ve rel_ve)))
    aform.c1_info_of_appr' aform.c1_info_of_appr'"
  unfolding aform.c1_info_of_appr'_def
  by transfer_prover

lemma c0_info_of_apprs_transfer[transfer_rule]:
  "((=) ===> rel_set rel_ve)
    aform.c0_info_of_apprs
    aform.c0_info_of_apprs"
  unfolding aform.c0_info_of_apprs_def
  by transfer_prover
lemma c0_info_of_appr'_transfer[transfer_rule]:
  "(rel_option (list_all2 (=)) ===> rel_set rel_ve)
    aform.c0_info_of_appr' aform.c0_info_of_appr'"
  unfolding aform.c0_info_of_appr'_def
  by transfer_prover

lemma aform_Csafe_vX[simp]: "aform.Csafe odo = (vX::'n rvec set)"
  by (simp add: vX_def cast_safe_set)

definition blinfuns_of_lvivl::"real list × real list ⇒ ('b ⇒⇩L 'b::executable_euclidean_space) set"
  where "blinfuns_of_lvivl x = blinfun_of_list ` list_interval (fst x) (snd x)"

lemma blinfun_of_list_transfer[transfer_rule]:
  "(list_all2 (=) ===> rel_blinfun rel_ve rel_ve) blinfun_of_list blinfun_of_list"
  unfolding blinfun_of_list_def
  by transfer_prover

lemma blinfuns_of_lvivl_transfer[transfer_rule]:
  "(rel_prod (list_all2 (=)) (list_all2 (=)) ===> rel_set (rel_blinfun rel_ve rel_ve))
    blinfuns_of_lvivl
    blinfuns_of_lvivl"
  unfolding blinfuns_of_lvivl_def
  by transfer_prover

definition "blinfuns_of_lvivl' x = (case x of None ⇒ UNIV | Some x ⇒ blinfuns_of_lvivl x)"

lemma blinfuns_of_lvivl'_transfer[transfer_rule]:
  "(rel_option (rel_prod (list_all2 (=)) (list_all2 (=))) ===> rel_set (rel_blinfun rel_ve rel_ve))
    blinfuns_of_lvivl'
    blinfuns_of_lvivl'"
  unfolding blinfuns_of_lvivl'_def
  by transfer_prover


lemma atLeastAtMost_transfer[transfer_rule]:
  "(A ===> A ===> rel_set A) atLeastAtMost atLeastAtMost"
  if [transfer_rule]: "(A ===> A ===> (=)) (≤) (≤)" "bi_total A" "bi_unique A"
  unfolding atLeastAtMost_def atLeast_def atMost_def
  by transfer_prover

lemma set_of_ivl_transfer[transfer_rule]:
  "(rel_prod A A ===> rel_set A) set_of_ivl set_of_ivl"
  if [transfer_rule]: "(A ===> A ===> (=)) (≤) (≤)" "bi_total A" "bi_unique A"
  unfolding set_of_ivl_def
  by transfer_prover

lemma set_of_lvivl_transfer[transfer_rule]:
  "(rel_prod (list_all2 (=)) (list_all2 (=)) ===> rel_set rel_ve) set_of_lvivl set_of_lvivl"
  unfolding set_of_lvivl_def
  by transfer_prover

lemma set_of_lvivl_eq: "set_of_lvivl I =
    (eucl_of_list ` list_interval (fst I) (snd I)::'b::executable_euclidean_space set)"
  if [simp]: "length (fst I) = DIM('b)" "length (snd I) = DIM('b)"
proof (auto simp: set_of_lvivl_def list_interval_def set_of_ivl_def, goal_cases)
  case (1 x)
  with lv_rel_le[where 'a='b, param_fo, OF lv_relI lv_relI, of  "fst I" "list_of_eucl x"]
    lv_rel_le[where 'a='b, param_fo, OF lv_relI lv_relI, of "list_of_eucl x" "snd I"]
  show ?case by force
next
  case (2 x)
  with lv_rel_le[where 'a='b, param_fo, OF lv_relI lv_relI, of  "fst I" "x"]
  show ?case by (auto simp: list_all2_lengthD)
next
  case (3 x)
  with lv_rel_le[where 'a='b, param_fo, OF lv_relI lv_relI, of "x" "snd I"]
  show ?case by (auto simp: list_all2_lengthD)
qed

lemma bounded_linear_matrix_vector_mul[THEN bounded_linear_compose, bounded_linear_intros]:
  "bounded_linear ((*v) x)" for x::"real^'x^'y"
  unfolding linear_linear
  by (rule matrix_vector_mul_linear)

lemma blinfun_of_list_eq: "blinfun_of_list x = blinfun_of_vmatrix (eucl_of_list x::((real, 'b) vec, 'b) vec)"
  if "length x = CARD('b::enum)*CARD('b)"
  unfolding blinfun_of_list_def
  apply (transfer fixing: x)
  apply (rule linear_eq_stdbasis)
  unfolding linear_conv_bounded_linear
    apply (auto intro!: bounded_linear_intros)
proof goal_cases
  case (1 b)
  have "(eucl_of_list x::((real, 'b) vec, 'b) vec) *v b = (eucl_of_list x::((real, 'b) vec, 'b) vec) *v eucl_of_list (list_of_eucl b)"
    by simp
  also have "… = (∑i<CARD('b).
        (∑j<CARD('b). x ! (i * CARD('b) + j) * (b ∙ Basis_list ! j)) *⇩R Basis_list ! i)"
    by (subst eucl_of_list_matrix_vector_mult_eq_sum_nth_Basis_list)
      (auto simp: that)
  also have "… = (∑i∈Basis.
       ∑j∈Basis. (b ∙ j * x ! (index Basis_list i * CARD('b) + index Basis_list j)) *⇩R i)"
    apply (subst sum_list_Basis_list[symmetric])+
    apply (subst sum_list_sum_nth)+
    by (auto simp add: atLeast0LessThan scaleR_sum_left intro!: sum.cong)
  finally show ?case by simp
qed

lemma blinfuns_of_lvivl_eq: "blinfuns_of_lvivl x =
    (blinfun_of_vmatrix ` set_of_lvivl x::((real, 'b) vec ⇒⇩L (real, 'b) vec) set)"
  if "length (fst x) = CARD('b::enum)*CARD('b)" "length (snd x) = CARD('b)*CARD('b)"
  apply (subst set_of_lvivl_eq)
  subgoal by (simp add: that)
  subgoal by (simp add: that)
  unfolding blinfuns_of_lvivl_def image_image
  by (auto simp: that blinfun_of_list_eq[symmetric] in_list_interval_lengthD cong: image_cong)

lemma range_blinfun_of_vmatrix[simp]: "range blinfun_of_vmatrix = UNIV"
  apply auto
  apply transfer
  subgoal for x by (rule image_eqI[where x="matrix x"]) auto
  done

lemma blinfun_of_vmatrix_image:
  "blinfun_of_vmatrix ` aform.set_of_lvivl' x =
    (blinfuns_of_lvivl' x::((real, 'b) vec ⇒⇩L (real, 'b) vec) set)"
  if "aform.lvivl'_invar (CARD('b)*CARD('b::enum)) x"
  using that
  by (auto simp: aform.set_of_lvivl'_def blinfuns_of_lvivl'_def blinfuns_of_lvivl_eq 
    aform.lvivl'_invar_def
      split: option.splits)

lemma one_step_result123:
  "solves_one_step_until_time_aform optns odo X0i t1 t2 E dE ⟹
    (x0, d0) ∈ aform.c1_info_of_appre X0i ⟹
    t ∈ {t1 .. t2} ⟹
    set_of_lvivl E ⊆ S ⟹
    blinfuns_of_lvivl' dE ⊆ dS ⟹
    length (fst E) = CARD('n) ⟹ length (snd E) = CARD('n) ⟹
    aform.lvivl'_invar (CARD('n) * CARD('n)) dE ⟹
    aform.c1_info_invare DIM('a) X0i ⟹
    aform.D odo = DIM('a) ⟹
      (t ∈ existence_ivl0 (x0::'a) ∧ flow0 x0 t ∈ S) ∧ Dflow x0 t o⇩L d0 ∈ dS"
  apply (transfer fixing: optns X0i t1 t2 t E dE)
  subgoal premises prems for x0 d0 S dS
  proof -
    have "t ∈ aform.existence_ivl0 odo x0 ∧ aform.flow0 odo x0 t ∈ S ∧ aform.Dflow odo x0 t o⇩L d0 ∈ dS"
      apply (rule one_step_in_ivl[of t t1 t2 x0 d0 X0i "fst E" "snd E" S dE dS odo optns])
      using prems
      by (auto simp: eucl_of_list_prod set_of_lvivl_def set_of_ivl_def blinfun_of_vmatrix_image aform.D_def
          solves_one_step_until_time_aform_def)
    with vflow_eq[of t x0] vDflow_eq[of t x0] vex_ivl_eq[symmetric, of t x0] 
    show ?thesis
      by simp
  qed
  done

lemmas one_step_result12 = one_step_result123[THEN conjunct1]
  and one_step_result3 = one_step_result123[THEN conjunct2]
lemmas one_step_result1 = one_step_result12[THEN conjunct1]
  and one_step_result2 = one_step_result12[THEN conjunct2]

lemma plane_of_transfer[transfer_rule]: "(rel_sctn A ===> rel_set A) plane_of plane_of"
  if [transfer_rule]: "(A ===> A ===> (=)) (∙) (∙)" "bi_total A"
  unfolding plane_of_def
  by transfer_prover

lemma below_halfspace_transfer[transfer_rule]: "(rel_sctn A ===> rel_set A) below_halfspace below_halfspace"
  if [transfer_rule]: "(A ===> A ===> (=)) (∙) (∙)" "bi_total A"
  unfolding below_halfspace_def le_halfspace_def
  by transfer_prover

definition "rel_nres A a b ⟷ (a, b) ∈ ⟨{(a, b). A a b}⟩nres_rel"

lemma FAILi_transfer[transfer_rule]: "(rel_nres B) FAILi FAILi"
  by (auto simp: rel_nres_def nres_rel_def)

lemma RES_transfer[transfer_rule]: "(rel_set B ===> rel_nres B) RES RES"
  by (auto simp: rel_nres_def nres_rel_def rel_set_def intro!: rel_funI RES_refine)

context includes autoref_syntax begin
lemma RETURN_dres_nres_relI:
  "(fi, f) ∈ A → B ⟹ (λx. dRETURN (fi x), (λx. RETURN (f x))) ∈ A → ⟨B⟩dres_nres_rel"
  by (auto simp: dres_nres_rel_def dest: fun_relD)
end

lemma br_transfer[transfer_rule]:
  "((B ===> C) ===> (B ===> (=)) ===> rel_set (rel_prod B C)) br br"
  if [transfer_rule]: "bi_total B"  "bi_unique C" "bi_total C"
  unfolding br_def
  by transfer_prover

lemma aform_appr_rel_transfer[transfer_rule]:
  "(rel_set (rel_prod (list_all2 (=)) (rel_set rel_ve))) aform.appr_rel aform.appr_rel"
  unfolding aform.appr_rel_br
  by (transfer_prover)

lemma appr1_rel_transfer[transfer_rule]: "(rel_set (rel_prod
  (rel_prod (list_all2 (=)) (rel_option (list_all2 (=))))
  (rel_set (rel_prod rel_ve (rel_blinfun rel_ve rel_ve))))) aform.appr1_rel aform.appr1_rel"
  unfolding aform.appr1_rel_internal
  by transfer_prover

lemma relAPP_transfer[transfer_rule]:
  "((rel_set (rel_prod B C) ===> D) ===> rel_set (rel_prod B C) ===> D) Relators.relAPP Relators.relAPP"
  unfolding relAPP_def
  by transfer_prover

lemma prod_rel_transfer[transfer_rule]:
  "(rel_set (rel_prod B C) ===> rel_set (rel_prod D E) ===> rel_set (rel_prod (rel_prod B D) (rel_prod C E)))
    prod_rel prod_rel"
  if [transfer_rule]:
  "bi_total B" "bi_unique B"
  "bi_unique C" "bi_total C"
  "bi_unique D" "bi_total D"
  "bi_unique E" "bi_total E"
  unfolding prod_rel_def_internal
  by transfer_prover

lemma Domain_transfer[transfer_rule]:
  "(rel_set (rel_prod A B) ===> rel_set A) Domain Domain"
  if [transfer_rule]:
  "bi_total A" "bi_unique A"
  "bi_total B" "bi_unique B"
  unfolding Domain_unfold
  by transfer_prover

lemma set_rel_transfer[transfer_rule]:
  "(rel_set (rel_prod B C) ===> rel_set (rel_prod (rel_set B) (rel_set C))) set_rel set_rel"
  if [transfer_rule]:
  "bi_total B" "bi_unique B"
  "bi_unique C" "bi_total C"
  unfolding set_rel_def_internal
  by transfer_prover

lemma relcomp_transfer[transfer_rule]:
  "(rel_set (rel_prod B C) ===> rel_set (rel_prod C D) ===> rel_set (rel_prod B D)) relcomp relcomp"
  if [transfer_rule]:
  "bi_total B" "bi_unique B"
  "bi_unique C" "bi_total C"
  "bi_unique D" "bi_total D"
  unfolding relcomp_unfold
  by transfer_prover

lemma Union_rel_transfer[transfer_rule]:
  "(rel_set (rel_prod B (rel_set C)) ===> rel_set (rel_prod C (rel_set D)) ===> rel_set (rel_prod B (rel_set D)))
    Union_rel Union_rel"
  if [transfer_rule]:
  "bi_total B" "bi_unique B"
  "bi_unique C" "bi_total C"
  "bi_unique D" "bi_total D"
  unfolding Union_rel_internal top_fun_def top_bool_def
  by transfer_prover

lemma fun_rel_transfer[transfer_rule]:
  "(rel_set (rel_prod B C) ===> rel_set (rel_prod D E) ===> rel_set (rel_prod (B ===> D) (C ===> E))) Relators.fun_rel Relators.fun_rel"
  if [transfer_rule]:
  "bi_unique B"
  "bi_unique C"
  "bi_unique D" "bi_total D"
  "bi_unique E" "bi_total E"
  unfolding fun_rel_def_internal
  by transfer_prover

lemma c1_info_of_apprse_transfer[transfer_rule]:
  "(list_all2 (rel_prod (rel_prod (=) (=)) (rel_prod (list_all2 (=)) (rel_option (list_all2 (=)))))
    ===> rel_set (rel_prod rel_ve (rel_blinfun rel_ve rel_ve)))
    aform.c1_info_of_apprse
    aform.c1_info_of_apprse"
  unfolding aform.c1_info_of_apprse_def
  by transfer_prover

term scaleR2_rel
(*
"scaleR2_rel"
  ::
  "('b × ('c × 'd) set) set
   ⇒ (((ereal × ereal) × 'b) × ('c × 'd) set) set"
*)
lemma scaleR2_rel_transfer[transfer_rule]:
  "(rel_set (rel_prod (=) (rel_set (rel_prod (=) (=)))) ===>
    rel_set (rel_prod (rel_prod (rel_prod (=) (=)) (=)) (rel_set (rel_prod (=) (=))))) scaleR2_rel scaleR2_rel"
  unfolding scaleR2_rel_internal
  by transfer_prover

lemma appr1_rele_transfer[transfer_rule]:
  "(rel_set (rel_prod
    (rel_prod (rel_prod (=) (=)) (rel_prod (list_all2 (=)) (rel_option (list_all2 (=)))))
    (rel_set (rel_prod rel_ve (rel_blinfun rel_ve rel_ve))))) aform.appr1e_rel aform.appr1e_rel"
  unfolding scaleR2_rel_internal
  by transfer_prover

lemma flow1_of_vec1_times: "flow1_of_vec1 ` (A × B) = A × blinfun_of_vmatrix ` B"
  by (auto simp: flow1_of_vec1_def vec1_of_flow1_def)

lemma stable_on_transfer[transfer_rule]:
  "(rel_set rel_ve ===> rel_set rel_ve ===> (=)) v.stable_on stable_on"
  unfolding stable_on_def v.stable_on_def
  by transfer_prover

theorem solves_poincare_map_aform:
  "solves_poincare_map_aform optns odo (λx. dRETURN (symstart x)) [S] guards ivl sctn roi XS RET dRET ⟹
    (symstart, symstarta) ∈ fun_rel (aform.appr1e_rel) (clw_rel aform.appr_rel ×⇩r clw_rel aform.appr1e_rel) ⟹
    (⋀X0. (λ(CX, X). flowsto (X0 - trap × UNIV) {0..} (CX × UNIV) X) (symstarta X0)) ⟹
    stable_on (aform.Csafe odo - set_of_lvivl ivl ∩ plane_of (map_sctn eucl_of_list sctn)) trap ⟹
    (⋀X. X ∈ set XS ⟹ aform.c1_info_invare DIM('a) X) ⟹
    aform.D odo = DIM('a) ⟹
    length (normal sctn) = DIM('a) ⟹
    length (fst ivl) = DIM('a) ⟹
    length (snd ivl) = DIM('a) ⟹
    length (normal S) = DIM('a) ⟹
    (⋀a xs b ba ro.
        (xs, ro) ∈ set guards ⟹
        ((a, b), ba) ∈ set xs ⟹
        length a = DIM('a) ∧
        length b = DIM('a) ∧ length (normal ba) = DIM('a)) ⟹
    length (fst RET) = CARD('n) ⟹ length (snd RET) = CARD('n) ⟹
    aform.lvivl'_invar (CARD('n) * CARD('n)) dRET ⟹
    poincare_mapsto
     ((set_of_lvivl ivl::('a set)) ∩ plane_of (map_sctn eucl_of_list sctn))
     (aform.c1_info_of_apprse XS - trap × UNIV)
     (below_halfspace (map_sctn eucl_of_list S))
     (aform.Csafe odo -
      set_of_lvivl ivl ∩ plane_of (map_sctn eucl_of_list sctn))
     (set_of_lvivl RET × blinfuns_of_lvivl' dRET)"
  apply (transfer fixing: optns symstart S guards ivl sctn roi XS RET dRET)
  subgoal premises prems for symstarta trap
  proof -
    have "aform.poincare_mapsto odo (set_of_lvivl ivl ∩ plane_of (map_sctn eucl_of_list sctn))
     (aform.c1_info_of_apprse XS - trap × UNIV) (below_halfspace (map_sctn eucl_of_list S))
     (aform.Csafe odo - set_of_lvivl ivl ∩ plane_of (map_sctn eucl_of_list sctn))
     (flow1_of_vec1 ` ({eucl_of_list (fst RET)..eucl_of_list (snd RET)} × aform.set_of_lvivl' dRET))"
      apply (rule solves_poincare_map[OF _ RETURN_dres_nres_relI RETURN_rule,
        of optns odo symstart S guards ivl sctn roi XS "fst RET" "snd RET" dRET symstarta trap])
      subgoal using prems(1) by (simp add: solves_poincare_map_aform_def)
      subgoal using prems(2) by (auto simp: fun_rel_def_internal)
      subgoal for X0
        using prems(3)[of X0] vflowsto_eq
        by auto
      subgoal
        unfolding aform.stable_on_def
      proof (safe, goal_cases)
        case (1 t x0)
        from 1 have a: "t ∈ v.existence_ivl0 x0" using vex_ivl_eq by blast
        with 1 have b: "v.flow0 x0 t ∈ trap" using vflow_eq by simp
        have c: "v.flow0 x0 s ∈ vX - set_of_lvivl ivl ∩ plane_of (map_sctn eucl_of_list sctn)"
          if s: "s ∈ {0<..t}"
          for s
          using 1(4)[rule_format, OF s]
          apply (subst (asm) vflow_eq)
          unfolding aform_Csafe_vX[symmetric]
          using s a
          by (auto dest!: a.v.ivl_subset_existence_ivl)
        from a b c show ?case 
          using prems(4)[unfolded v.stable_on_def, rule_format, OF b a 1(3) c]
          by simp
      qed
      subgoal using prems by auto
      subgoal using prems by (auto simp: aform.D_def)
      subgoal using prems by auto
      subgoal using prems by auto
      subgoal using prems by auto
      subgoal using prems by auto
      subgoal using prems by auto
      subgoal using prems by auto
      subgoal using prems by auto
      subgoal using prems by auto
      done
    then show ?thesis
      using vflow_eq vex_ivl_eq vflowsto_eq prems
      apply (subst vpoincare_mapsto_eq[symmetric])
      by (auto simp: set_of_lvivl_def set_of_ivl_def blinfun_of_vmatrix_image flow1_of_vec1_times)
  qed
  done

theorem solves_poincare_map_aform':
  "solves_poincare_map_aform' optns odo S guards ivl sctn roi XS RET dRET⟹
    (⋀X. X ∈ set XS ⟹ aform.c1_info_invare DIM('a) X) ⟹
    aform.D odo = DIM('a) ⟹
    length (normal sctn) = DIM('a) ⟹
    length (fst ivl) = DIM('a) ⟹
    length (snd ivl) = DIM('a) ⟹
    length (normal S) = DIM('a) ⟹
    (⋀a xs b ba ro.
        (xs, ro) ∈ set guards ⟹
        ((a, b), ba) ∈ set xs ⟹
        length a = DIM('a) ∧
        length b = DIM('a) ∧ length (normal ba) = DIM('a)) ⟹
    length (fst RET) = CARD('n) ⟹ length (snd RET) = CARD('n) ⟹
    aform.lvivl'_invar (CARD('n) * CARD('n)) dRET ⟹
    poincare_mapsto
     ((set_of_lvivl ivl::('a set)) ∩ plane_of (map_sctn eucl_of_list sctn))
     (aform.c1_info_of_apprse XS)
     (below_halfspace (map_sctn eucl_of_list S))
     (aform.Csafe odo -
      set_of_lvivl ivl ∩ plane_of (map_sctn eucl_of_list sctn))
     (set_of_lvivl RET × blinfuns_of_lvivl' dRET)"
  apply (transfer fixing: optns S guards ivl sctn roi XS RET dRET)
  subgoal
    using solves_poincare_map'[of optns odo S guards ivl sctn roi XS "fst RET" "snd RET" dRET]
    using vflow_eq vex_ivl_eq vflowsto_eq
    apply (subst vpoincare_mapsto_eq[symmetric])
    by (auto intro!: closed_Int simp: set_of_lvivl_def set_of_ivl_def blinfun_of_vmatrix_image
        flow1_of_vec1_times aform.D_def solves_poincare_map_aform'_def)
  done

theorem solves_poincare_map_onto_aform:
  "solves_poincare_map_onto_aform optns odo guards ivl sctn roi XS RET dRET⟹
    (⋀X. X ∈ set XS ⟹ aform.c1_info_invare DIM('a) X) ⟹
    aform.D odo = DIM('a) ⟹
    length (normal sctn) = DIM('a) ⟹
    length (fst ivl) = DIM('a) ⟹
    length (snd ivl) = DIM('a) ⟹
    (⋀a xs b ba ro.
        (xs, ro) ∈ set guards ⟹
        ((a, b), ba) ∈ set xs ⟹
        length a = DIM('a) ∧
        length b = DIM('a) ∧ length (normal ba) = DIM('a)) ⟹
    length (fst RET) = CARD('n) ⟹ length (snd RET) = CARD('n) ⟹
    aform.lvivl'_invar (CARD('n) * CARD('n)) dRET ⟹
    poincare_mapsto
     ((set_of_lvivl ivl::('a set)) ∩ plane_of (map_sctn eucl_of_list sctn))
     (aform.c1_info_of_apprse XS)
     UNIV
     (aform.Csafe odo -
      set_of_lvivl ivl ∩ plane_of (map_sctn eucl_of_list sctn))
     (set_of_lvivl RET × blinfuns_of_lvivl' dRET)"
  apply (transfer fixing: optns guards ivl sctn roi XS RET dRET)
  subgoal
    using solves_poincare_map_onto[of optns odo guards ivl sctn roi XS "fst RET" "snd RET" dRET, where 'n='n,
          unfolded aform.poincare_maps_onto_def]
    using vflow_eq vex_ivl_eq vflowsto_eq
    apply (subst vpoincare_mapsto_eq[symmetric])
    by (auto intro!: closed_Int simp: set_of_lvivl_def set_of_ivl_def blinfun_of_vmatrix_image
        flow1_of_vec1_times aform.D_def solves_poincare_map_onto_aform_def)
  done

end

end

subsection ‹Example Utilities!›

hide_const floatarith.Max floatarith.Min

lemma degree_sum_pdevs_scaleR_Basis:
  "degree (sum_pdevs (λi. pdevs_scaleR (a i) i) (Basis::'b::euclidean_space set)) = Max ((λi. degree (a i)) ` Basis)"
  apply (rule antisym)
  subgoal apply (rule degree_le)
    by (auto )
  subgoal
    apply (rule Max.boundedI)
      apply simp
     apply simp
    apply (auto simp: intro!: degree_leI)
    by (auto simp: euclidean_eq_iff[where 'a='b])
  done

lemma Inf_aform_eucl_of_list_aform:
  assumes "length a = DIM('b::executable_euclidean_space)"
  shows "Inf_aform (eucl_of_list_aform a::'b aform) = eucl_of_list (map Inf_aform a)"
  using assms
  apply (auto simp: eucl_of_list_aform_def Inf_aform_def[abs_def] algebra_simps
      eucl_of_list_inner inner_sum_left
      intro!: euclidean_eqI[where 'a='b])
  apply (auto simp: tdev_def inner_sum_left abs_inner inner_Basis if_distrib cong: if_cong)
  apply (rule sum.mono_neutral_cong_left)
     apply simp
  by (auto simp: degree_sum_pdevs_scaleR_Basis)

lemma Sup_aform_eucl_of_list_aform:
  assumes "length a = DIM('b::executable_euclidean_space)"
  shows "Sup_aform (eucl_of_list_aform a::'b aform) = eucl_of_list (map Sup_aform a)"
  using assms
  apply (auto simp: eucl_of_list_aform_def Sup_aform_def[abs_def] algebra_simps
      eucl_of_list_inner inner_sum_left
      intro!: euclidean_eqI[where 'a='b])
  apply (auto simp: tdev_def inner_sum_left abs_inner inner_Basis if_distrib cong: if_cong)
  apply (rule sum.mono_neutral_cong_right)
     apply simp
  by (auto simp: degree_sum_pdevs_scaleR_Basis)

lemma
  eucl_of_list_map_Inf_aform_leI:
  assumes "x ∈ Affine (eucl_of_list_aform a::'b::executable_euclidean_space aform)"
  assumes "length a = DIM('b)"
  shows "eucl_of_list (map Inf_aform a) ≤ x"
  using Inf_aform_le_Affine[OF assms(1)] assms(2)
  by (auto simp: Inf_aform_eucl_of_list_aform)

lemma
  eucl_of_list_map_Sup_aform_geI:
  assumes "x ∈ Affine (eucl_of_list_aform a::'b::executable_euclidean_space aform)"
  assumes "length a = DIM('b)"
  shows "x ≤ eucl_of_list (map Sup_aform a)"
  using Sup_aform_ge_Affine[OF assms(1)] assms(2)
  by (auto simp: Sup_aform_eucl_of_list_aform)

lemma
  mem_Joints_appendE:
  assumes "x ∈ Joints (xs @ ys)"
  obtains x1 x2 where "x = x1 @ x2" "x1 ∈ Joints xs" "x2 ∈ Joints ys"
  using assms
  by (auto simp: Joints_def valuate_def)

lemma c1_info_of_appr_subsetI1:
  fixes X1::"'b::executable_euclidean_space set"
  assumes subset: "{eucl_of_list (map Inf_aform (fst R)) .. eucl_of_list (map Sup_aform (fst R))} ⊆ X1"
  assumes len: "length (fst R) = DIM('b)"
  shows "aform.c1_info_of_appr R ⊆ X1 × UNIV"
  using len
  apply (auto simp: aform.c1_info_of_appr_def flow1_of_list_def
      split: option.splits
      intro!: subsetD[OF subset] elim!: mem_Joints_appendE)
  subgoal
    by (auto intro!: eucl_of_list_mem_eucl_of_list_aform eucl_of_list_map_Inf_aform_leI eucl_of_list_map_Sup_aform_geI)
  subgoal
    by (auto intro!: eucl_of_list_mem_eucl_of_list_aform eucl_of_list_map_Inf_aform_leI eucl_of_list_map_Sup_aform_geI)
  subgoal
    apply (subst (2) eucl_of_list_take_DIM[symmetric, OF refl])
      apply (auto simp: min_def)
    apply (simp add: Joints_imp_length_eq eucl_of_list_map_Inf_aform_leI eucl_of_list_mem_eucl_of_list_aform)
    apply (simp add: Joints_imp_length_eq eucl_of_list_map_Inf_aform_leI eucl_of_list_mem_eucl_of_list_aform)
      done
  subgoal
    apply (subst (2) eucl_of_list_take_DIM[symmetric, OF refl])
      apply (auto simp: min_def)
    by (simp add: Joints_imp_length_eq eucl_of_list_map_Sup_aform_geI eucl_of_list_mem_eucl_of_list_aform)
  done

lemmas [simp] = compute_tdev

syntax product_aforms::"(real aform) list ⇒ (real aform) list ⇒ (real aform) list"
  (infixr "×⇩a" 70)

lemma matrix_inner_Basis_list:
  includes vec_syntax
  assumes "k < CARD('n) * CARD('m)"
  shows "(f::(('n::enum rvec, 'm::enum) vec)) ∙ Basis_list ! k =
    vec_nth (vec_nth f (enum_class.enum ! (k div CARD('n)))) (enum_class.enum ! (k mod CARD('n)))"
proof -
  have "f ∙ Basis_list ! k =
    (∑x∈UNIV.
       ∑xa∈UNIV.
         if enum_class.enum ! (k mod CARD('n)) = xa ∧ enum_class.enum ! (k div CARD('n)) = x then f $ x $ xa else 0)"
    using assms
    unfolding inner_vec_def
    apply (auto simp: Basis_list_vec_def concat_map_map_index)
    apply (subst (2) sum.cong[OF refl])
     apply (subst sum.cong[OF refl])
      apply (subst (2) vec_nth_Basis2)
       apply (force simp add: Basis_vec_def Basis_list_real_def)
      apply (rule refl)
     apply (rule refl)
    apply (auto simp: if_distribR if_distrib axis_eq_axis Basis_list_real_def cong: if_cong)
    done
  also have "… = f $ enum_class.enum ! (k div CARD('n)) $ enum_class.enum ! (k mod CARD('n))"
    apply (subst if_conn)
    apply (subst sum.delta')
     apply simp
    by (simp add: sum.delta')
  finally show ?thesis by simp
qed

lemma list_of_eucl_matrix:
  includes vec_syntax
  shows "(list_of_eucl (M::(('n::enum rvec, 'm::enum) vec))) =
    concat (map (λi. map (λj. M $ (enum_class.enum ! i)$ (enum_class.enum ! j) )
      [0..<CARD('n)]) [0..<CARD('m)])"
  by (auto intro!: nth_equalityI simp: length_concat o_def sum_list_distinct_conv_sum_set ac_simps
      concat_map_map_index matrix_inner_Basis_list)

lemma axis_eq_eucl_of_list:
  "(axis i 1::'n::enum rvec) = eucl_of_list ((replicate CARD('n) 0)[index enum_class.enum i := 1])"
  apply (auto intro!: euclidean_eqI[where 'a="'n rvec"]
      simp: eucl_of_list_inner nth_list_update )
   apply (auto simp: index_Basis_list_axis1[symmetric])
  by (simp add: inner_axis inner_commute vec_nth_Basis)

lemma eucl_of_list_012: "eucl_of_list [vec_nth A 0, vec_nth A 1, vec_nth A 2] = A" for A::"3 rvec"
  apply vector
  apply (auto simp: vec_nth_eucl_of_list_eq index_Basis_list_axis1)
  using exhaust_3 three_eq_zero by blast


definition "ldec x = (case quotient_of x of (i, j) ⇒ real_of_float (lapprox_rat 20 i j))"
definition "udec x = (case quotient_of x of (i, j) ⇒ real_of_float (rapprox_rat 20 i j))"

lemma ldec: "ldec x ≤ real_of_rat x"
  and udec: "real_of_rat x ≤ udec x"
   apply (auto simp: ldec_def udec_def split: prod.splits
      intro!: lapprox_rat[le] rapprox_rat[ge])
   apply (metis Fract_of_int_quotient less_eq_real_def less_int_code(1) of_rat_rat
      quotient_of_denom_pos quotient_of_div)
  apply (metis Fract_of_int_quotient less_eq_real_def less_int_code(1) of_rat_rat
      quotient_of_denom_pos quotient_of_div)
  done

context includes floatarith_notation begin

definition "matrix_of_degrees⇩e =
  (let
    ur = Rad_of (Var 0);
    vr = Rad_of (Var 1)
   in
    [Cos ur, Cos vr, 0,
     Sin ur, Sin vr, 0,
     0, 0, 0])"

definition "matrix_of_degrees u v =
  approx_floatariths 30 matrix_of_degrees⇩e (aforms_of_point ([u, v, 0]))"

lemma interpret_floatariths_matrix_of_degrees:
  "interpret_floatariths matrix_of_degrees⇩e
    (([u::real, v::real, 0]))
   =
  [cos (rad_of u), cos (rad_of v), 0,
   sin (rad_of u), sin (rad_of v), 0,
   0, 0, 0]"
  by (auto simp: matrix_of_degrees⇩e_def Let_def inverse_eq_divide)

definition "num_options p sstep m N a projs print_fun =
  ⦇
    precision = p,
    adaptive_atol = FloatR 1 (- a),
    adaptive_rtol = FloatR 1 (- a),
    method_id = 2,
    start_stepsize  = FloatR 1 (- sstep),
    iterations = 40,
    halve_stepsizes = 40,
    widening_mod = 10,
    rk2_param = FloatR 1 0,
    default_reduce = correct_girard (p) (m) (N),
    printing_fun = (λa b.
        let
           _ = fold (λ(x, y, c) _.
               print_fun (String.implode (shows_segments_of_aform (x) (y) b c ''⏎''))) projs ();
           _ = print_fun (String.implode (''# '' @ shows_box_of_aforms_hr (b) '''' @ ''⏎''))
        in
         ()
    ),
    tracing_fun = (λa b.
      let
        _ = print_fun (String.implode (''# '' @ a @ ''⏎''))
      in case b of Some b ⇒
        (let
          _ = ()
        in print_fun (String.implode (''# '' @ shows_box_of_aforms_hr (b) '''' @ ''⏎'')))
        | None ⇒ ())
  ⦈"

definition "num_options_c1 p sstep m N a projs dcolors print_fun =
  (let
    no = num_options p sstep m N a (map (λ(x, y, c, ds). (x, y, c)) projs) print_fun;
    D = length dcolors
  in no
    ⦇printing_fun:=
      (λa b. let _ = printing_fun no a b
          in if a then ()
          else fold (λ(x, y, c, ds) _. print_fun
              (String.implode (shows_aforms_vareq D [(x, y)] ds
                dcolors
                dcolors
            (FloatR 1 (-1)) ''# no C1 info'' b ''''))) projs ()
      )
    ⦈)"

definition "num_options_code p sstep m N a projs print_fun =
  num_options (nat_of_integer p) (int_of_integer sstep) (nat_of_integer m) (nat_of_integer N)
    (int_of_integer a) (map (λ(i, j, k). (nat_of_integer i, nat_of_integer j, k)) projs) print_fun"

definition "ro s n M g0 g1 inter_step =
  ⦇max_tdev_thres = FloatR 1 s,
      pre_split_reduce = correct_girard 30 n M,
      pre_inter_granularity = FloatR 1 g0,
      post_inter_granularity = (FloatR 1 g1),
      pre_collect_granularity = FloatR 1 g0,
      max_intersection_step = FloatR 1 inter_step⦈"

definition "code_ro s n m g0 g1 inter_step =
  ro (int_of_integer s) (nat_of_integer n) (nat_of_integer m) (int_of_integer g0) (int_of_integer g1) (int_of_integer inter_step)"

fun xsec:: "real ⇒ real × real ⇒ real × real ⇒ (real list × real list) × real list sctn"
  where "xsec x (y0, y1) (z0, z1) = (([x, y0, z0], [x, y1, z1]), Sctn [1,0,0] x)"
fun xsec':: "real ⇒ real × real ⇒ real × real ⇒ (real list × real list) × real list sctn"
  where "xsec' x (y0, y1) (z0, z1) = (([x, y0, z0], [x, y1, z1]), Sctn [-1,0,0] (-x))"

fun ysec:: "real × real ⇒ real ⇒ real × real ⇒ (real list × real list) × real list sctn"
  where "ysec (x0, x1) y (z0, z1) = (([x0, y, z0], [x1, y, z1]), Sctn [0, 1,0] y)"
fun ysec':: "real × real ⇒ real ⇒ real × real ⇒ (real list × real list) × real list sctn"
  where "ysec' (x0, x1) y (z0, z1) = (([x0, y, z0], [x1, y, z1]), Sctn [0, -1,0] (-y))"

fun zsec:: "real × real ⇒ real × real ⇒ real ⇒ (real list × real list) × real list sctn"
  where "zsec (x0, x1) (y0, y1) z = (([x0, y0, z], [x1, y1, z]), Sctn [0, 0, 1] z)"
fun zsec':: "real × real ⇒ real × real ⇒ real ⇒ (real list × real list) × real list sctn"
  where "zsec' (x0, x1) (y0, y1) z = (([x0, y0, z], [x1, y1, z]), Sctn [0, 0, -1] (-z))"


fun xsec2:: "real ⇒ real × real ⇒ (real list × real list) × real list sctn"
  where "xsec2 x (y0, y1) = (([x, y0], [x, y1]), Sctn [1,0] x)"
fun xsec2':: "real ⇒ real × real ⇒(real list × real list) × real list sctn"
  where "xsec2' x (y0, y1) = (([x, y0], [x, y1]), Sctn [-1,0] (-x))"

fun ysec2:: "real × real ⇒ real ⇒ (real list × real list) × real list sctn"
  where "ysec2 (x0, x1) y = (([x0, y], [x1, y]), Sctn [0, 1] y)"
fun ysec2':: "real × real ⇒ real ⇒ (real list × real list) × real list sctn"
  where "ysec2' (x0, x1) y = (([x0, y], [x1, y]), Sctn [0, -1] (-y))"

fun ysec4:: "real × real ⇒ real ⇒ real × real ⇒ real × real ⇒ (real list × real list) × real list sctn"
  where "ysec4 (x0, x1) y (z0, z1) (m0, m1) = (([x0, y, z0, m0], [x1, y, z1, m1]), Sctn [0, 1,0, 0] (y))"

fun ysec4':: "real × real ⇒ real ⇒ real × real ⇒ real × real ⇒ (real list × real list) × real list sctn"
  where "ysec4' (x0, x1) y (z0, z1) (m0, m1) = (([x0, y, z0, m0], [x1, y, z1, m1]), Sctn [0, -1,0, 0] (-y))"

definition "code_sctn N n c = Sctn ((replicate (nat_of_integer N) (0::real))[nat_of_integer n := 1]) c"
definition "code_sctn' N n c = Sctn ((replicate (nat_of_integer N) (0::real))[nat_of_integer n := -1]) (-c)"

definition "lrat i j = real_of_float (lapprox_rat 30 (int_of_integer i) (int_of_integer j))"
definition "urat i j = real_of_float (lapprox_rat 30 (int_of_integer i) (int_of_integer j))"

definition [simp]: "TAG_optns (optns::string × ((String.literal ⇒ unit) ⇒ (real aform) numeric_options)) = True"
lemma TAG_optns: "P ⟹ (TAG_optns optns ⟹ P)"
  by auto

definition [simp]: "TAG_reach_optns (roi::real aform reach_options) = True"
lemma TAG_reach_optns: "P ⟹ (TAG_reach_optns optns ⟹ P)"
  by auto

definition [simp]: "TAG_sctn (b::bool) = True"
lemma TAG_sctn: "P ⟹ (TAG_sctn optns ⟹ P)"
  by auto

subsection ‹Integrals and Computation›

lemma has_vderiv_on_PairD:
  assumes "((λt. (f t, g t)) has_vderiv_on fg') T"
  shows "(f has_vderiv_on (λt. fst (fg' t))) T" "(g has_vderiv_on (λt. snd (fg' t))) T"
proof -
  from assms have "((λx. (f x, g x)) has_derivative (λxa. xa *⇩R fg' t)) (at t within T)"
     if "t ∈ T" for t
    by (auto simp: has_vderiv_on_def has_vector_derivative_def that
        intro!: derivative_eq_intros)
  from diff_chain_within[OF this has_derivative_fst[OF has_derivative_ident]]
    diff_chain_within[OF this has_derivative_snd[OF has_derivative_ident]]
  show "(f has_vderiv_on (λt. fst (fg' t))) T" "(g has_vderiv_on (λt. snd (fg' t))) T"
    by (auto simp: has_vderiv_on_def has_vector_derivative_def o_def)
qed

lemma solves_autonomous_odeI:
  assumes "((λt. (t, phi t)) solves_ode (λt x. (1, f (fst x) (snd x)))) S (T × X)"
  shows "(phi solves_ode f) S X"
proof (rule solves_odeI)
  from solves_odeD[OF assms]
  have
    "((λt. (t, phi t)) has_vderiv_on (λt. (1, f (fst (t, phi t)) (snd (t, phi t))))) S"
    "⋀t. t ∈ S ⟹ (phi t) ∈ X"
    by auto
  from has_vderiv_on_PairD(2)[OF this(1)] this(2)
  show "(phi has_vderiv_on (λt. f t (phi t))) S" "⋀t. t ∈ S ⟹ phi t ∈ X"
    by auto
qed

lemma integral_solves_autonomous_odeI:
  fixes f::"real ⇒ 'a::executable_euclidean_space"
  assumes "(phi solves_ode (λt _. f t)) {a .. b} X" "phi a = 0"
  assumes "a ≤ b"
  shows "(f has_integral phi b) {a .. b}"
proof -
  have "(f has_integral phi b - phi a) {a..b}"
    apply (rule fundamental_theorem_of_calculus[of a b phi f])
    unfolding has_vderiv_on_def[symmetric]
     apply fact
    using solves_odeD[OF assms(1)]
    by (simp add: has_vderiv_on_def)
  then show ?thesis by (simp add: assms)
qed

lemma zero_eq_eucl_of_list_rep_DIM: "(0::'a::executable_euclidean_space) = eucl_of_list (replicate (DIM('a)) 0)"
  by (auto intro!: euclidean_eqI[where 'a='a] simp: eucl_of_list_inner)

lemma zero_eq_eucl_of_list_rep: "(0::'a::executable_euclidean_space) = eucl_of_list (replicate D 0)"
  if "D ≥ DIM('a)"
proof -
  from that have "replicate D (0::real) = replicate (DIM('a)) 0 @ replicate (D - DIM('a)) 0"
    by (auto simp: replicate_add[symmetric])
  also have "eucl_of_list … = (eucl_of_list (replicate DIM('a) 0)::'a)"
    by (rule eucl_of_list_append_zeroes)
  also have "… = 0"
    by (rule zero_eq_eucl_of_list_rep_DIM[symmetric])
  finally show ?thesis by simp
qed

lemma one_has_ivl_integral: "((λx. 1::real) has_ivl_integral b - a) a b"
  using has_integral_const_real[of "1::real" a b] has_integral_const_real[of "1::real" b a]
  by (auto simp: has_ivl_integral_def split: if_splits)

lemma Joints_aforms_of_point_self[simp]: "xs ∈ Joints (aforms_of_point xs)"
  by (simp add: aforms_of_point_def)

lemma bind_eq_dRETURN_conv:
  "(f ⤜ g = dRETURN S) ⟷ (∃R. f = dRETURN R ∧ g R = dRETURN S)"
  by (cases f) auto

end

lemma list_of_eucl_memI: "list_of_eucl (x::'x::executable_euclidean_space) ∈ S"
  if "x ∈ eucl_of_list ` S" "⋀x. x ∈ S ⟹ length x = DIM('x)"
  using that
  by auto

lemma Joints_aforms_of_ivls_append_point:
  "Joints (xs @ aforms_of_ivls p p) = (λx. x @ p) ` Joints xs"
  using aform.set_of_appr_of_ivl_append_point[unfolded aform_ops_def approximate_set_ops.simps] .


context ode_interpretation begin

theorem solves_one_step_ivl:
  assumes T: "T ⊆ {t1 .. t2}"
  assumes X: "X ⊆ {eucl_of_list lx .. eucl_of_list ux}" "length lx = DIM('a)" "length ux = DIM('a)"
  assumes S: "{eucl_of_list ls::'a .. eucl_of_list us} ⊆ S"
  assumes lens: "length ls = DIM('a)" "length us = DIM('a)" ― ‹TODO: this could be verified›
  assumes [simp]: "aform.D odo = DIM('a)"
  assumes r: "solves_one_step_until_time_aform optns odo ((1,1), aforms_of_ivls lx ux, None) t1 t2 (ls, us) None"
  shows "t ∈ T ⟶ x0 ∈ X ⟶ t ∈ existence_ivl0 x0 ∧ flow0 x0 t ∈ S"
proof (intro impI)
  assume t: "t ∈ T" and x0: "x0 ∈ X"
  from S have S: "set_of_lvivl (ls, us) ⊆ S"
    by (auto simp: set_of_lvivl_def set_of_ivl_def)
  have lens: "length (fst (ls, us)) = CARD('n)" "length (snd (ls, us)) = CARD('n)"
    by (auto simp: lens)
  have x0: "list_of_eucl x0 ∈ Joints (aforms_of_ivls lx ux)"
    apply (rule aforms_of_ivls)
    subgoal by (simp add: X)
    subgoal by (simp add: X)
    using subsetD[OF X(1) x0]
    apply (auto simp: eucl_le[where 'a='a] X)
    apply (metis assms(3) dim length_Basis_list list_of_eucl_eucl_of_list list_of_eucl_nth nth_Basis_list_in_Basis)
    apply (metis assms(4) dim length_Basis_list list_of_eucl_eucl_of_list list_of_eucl_nth nth_Basis_list_in_Basis)
    done
  from t T have t: "t ∈ {t1 .. t2}" by auto
  show "t ∈ existence_ivl0 x0 ∧ flow0 x0 t ∈ S"
    by (rule one_step_result12[OF r aform.c1_info_of_appre_c0_I[OF x0] t S subset_UNIV lens])
      (auto simp: aform.c1_info_invare_None lens X)
qed

theorem solves_one_step_ivl':
  assumes T: "T ⊆ {t1 .. t2}"
  assumes X: "X ⊆ {eucl_of_list lx .. eucl_of_list ux}"
    "length lx = DIM('a)" "length ux = DIM('a)"
  assumes DS: "list_interval lds uds ⊆ list_interval ld ud"
    and lends: "length lds = DIM('a)*DIM('a)" "length uds = DIM('a)*DIM('a)"
  assumes S: "{eucl_of_list ls::'a .. eucl_of_list us} ⊆ S"
  assumes lens0: "length ls = DIM('a)" "length us = DIM('a)" ― ‹TODO: this could be verified›
    "length dx0s = DIM('a)*DIM('a)"
  assumes [simp]: "aform.D odo = DIM('a)"
  assumes r: "solves_one_step_until_time_aform optns odo
    ((1,1), aforms_of_ivls lx ux, Some (aforms_of_point dx0s)) t1 t2 (ls, us) (Some (lds, uds))"
  shows "t ∈ T ⟶ x0 ∈ X ⟶ t ∈ existence_ivl0 x0 ∧ flow0 x0 t ∈ S ∧
    Dflow x0 t o⇩L blinfun_of_list dx0s ∈ blinfuns_of_lvivl (ld, ud)"
proof (intro impI)
  assume t: "t ∈ T" and x0: "x0 ∈ X"
  from S have S: "set_of_lvivl (ls, us) ⊆ S"
    by (auto simp: set_of_lvivl_def set_of_ivl_def)
  have lens: "length (fst (ls, us)) = CARD('n)" "length (snd (ls, us)) = CARD('n)"
    by (auto simp: lens0)
  have x0: "list_of_eucl x0 ∈ Joints (aforms_of_ivls lx ux)"
    apply (rule aforms_of_ivls)
    subgoal by (simp add: X)
    subgoal by (simp add: X)
    using subsetD[OF X(1) x0]
    apply (auto simp: eucl_le[where 'a='a] X)
    apply (metis assms(3) dim length_Basis_list list_of_eucl_eucl_of_list list_of_eucl_nth nth_Basis_list_in_Basis)
    apply (metis assms(4) dim length_Basis_list list_of_eucl_eucl_of_list list_of_eucl_nth nth_Basis_list_in_Basis)
    done
  have x0dx0: "(x0, blinfun_of_list dx0s) ∈
      aform.c1_info_of_appre ((1, 1), aforms_of_ivls lx ux, Some (aforms_of_point dx0s))"
    apply (auto simp: aform.c1_info_of_appre_def aform.c1_info_of_appr_def)
    apply (rule image_eqI[where x="list_of_eucl x0