section {* Euler Method on Affine Forms: Code *}
theory Euler_Affine_Code
imports
  Print
  "../../Collections/Iterator/SetIteratorOperations"
  "../../Collections/ICF/impl/ListSetImpl"
  "../../Collections/ICF/ICF_Refine_Monadic"
  "~~/src/HOL/Library/AList"
  "~~/src/HOL/Library/Monad_Syntax"
  "~~/src/HOL/Library/While_Combinator"
  "../Numerics/Runge_Kutta"
  "../../Affine_Arithmetic/Affine_Arithmetic"
begin
text {*\label{sec:euleraformcode}*}
lemma (in conditionally_complete_lattice) cInf_insert2:
  "X \<noteq> {} \<Longrightarrow> bdd_below X \<Longrightarrow> Inf (insert a (insert b X)) = inf (inf a b) (Inf X)"
  using inf_le1[of "inf a b" "Inf X"] inf_le1[of a b]
  by (intro cInf_eq_non_empty) (auto intro: le_infI2 cInf_lower cInf_greatest)

lemma (in conditionally_complete_lattice) cSup_insert2:
  "X \<noteq> {} \<Longrightarrow> bdd_above X \<Longrightarrow> Sup (insert a (insert b X)) = sup (sup a b) (Sup X)"
  using order_trans[OF sup_ge1[of a b] sup_ge1[of "sup a b" "Sup X"]]
    order_trans[OF sup_ge2[of b a] sup_ge1[of "sup a b" "Sup X"]]
  by (intro cSup_eq_non_empty) (auto intro: le_supI2 cSup_upper cSup_least)

lemma (in conditionally_complete_lattice) Inf_set_fold_inf:
  shows "Inf (set (x#xs)) = fold inf xs x"
proof (induct xs arbitrary: x)
  case Nil thus ?case by simp
next
  case (Cons x xs y)
  from Cons[of "inf x y"]
  show ?case
    apply (cases "xs = []")
    apply (simp add: cInf_insert inf_commute)
    apply (simp)
    apply (auto simp: cInf_insert2)
    apply (auto simp: cInf_insert inf_commute)
    done
qed

lemma (in conditionally_complete_lattice) Sup_set_fold_sup:
  shows "Sup (set (x#xs)) = fold sup xs x"
proof (induct xs arbitrary: x)
  case Nil thus ?case by simp
next
  case (Cons x xs y)
  from Cons[of "sup x y"]
  show ?case
    apply (cases "xs = []")
    apply (simp add: cSup_insert sup_commute)
    apply (simp)
    apply (auto simp: cSup_insert2)
    apply (auto simp: cSup_insert sup_commute)
    done
qed


definition "THE_NRES x = (case x of Some r \<Rightarrow> RETURN r | None \<Rightarrow> SUCCEED)"

lemma THE_NRES_rule[refine_vcg]: "(\<And>r. x = Some r \<Longrightarrow> f r) \<Longrightarrow> THE_NRES x \<le> SPEC f"
  by (auto simp: THE_NRES_def split: option.split)

schematic_lemma transfer_THE_NRES[refine_transfer]: "nres_of ?X \<le> THE_NRES x"
  unfolding THE_NRES_def
  by (refine_transfer)


text{*\label{sec:euleraformcode}*}

(* TODO: not interesting anymore *)
datatype 'a sctn = Sctn (normal: 'a) (pstn: real)

definition "le_halfspace sctn x \<longleftrightarrow> x \<bullet> normal sctn \<le> pstn sctn"
definition "plane_of sctn = {x. x \<bullet> normal sctn = pstn sctn}"

record ('a, 'b, 'c) options =
  precision :: nat
  tolerance :: real
  stepsize :: real
  iterations :: nat
  halve_stepsizes :: nat
  widening_mod :: nat
  rk2_param :: real
  max_tdev_thres :: real
  presplit_summary_tolerance :: real
  collect_granularity :: real
  start_section :: "'a sctn"
  checkpoint_sections :: "'a sctn list"
  stop_sections :: "'a sctn list"
  safe_appr ::"'b \<Rightarrow> bool"
  snap_sections :: "real"
  min_section_distance :: "real"
  next_section_distance_factor :: "real"
  next_section_turn_distance_factor :: "real"
  printing_fun :: "'b \<Rightarrow> unit"
  irect_mod :: "nat"
  result_fun :: "nat * real * 'b * (real * 'b * real * 'b) list \<Rightarrow> 'c"
  max_intersection_step :: real
  adaptive_atol :: real --"TODO: implement relative error tolerance?!"
  method_id :: nat

text {* TODO: options as parameter to the locale!!! *}

locale approximate_sets0 =
  fixes appr_of_ivl::"'a::executable_euclidean_space \<Rightarrow> 'a \<Rightarrow> 'b"
  fixes msum_appr::"'b \<Rightarrow> 'b \<Rightarrow> 'b"
  fixes set_of_appr::"'b \<Rightarrow> 'a set"
  fixes set_of_apprs::"'b list \<Rightarrow> 'a list set"
  fixes inf_of_appr::"'b \<Rightarrow> 'a"
  fixes sup_of_appr::"'b \<Rightarrow> 'a"
  fixes add_appr::"('a, 'b, 'c) options \<Rightarrow> 'b \<Rightarrow> 'b \<Rightarrow> 'b list \<Rightarrow> 'b option"
  fixes scale_appr::"('a, 'b, 'c) options \<Rightarrow> real \<Rightarrow> real \<Rightarrow> 'b \<Rightarrow> 'b list \<Rightarrow> 'b option"
  fixes scale_appr_ivl::"('a, 'b, 'c) options \<Rightarrow> real \<Rightarrow> real \<Rightarrow> 'b \<Rightarrow> 'b list \<Rightarrow> 'b option"
  fixes split_appr::"('a, 'b, 'c) options \<Rightarrow> 'b \<Rightarrow> 'b list"
  fixes support_appr::"('a, 'b, 'c) options \<Rightarrow> 'a \<Rightarrow> 'b \<Rightarrow> real"
  fixes appr_inf_inner::"('a, 'b, 'c) options \<Rightarrow> 'b \<Rightarrow> 'a \<Rightarrow> real"
  fixes appr_sup_inner::"('a, 'b, 'c) options \<Rightarrow> 'b \<Rightarrow> 'a \<Rightarrow> real"
  fixes inter_appr_plane::"('a, 'b, 'c) options \<Rightarrow> 'b \<Rightarrow> 'a sctn \<Rightarrow> 'b option"
  fixes reduce_appr::"('a, 'b, 'c) options \<Rightarrow> 'b \<Rightarrow> 'b"
  fixes independent_appr::"'b \<Rightarrow> 'b \<Rightarrow> 'b"
  fixes width_appr::"('a, 'b, 'c) options \<Rightarrow> 'b \<Rightarrow> real"
begin

definition ivl_appr_of_appr::"'b \<Rightarrow> 'b" where
  "ivl_appr_of_appr x = (appr_of_ivl (inf_of_appr x) (sup_of_appr x))"

definition
  "ivl_of_apprs xs = appr_of_ivl (Inf (set (map inf_of_appr xs))) (Sup (set (map sup_of_appr xs)))"

lemma compute_ivl_of_apprs:
  "ivl_of_apprs [] = appr_of_ivl (Inf {}) (Sup {})"
  "ivl_of_apprs (x#xs) =
    appr_of_ivl (fold inf (map inf_of_appr xs) (inf_of_appr x)) (fold sup (map sup_of_appr xs) (sup_of_appr x))"
  by (simp_all add: ivl_of_apprs_def Inf_set_fold_inf[symmetric] Sup_set_fold_sup[symmetric])

definition "appr_above_sctn optns X sctn = (appr_inf_inner optns X (normal sctn) > pstn sctn)"

definition "appr_below_sctn optns X sctn = (appr_sup_inner optns X (normal sctn) \<le> pstn sctn)"

definition "intersects_appr optns X sctn = (appr_inf_inner optns X (normal sctn) \<le> pstn sctn \<and> appr_sup_inner optns X (normal sctn) \<ge> pstn sctn)"

definition disjoint_apprs where
  "disjoint_apprs X Y =
    (let iX = inf_of_appr X; sX = sup_of_appr X; iY = inf_of_appr Y; sY = sup_of_appr Y
    in list_ex (\<lambda>i. sX \<bullet> i < iY \<bullet> i \<or> sY \<bullet> i < iX \<bullet> i) Basis_list)"


text {* irects are for easy and exact exchange and management of sets as small hypercubes *}

definition appr_of_irect::"('a, 'b, 'c) options \<Rightarrow> int list \<Rightarrow> 'b"
  where "appr_of_irect optns is =
    appr_of_ivl
      (\<Sum>(i, b)\<leftarrow>zip is Basis_list. FloatR i (-irect_mod optns) *\<^sub>R b)
      (\<Sum>(i, b)\<leftarrow>zip is Basis_list. FloatR (i + 1) (-irect_mod optns) *\<^sub>R b)"

definition "project_appr X b y =
  (let i = inf_of_appr X; s = sup_of_appr X in
   appr_of_ivl (i + (y - i \<bullet> b) *\<^sub>R b) (s + (y - s \<bullet> b) *\<^sub>R b))"

definition "irects_of_appr optns X =
  (let
    I = inf_of_appr X;
    S = sup_of_appr X;
    i0s = [floor ((I\<bullet>b) * 2^(irect_mod optns)). b\<leftarrow>Basis_list];
    s0s = [floor ((S\<bullet>b) * 2^(irect_mod optns)). b\<leftarrow>Basis_list]
  in product_lists (map (\<lambda>(i, s). [i..s]) (zip i0s s0s)))"

definition "irects_of_apprs optns XS = remdups (concat (map (irects_of_appr optns) XS))"

definition "print_irects irects =
  (let _ = map (\<lambda>is.
    let _ = map (\<lambda>j. let _ = print (int_to_string j) in print (STR '' '')) is in println (STR '''')) irects
  in ())"

end

declare approximate_sets0.appr_of_irect_def[code]
declare approximate_sets0.project_appr_def[code]
declare approximate_sets0.irects_of_appr_def[code]
declare approximate_sets0.irects_of_apprs_def[code]
declare approximate_sets0.print_irects_def[code]
declare approximate_sets0.ivl_appr_of_appr_def[code]
declare approximate_sets0.appr_above_sctn_def[code]
declare approximate_sets0.appr_below_sctn_def[code]
declare approximate_sets0.intersects_appr_def[code]
declare approximate_sets0.disjoint_apprs_def[code]
declare approximate_sets0.compute_ivl_of_apprs[code]
declare approximate_sets0.project_appr_def[code]

type_synonym 'a enclosure = "nat * real * 'a * (real * 'a * real * 'a) list"

locale approximate_ivp0 = approximate_sets0
  appr_of_ivl msum_appr set_of_appr set_of_apprs inf_of_appr sup_of_appr add_appr
  scale_appr scale_appr_ivl
  split_appr support_appr appr_inf_inner appr_sup_inner
  inter_appr_plane reduce_appr independent_appr width_appr
for appr_of_ivl msum_appr set_of_appr set_of_apprs inf_of_appr
  and sup_of_appr::"'b \<Rightarrow> 'a::executable_euclidean_space"
  and add_appr:: "('a, 'b, 'c) options \<Rightarrow> 'b \<Rightarrow> 'b \<Rightarrow> 'b list \<Rightarrow> 'b option"
  and scale_appr scale_appr_ivl split_appr support_appr appr_inf_inner appr_sup_inner
    inter_appr_plane reduce_appr independent_appr width_appr +
  fixes ode_approx::"('a, 'b, 'c) options \<Rightarrow> 'b list \<Rightarrow> 'b option"
  fixes ode_d_approx:: "('a, 'b, 'c) options \<Rightarrow> 'b list \<Rightarrow> 'b option"
  fixes ode_d2_approx:: "('a, 'b, 'c) options \<Rightarrow> 'b list \<Rightarrow> 'b option"
begin

abbreviation "extend_appr \<equiv> \<lambda>x l u. msum_appr x (appr_of_ivl l u)"

definition P_appr::"('a, 'b, 'c) options \<Rightarrow> 'b \<Rightarrow> real \<Rightarrow> 'b \<Rightarrow> 'b option" where
  "P_appr optns X0 h X = map_option (\<lambda>Y.
    extend_appr X0 (inf 0 (h *\<^sub>R inf_of_appr Y))
                   (sup  0 (h *\<^sub>R sup_of_appr Y)))
    (ode_approx optns [X])"

fun P_iter::"('a, 'b, 'c) options \<Rightarrow> 'b \<Rightarrow> real \<Rightarrow> nat \<Rightarrow> 'b \<Rightarrow> 'b option" where
  "P_iter optns X0 h 0 X =
    (let _ = print (STR ''=P_iter failed: '');
      _ = print_eucl (inf_of_appr X);
      _ = print (STR '' - '');
      _ = print_eucl (sup_of_appr X);
      _ = println (STR '''') in None)"
| "P_iter optns X0 h (Suc i) X =
    bind_err (STR ''=P_appr failed'') (P_appr optns X0 h X) (\<lambda>X'.
      let (l', u') = (inf_of_appr X', sup_of_appr X') in
      let (l, u) = (inf_of_appr X, sup_of_appr X) in
      if l \<le> l' \<and> u' \<le> u then Some X
      else P_iter optns X0 h i (appr_of_ivl (inf l' l - (if i mod (widening_mod optns) = 0 then abs (l' - l) else 0)) (sup u' u + (if i mod widening_mod optns = 0 then abs (u' - u) else 0))))"

fun cert_stepsize::"('a, 'b, 'c) options \<Rightarrow> 'b \<Rightarrow> real \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> (real * 'b) option" where
  "cert_stepsize optns X0 h n 0 = (let _ = println (STR ''=cert_stepsize failed'') in None)"
| "cert_stepsize optns X0 h n (Suc i) =
    (case P_iter optns (ivl_appr_of_appr X0) h n (ivl_appr_of_appr X0) of
      Some X' \<Rightarrow> Some (h, X')
    | None \<Rightarrow>
      let
        _ = print (STR ''=cert_stepsize failed on: '');
        _ = print_eucl (inf_of_appr X0);
        _ = print (STR '' - '');
        _ = print_eucl (sup_of_appr X0);
        _ = println (STR '''')
      in cert_stepsize optns X0 (h / 2) n i)"

lemma cert_stepsize_pos: "cert_stepsize optns X0 h n i = Some (h', cx) \<Longrightarrow> h > 0 \<Longrightarrow> h' > 0"
  by (induct i arbitrary: h h') (auto split: option.split_asm)

definition "euler_step optns X0 =
  bind_err (STR ''=certify stepsize failed'') (cert_stepsize optns X0 (stepsize optns) (iterations optns) (halve_stepsizes optns))
  (\<lambda>(h, CX).
      bind_err (STR ''=ode_approx X0 failed'') (ode_approx optns [X0])
      (\<lambda>X0'.
        bind_err (STR ''=ode_approx CX failed'') (ode_approx optns [CX])
        (\<lambda>F.
          bind_err (STR ''=ode_d_approx failed'') (ode_d_approx optns [CX, F])
          (\<lambda>D.
            bind_err (STR ''=scale_appr err failed'') (scale_appr optns (h*h) 2 (ivl_appr_of_appr D) [])
            (\<lambda>ERR.
                bind_err (STR ''=scale_appr euler failed'') (scale_appr optns h 1 X0' [X0])
                (\<lambda>S.
                  bind_err (STR ''=scale_appr_ivl euler failed'') (scale_appr_ivl optns 0 h X0' [X0])
                  (\<lambda>S'.
                    bind_err (STR ''=add_appr euler failed'') (add_appr optns X0 S [])
                    (\<lambda>X1.
                      bind_err (STR ''=add_appr euler_ivl failed'') (add_appr optns X0 S' [])
                      (\<lambda>CX1.
                        let
                          res = msum_appr X1 (ivl_appr_of_appr ERR);
                          res_ivl = msum_appr CX1 (appr_of_ivl (inf 0 (inf_of_appr ERR)) (sup 0 (sup_of_appr ERR)))
                        in
                        Some (h, width_appr optns ERR, res_ivl, res))))))))))"

definition "euler_step_nres optns X0 =
  THE_NRES (euler_step optns X0)"

schematic_lemma euler_step_nres_code_aux: "nres_of ?euler_step_dres \<le> euler_step_nres optns X"
  unfolding euler_step_nres_def
  by refine_transfer

concrete_definition euler_step_dres for optns X uses euler_step_nres_code_aux


definition "rk2_incr_nres optns h X =
  do {
    let p = rk2_param optns;
    FX \<leftarrow> THE_NRES (ode_approx optns [X]);
    A \<leftarrow> THE_NRES (scale_appr optns (p * 2 - 1) (p * 2) FX [X]); (* TODO: approximate calculation of coeff?! *)

    S \<leftarrow> THE_NRES (scale_appr optns (h * p) 1 FX [A, X]);
    XpS \<leftarrow> THE_NRES (add_appr optns X S [A, X]);
    F2 \<leftarrow> THE_NRES (ode_approx optns [XpS, X, A]);
    B \<leftarrow> THE_NRES (scale_appr optns 1 (p * 2) F2 [X, A]);

    S' \<leftarrow> THE_NRES (scale_appr_ivl optns 0 (h * p) FX [A, X]);
    XpS' \<leftarrow> THE_NRES (add_appr optns X S' [A, X]);
    F2' \<leftarrow> THE_NRES (ode_approx optns [XpS', X, A]);
    B' \<leftarrow> THE_NRES (scale_appr optns 1 (p * 2) F2' [X, A]);

    I1 \<leftarrow> THE_NRES (add_appr optns A B [X]);
    I1' \<leftarrow> THE_NRES (add_appr optns A B' [X]);
    RETURN (I1, I1')
  }"

schematic_lemma rk2_incr_dres: "nres_of ?rk2_incr_nres \<le> rk2_incr_nres optns h X"
  unfolding rk2_incr_nres_def
  by refine_transfer

concrete_definition rk2_incr_dres for optns h X uses rk2_incr_dres

definition "rk2_remainder_nres optns h xt xttt fxttt =
  do {
    let p = rk2_param optns;
    (* first term *)
    a \<leftarrow> THE_NRES (ode_d2_approx optns [xttt,fxttt,fxttt,xt]);
    b \<leftarrow> THE_NRES (ode_d_approx optns [xttt, fxttt, a, xt]);
    c \<leftarrow> THE_NRES (ode_d_approx optns [xttt, b, a, xt, fxttt]);
    d \<leftarrow> THE_NRES (add_appr optns a c [xt, xttt, fxttt]);
    e \<leftarrow> THE_NRES (scale_appr optns 1 6 d [xt, xttt, fxttt]);
    (* second term *)
    fxt \<leftarrow> THE_NRES (ode_approx optns [xt, e, xttt, fxttt]);
    f \<leftarrow> THE_NRES (scale_appr_ivl optns 0 (h * p) fxt [e, xt, xttt, fxttt]);
    fa \<leftarrow> THE_NRES (add_appr optns xt f [fxt, e, xttt, fxttt]);
    g \<leftarrow> THE_NRES (ode_d2_approx optns [fa,fxt,fxt,e, xt, xttt, fxttt]);
    h' \<leftarrow> THE_NRES (scale_appr optns (-p) 4 g [e, xt, xttt, fxttt]);
    i \<leftarrow> THE_NRES (add_appr optns h' e [xt, xttt, fxttt]);
    j \<leftarrow> THE_NRES (scale_appr optns (h^3) 1 i [xt, xttt, fxttt]);
    j' \<leftarrow> THE_NRES (scale_appr_ivl optns 0 (h^3) i [xt, xttt, fxttt]);
    RETURN (j, j')
  }"

schematic_lemma rk2_remainder_dres: "nres_of ?rk2_remainder_nres \<le> rk2_remainder_nres optns h xt xttt fxttt"
  unfolding rk2_remainder_nres_def
  by refine_transfer

concrete_definition rk2_remainder_dres for optns xt xttt fxttt uses rk2_remainder_dres
declare rk2_remainder_dres.refine[refine_transfer]

definition "rk2_step optns X0 =
  do {
    (h, CX0) \<leftarrow> THE_NRES (cert_stepsize optns X0 (stepsize optns) (iterations optns) (halve_stepsizes optns));
    let CX = independent_appr X0 CX0;
    X0' \<leftarrow> THE_NRES (ode_approx optns [X0]);
    F \<leftarrow> THE_NRES (ode_approx optns [CX, X0]);
    (R, R') \<leftarrow> rk2_remainder_nres optns h X0 CX F;
    (I, I') \<leftarrow> rk2_incr_nres optns h X0;
    S \<leftarrow> THE_NRES (scale_appr optns h 1 I [X0]);
    S' \<leftarrow> THE_NRES (scale_appr_ivl optns 0 h I' [X0]);
    X1 \<leftarrow> THE_NRES (add_appr optns X0 S []);
    CX1 \<leftarrow> THE_NRES (add_appr optns X0 S' []);
    let res = msum_appr X1 (ivl_appr_of_appr R);
    let res_ivl = msum_appr CX1 (ivl_appr_of_appr R');
    RETURN (h, width_appr optns R, res_ivl, res)
  }"

schematic_lemma rk2_step_nres_code_aux:
  notes rk2_remainder_dres.refine[refine_transfer]
    and rk2_incr_dres.refine[refine_transfer]
  shows "nres_of ?rk2_step_dres \<le> rk2_step optns X"
  unfolding rk2_step_def
  by refine_transfer

concrete_definition rk2_step_dres for optns X uses rk2_step_nres_code_aux
declare rk2_step_dres.refine[refine_transfer]

fun advance_euler::"('a, 'b, 'c) options \<Rightarrow> 'b enclosure option \<Rightarrow> 'b enclosure option"
where
  "advance_euler optns None = None"
| "advance_euler optns (Some (i, t, X, XS)) =
    (case euler_step optns X of
      Some (h, _, CX, X1) \<Rightarrow>
        let _ = printing_fun optns X1
        in Some (Suc i, t + h, X1, (t, CX, t + h, X1)#XS)
    | None \<Rightarrow> None)"

primrec euler_series where
  "euler_series optns t0 X0 0 = Some (0, t0, X0, [])"
| "euler_series optns t0 X0 (Suc i) = advance_euler optns (euler_series optns t0 X0 i)"

definition "split_appr_while optns =
  WHILE (\<lambda>XS. \<exists>X\<in>set XS. width_appr optns X > max_tdev_thres optns) (\<lambda>XS.
    RETURN (List.maps (\<lambda>x. if width_appr optns x > max_tdev_thres optns then split_appr optns x else [x]) XS)
  )"

schematic_lemma split_appr_while_dres: "nres_of ?split_appr_while_dres \<le> split_appr_while optns X"
  unfolding split_appr_while_def
  by refine_transfer

concrete_definition split_appr_while_dres for optns X uses split_appr_while_dres
declare split_appr_while_dres.refine[refine_transfer]

definition partition::"('a, 'b, 'c) options \<Rightarrow> 'b list \<Rightarrow> 'b list nres"
  where
 "partition optns xs = do {
    let r = ivl_of_apprs xs;
    (rs, ps) \<leftarrow>
      WHILE (\<lambda>(rs, ps). rs \<noteq> []) (\<lambda>(rs, ps). do {
        r \<leftarrow> SPEC (\<lambda>r. r \<in> set rs);
        rs' \<leftarrow> SPEC (\<lambda>rs'. r \<notin> set rs' \<and> set rs = insert r (set rs'));
        let is = filter (\<lambda>a. \<not> disjoint_apprs a r) xs;
        if is = [] then RETURN (rs', ps)
        else
          if width_appr optns r \<le> collect_granularity optns then
            RETURN (rs', r#ps)
          else do {
            let ivl = ivl_of_apprs is;
            let r' = appr_of_ivl (sup (inf_of_appr r) (inf_of_appr ivl)) (inf (sup_of_appr r) (sup_of_appr ivl));
            let r'' = split_appr optns r';
            RETURN (r''@rs', ps)
          }
      }) ([r], []);
    RETURN ps
  }"

definition partition_impl::"('a, 'b, 'c) options \<Rightarrow> 'b list \<Rightarrow> 'b list nres"
  where
 "partition_impl optns xs = do {
    let r = ivl_of_apprs xs;
    (rs, ps) \<leftarrow>
      WHILE (\<lambda>(rs, ps). rs \<noteq> []) (\<lambda>(rs, ps). do {
        let r = hd rs;
        let rs' = List.removeAll r rs;
        let is = filter (\<lambda>a. \<not> disjoint_apprs a r) xs;
        if is = [] then do {
          RETURN (rs', ps)
        } else
          if width_appr optns r \<le> collect_granularity optns then do {
            RETURN (rs', r#ps)
          } else do {
            let ivl = ivl_of_apprs is;
            let r' = appr_of_ivl (sup (inf_of_appr r) (inf_of_appr ivl)) (inf (sup_of_appr r) (sup_of_appr ivl));
            let r'' = split_appr optns r';
            RETURN (r''@rs', ps)
          }
      }) ([r], []);
    RETURN ps
  }"

lemma partition_impl: "partition_impl optns xs \<le> \<Down>Id (partition optns xs)"
  unfolding partition_impl_def partition_def
    Let_def[where s="println x" for x]
    Let_def[where s="print x" for x]
  apply (refine_rcg)
  apply refine_dref_type
  apply auto
  done

schematic_lemma partition_aux: "nres_of ?partition \<le> partition_impl optns xs"
  unfolding partition_impl_def
  by refine_transfer

concrete_definition partition_dres for optns xs uses partition_aux
declare partition_dres.refine[refine_transfer]

lemmas partition_partition_impl[refine_transfer] =
  order_trans[OF partition_dres.refine refine_IdD[OF partition_impl]]


definition euler_series_result::
  "('a, 'b, 'c) options \<Rightarrow> real \<Rightarrow> 'b \<Rightarrow> nat \<Rightarrow> 'c option"
  where [simp]: "euler_series_result optns t0 X0 i =
    map_option (result_fun optns) (euler_series optns t0 X0 i)"

lemma euler_series_print:
  "euler_series optns t0 X0 i =
    fold (\<lambda>a b.
      case b of
        None \<Rightarrow> None
      | Some (a', t0', X0', ress) \<Rightarrow>
        (case euler_step optns X0' of
          None \<Rightarrow> None
        | Some (h, _, CX, X1) \<Rightarrow>
          let
            _ = printing_fun optns X1
          in Some (Suc a', t0' + h, X1, (t0', CX, t0' + h, X1)#ress))) [0..<i] (Some (0, t0, X0, []))"
  unfolding Let_def
  by (induct i) (auto split: option.split)

definition "middle_appr X = (inf_of_appr X + sup_of_appr X) /\<^sub>R 2"

definition "transversal_directions f =
  (let
    I = inf_of_appr f;
    S = sup_of_appr f
  in listsum (map (\<lambda>b. (if I \<bullet> b \<le> 0 then if S \<bullet> b \<le> 0 then S \<bullet> b else 0 else if S \<bullet> b \<ge> 0 then I \<bullet> b else 0) *\<^sub>R b) (Basis_list::'a list)))"

definition "strongest_direction T =
  (let
    b = fold (\<lambda>b' b. if abs (T \<bullet> b') \<ge> abs (T \<bullet> b) then b' else b) (Basis_list::'a list) 0
  in
    (sgn (T \<bullet> b) *\<^sub>R b, abs (T \<bullet> b)))"

definition "strongest_transversal f = strongest_direction (transversal_directions f)"

definition "intersection_iteration_impl optns sctn X = do {
  F \<leftarrow> THE_NRES (ode_approx optns [X]);
  WHILE (\<lambda>(X, F, PS, CXS). \<not> appr_above_sctn optns X sctn \<and> fst (strongest_transversal F) = normal sctn)
    (\<lambda>(X, F, PS, CXS).
      do {
        (_, _, CX', X') \<leftarrow> (if method_id optns = 2 then rk2_step else euler_step_nres) optns X;
        P \<leftarrow> THE_NRES (inter_appr_plane optns (reduce_appr optns CX') sctn);
        F' \<leftarrow> THE_NRES (ode_approx optns [X']);
        RETURN (X', F', P#PS, CX'#CXS)
      })
    (X, F, [], [X])
  }"

definition "do_intersection optns sctn X =
  do {
    (Y, F, PS, CXS) \<leftarrow> intersection_iteration_impl optns sctn X;
    RETURN (PS, CXS, if appr_above_sctn optns Y sctn then [] else [Y])}"

schematic_lemma do_intersection_dres: "nres_of ?do_intersection \<le> do_intersection optns sctn X"
  unfolding intersection_iteration_impl_def do_intersection_def euler_step_nres_def
    cond_application_beta
  by refine_transfer

concrete_definition do_intersection_dres for optns sctn X uses do_intersection_dres
declare do_intersection_dres.refine[refine_transfer]

definition "resolve_sctn optns sctn XS =
  WHILE (\<lambda>(XS, PS, RS). XS \<noteq> []) (\<lambda>(XS, PS, RS). do {
    (h, X) \<leftarrow> SPEC (\<lambda>hX. hX \<in> set XS);
    XS' \<leftarrow> SPEC (\<lambda>XS'. (h, X) \<notin> set XS' \<and> set XS = insert (h, X) (set XS'));
    if width_appr optns X > max_tdev_thres optns
    (* TODO: use queue, put sets to R, which are longest in transversal direction *)
    then do {
      let splits = split_appr optns X;
      ASSUME (list_all (safe_appr optns) splits);
      RETURN (map (\<lambda>X. (h, X))(splits)@XS', PS, RS)
    } else do {
      (h', error, CY, Y) \<leftarrow> (if method_id optns = 2 then
        rk2_step else euler_step_nres) (optns\<lparr>stepsize:=h\<rparr>) X;
      ASSUME (safe_appr optns Y);
      ASSUME (safe_appr optns CY);
      if error > adaptive_atol optns then do {
        RETURN ((h'/2, X)#map (\<lambda>(h,X). ((h/2), X)) XS', PS, RS)
      } else if intersects_appr optns CY sctn \<or> intersects_appr optns Y sctn then (
        if h' > max_intersection_step optns then do {
          RETURN ((h'/2, X)#XS', PS, RS)
        } else do {
          (PS', CXS', RS') \<leftarrow> do_intersection (optns\<lparr>stepsize:=h'\<rparr>) sctn X;
          ASSUME (list_all (\<lambda>CX. safe_appr optns CX) PS');
          ASSUME (list_all (\<lambda>CX. safe_appr optns CX) CXS');
          ASSUME (list_all (\<lambda>CX. safe_appr optns CX) RS');
          RETURN (XS', (X, PS')#PS, RS'@RS)
        }
      ) else do {
        F \<leftarrow> THE_NRES (ode_approx optns [Y]);
        if fst (strongest_transversal F) = normal sctn then
          RETURN (XS'@[(min (3*h'/2) (max h' (truncate_down (precision optns) (h' * real_divl (precision optns) (adaptive_atol optns) (error)))), Y)], PS, RS)
        else
          RETURN (XS', PS, Y#RS)
      }
    }
  })
  (XS, [], [])"

definition "resolve_sctn_impl optns sctn XS =
  WHILE (\<lambda>(XS, PS, RS). XS \<noteq> []) (\<lambda>(XS, PS, RS). do {
    let (h, X) = hd XS;
    let XS' = List.removeAll (h, X) XS;
    let _ = print (STR ''= XS: '');
    let _ = println (int_to_string (length XS));
    let _ = print (STR ''= PS: '');
    let _ = println (int_to_string (length PS));
    let _ = print (STR ''= RS:  '');
    let _ = println (int_to_string (length RS));
    if width_appr optns X > max_tdev_thres optns
    then do {
      let _ = print (STR ''= splitting: '');
      let _ = print_real (width_appr optns X);
      let splits = split_appr optns X;
      ASSUME (list_all (safe_appr optns) splits);
      let _ = print (STR ''-->'');
      let _ = map ((\<lambda>x. let _ = print (STR '' '') in print_real x) o width_appr optns) splits;
      RETURN (map (\<lambda>X. (h, X))splits@XS', PS, RS)
    } else do {
      let _ = print (STR ''= euler step: stepsize = '');
      let _ = print_real h;
      let _ = println (STR '''');
      (* TODO: optimization: do not compute interval bounds except for intersections! *)
      (h', error, CY, Y) \<leftarrow> (if method_id optns = 2 then rk2_step else euler_step_nres)
        (optns\<lparr>stepsize:=h\<rparr>) X;
      ASSUME (safe_appr optns CY);
      ASSUME (safe_appr optns Y);
      let _ = println (STR ''= interpolated step:'');
      let _ = printing_fun optns CY;
      let _ = println (STR ''= discrete step:'');
      let _ = printing_fun optns Y;
      let _ = print (STR ''= error: '');
      let _ = print_real error;
      let _ = println (STR '''');
      if error > adaptive_atol optns then do {
        let _ = println (STR ''= revoking step because of too large error '');
        RETURN ((h'/2, X)#map (\<lambda>(h,X). ((h/2), X)) XS', PS, RS)
      } else if intersects_appr optns CY sctn \<or> intersects_appr optns Y sctn then do {
        let _ = println (STR ''= error OK, intersection '');
        if h' > max_intersection_step optns then do {
          let _ = println (STR ''= revoke_intersection, can try with smaller steps'');
          RETURN ((h'/2, X)#XS', PS, RS)
        } else do {
          let _ = println (STR ''= do_intersection'');
          (PS', CXS', RS') \<leftarrow> do_intersection (optns\<lparr>stepsize:=h'\<rparr>) sctn X;
          ASSUME (list_all (\<lambda>CX. safe_appr optns CX) PS');
          ASSUME (list_all (\<lambda>CX. safe_appr optns CX) CXS');
          ASSUME (list_all (\<lambda>CX. safe_appr optns CX) RS');
          let _ = println (STR ''= done_intersection'');
          let _ = map (\<lambda>X. let _ = println (STR ''= interpolated step:'') in printing_fun optns X) (CXS');
          let _ = map (\<lambda>X. let _ = println (STR ''= discrete step:'') in printing_fun optns X) (RS' @ PS');
          RETURN (XS', (X, PS')#PS, RS'@RS)
        }
      } else do {
        let _ = println (STR ''= error OK, euler stepped '');
        F \<leftarrow> THE_NRES (ode_approx optns [Y]);
        if fst (strongest_transversal F) = normal sctn then
          RETURN (XS'@[(min (3*h'/2) (max h' (truncate_down (precision optns) (h' * real_divl (precision optns) (adaptive_atol optns) (error)))), Y)], PS, RS)
        else
          RETURN (XS', PS, Y#RS)
      }
    }
  })
  (XS, [], [])"

lemma resolve_sctn_impl:
  "resolve_sctn_impl optns sctn XS \<le> \<Down>Id (resolve_sctn optns sctn XS)"
  unfolding resolve_sctn_def resolve_sctn_impl_def
    Let_def[where s="println x" for x]
    Let_def[where s="print x" for x]
    Let_def[where s="print_real x" for x]
  apply refine_rcg
  apply refine_dref_type
  apply (auto simp: neq_Nil_conv)
  done

(* run-time assertions: TODO: why not assert? *)
declare det_assert.transfer_ASSUME_remove[refine_transfer del]
lemma [refine_transfer]: "nres_of (if \<phi> then dRETURN () else dSUCCEED) \<le> ASSUME \<phi>"
  by auto

schematic_lemma resolve_sctn_impl_aux:
  "nres_of ?resolve_sctn_impl \<le> resolve_sctn_impl optns sctn XS"
  unfolding resolve_sctn_impl_def euler_step_nres_def cond_application_beta
  by refine_transfer

concrete_definition resolve_sctn_dres for optns sctn XS uses resolve_sctn_impl_aux
declare resolve_sctn_dres.refine[refine_transfer]
lemmas resolve_sctn_impl_refine[refine_transfer] = order_trans[OF resolve_sctn_dres.refine refine_IdD[OF resolve_sctn_impl]]

definition "resolve_collect optns sctn X0s =
  do {
    (_, PS, RS) \<leftarrow> resolve_sctn optns sctn (map (\<lambda>X. (stepsize optns, X)) X0s);
    if (PS = []) then
      RETURN ([], [], RS)
    else do {
      let PS0 = map fst PS;
      let PS' = concat (map snd PS);
      let _ = print (STR ''= partitioning '');
      let _ = print (int_to_string (length PS'));
      cPS \<leftarrow> partition optns PS';
      let _ = print (STR ''partitioned to '');
      let _ = print (int_to_string (length cPS));
      let _ = println (STR '''');
      let _ = map (printing_fun optns) cPS;
      if 2 * length cPS \<le> length PS'
      then RETURN (cPS, cPS, RS)
      else RETURN (PS0, cPS, RS)
    }
  }"

schematic_lemma resolve_collect_aux: "nres_of ?rresolve_collect \<le> resolve_collect optns sctn X0S"
  unfolding resolve_collect_def
  by refine_transfer

concrete_definition resolve_collect_dres for optns sctn X0S uses resolve_collect_aux
declare resolve_collect_dres.refine[refine_transfer]

definition "max_height optns b Xs = (fold max (map (\<lambda>(X). appr_sup_inner optns X b) (tl Xs)) ((\<lambda>(X). appr_sup_inner optns X b) (hd Xs)))"

definition "next_sections optns d Xs =
  do {
    m \<leftarrow>
      FOREACH (set (remdups Xs)) (\<lambda>x m. do {
        f \<leftarrow> THE_NRES (ode_approx optns [x]);
        RETURN ((x, strongest_transversal f)#m)
      }) [];
    let _ = print (STR ''= next sections with d = '');
    let _ = print_real d;
    let _ = println (STR '''');
    RETURN (map (\<lambda>b.
      let
        xs = map (\<lambda>(x, _, s). (x, s)) (filter (\<lambda>(x, b', s). b = b') m);
        h' = max_height optns b (map fst xs);
        f = fold max (map snd xs) 0;
        h'' = h' + abs (min_section_distance optns) + abs (appr_sup_inner optns (appr_of_ivl ((abs d * f)*\<^sub>Rb) ((abs d * f)*\<^sub>Rb)) b);
        sctn' = List.find (\<lambda>sctn. normal sctn = b \<and> h' \<le> pstn sctn \<and>
          (abs (h'' - pstn sctn) \<le> snap_sections optns \<or> h'' \<ge> pstn sctn)) (stop_sections optns@checkpoint_sections optns);
        sctn = the_default (Sctn b h'') sctn'
      in (map fst xs, sctn)) (remdups (map (fst o snd) m)))
  }"

schematic_lemma next_sections_aux: "nres_of ?next_sections \<le> next_sections optns d Xs"
  unfolding next_sections_def
  using set_iterator_foldli_correct[refine_transfer] distinct_remdups[refine_transfer]
  by refine_transfer

concrete_definition next_sections_dres for optns d Xs uses next_sections_aux
declare next_sections_dres.refine[refine_transfer]

definition "poincare_distance_d optns X0s =
  do {
    (_, _, stops) \<leftarrow>
      WHILE (\<lambda>(unres, cps, stops). unres \<noteq> [] \<or> cps \<noteq> []) (\<lambda>(unres, cps, stops). do {
        if (unres \<noteq> [])
        then do {
          let (XS, sctn) = hd unres;
          let _ = print (STR ''= unres: '');
          let _ = print (int_to_string (length unres));
          let _ = print (STR '' cps: '');
          let _ = print (int_to_string (length cps));
          let _ = print (STR '' stops: '');
          let _ = print (int_to_string (length stops));
          let _ = print (STR '' XS: '');
          let _ = print (int_to_string (length XS));
          let _ = print (STR '' Flowing towards: '');
          let _ = print_eucl (normal sctn);
          let _ = print (STR '' -- '');
          let _ = print_real (pstn sctn);
          let _ = println (STR '''');
          ASSUME (list_all (safe_appr optns) XS);
          (XS, XSg, RS) \<leftarrow> resolve_collect optns sctn XS;
          let _ = println (STR ''= resolved and collected'');
          RS' \<leftarrow> next_sections optns (next_section_turn_distance_factor optns * stepsize optns) RS;
          if sctn \<in> set (stop_sections optns) then
            let _ = print (STR ''= stopping '')
            in RETURN (tl unres@RS', cps, (XSg, sctn)#stops)
          else if sctn \<in> set (checkpoint_sections optns) then
            let _ = print (STR ''= checkpointing '')
            in RETURN (tl unres@RS', (XSg, sctn)#cps, stops)
          else do {
            let _ = print (STR ''= adding '');
            let _ = print (STR '' XS: '');
            let _ = print (int_to_string (length XS));
            let _ = print (STR '' RS: '');
            let _ = print (int_to_string (length RS));
            let _ = println (STR '''');
            XS' \<leftarrow> next_sections optns (next_section_distance_factor optns * stepsize optns) XS;
            RETURN (tl unres@XS'@RS', cps, stops)
          }
        } else do {
          let _ = print (STR '' We are at a checkpoint, here we can reunite different groups again'');
          let sctns = map snd cps;
          let unres' = map (\<lambda>sctn. (List.maps fst (filter (\<lambda>(_, sctn'). sctn' = sctn) cps),sctn)) sctns;
          unres'' \<leftarrow> FOREACH (set (remdups unres')) (\<lambda>(XS, sctn) unres''. do {
            XS' \<leftarrow> next_sections optns (next_section_distance_factor optns * stepsize optns) XS;
            RETURN (XS'@unres'')}) [];
          RETURN (unres'', [], stops)
        }
      }) ([(X0s, start_section optns)], [], []);
    RETURN stops
  }"

schematic_lemma poincare_distance_d_aux: "nres_of ?poincare_distance_d \<le> poincare_distance_d optns X0s"
  using set_iterator_foldli_correct[refine_transfer] distinct_remdups[refine_transfer]
  unfolding poincare_distance_d_def
  by refine_transfer

concrete_definition poincare_distance_d_dres for optns X0s uses poincare_distance_d_aux
declare poincare_distance_d_dres.refine[refine_transfer]

definition poincare_irects::"_ \<Rightarrow> int list \<Rightarrow> 'b list nres"
  where "poincare_irects optns is = do {
      XS \<leftarrow> poincare_distance_d optns [appr_of_irect optns is];
      RETURN (List.maps (\<lambda>(a, b). map (\<lambda>a. project_appr a (normal b)(pstn b)) a) XS)
    }"

schematic_lemma poincare_irects_aux: "nres_of ?poincare_irects \<le> poincare_irects optns IS"
  using set_iterator_foldli_correct[refine_transfer] distinct_remdups[refine_transfer]
  unfolding poincare_irects_def
  by refine_transfer

concrete_definition poincare_irects_dres for optns IS uses poincare_irects_aux
declare poincare_irects_dres.refine[refine_transfer]

end

declare approximate_ivp0.split_appr_while_dres_def[code]
declare approximate_ivp0.euler_step_dres_def[code]
declare approximate_ivp0.do_intersection_dres_def[code]
declare approximate_ivp0.resolve_sctn_dres_def[code]
declare approximate_ivp0.strongest_direction_def[code]
declare approximate_ivp0.middle_appr_def[code]
declare approximate_ivp0.transversal_directions_def[code]
declare approximate_ivp0.strongest_transversal_def[code]
declare approximate_ivp0.poincare_distance_d_def[code]
declare approximate_ivp0.next_sections_def[code]
declare approximate_ivp0.P_appr_def[code]
declare approximate_ivp0.P_iter.simps[code]
declare approximate_ivp0.cert_stepsize.simps[code]
declare approximate_ivp0.euler_step_def[code]
declare approximate_ivp0.advance_euler.simps[code]
declare approximate_ivp0.euler_series_result_def[code]
declare approximate_ivp0.euler_series_print[code]
declare approximate_ivp0.partition_dres_def[code]
declare approximate_ivp0.poincare_distance_d_dres_def[code]
declare approximate_ivp0.next_sections_dres_def[code]
declare approximate_ivp0.max_height_def[code]
declare approximate_ivp0.resolve_collect_dres_def[code]
declare approximate_ivp0.poincare_irects_dres_def[code]
declare approximate_ivp0.rk2_incr_dres_def[code]
declare approximate_ivp0.rk2_remainder_dres_def[code]
declare approximate_ivp0.rk2_step_dres_def[code]

abbreviation "msum_aform' \<equiv> \<lambda>X. msum_aform (degree_aform X) X"

abbreviation "uncurry_options \<equiv> \<lambda>f x. f (precision x) (tolerance x)"

definition "split_aform_largest_uncond X =
    (let (i, x) = max_pdev (snd X) in split_aform X i)"

definition "split_aform_largest optns X =
  (let (a, b) =
      split_aform_largest_uncond (fst X,
        summarize_threshold (precision optns) (presplit_summary_tolerance optns)
          (degree_aform X) (snd X))
    in [a, b])"

text {* intersection with plane *}

definition inter_aform_plane
  where "inter_aform_plane optns Xs sctn =
    inter_aform_plane_ortho (precision optns) Xs (normal sctn) (pstn sctn)"

definition "reduce_aform optns X = apsnd (summarize_threshold (precision optns) (presplit_summary_tolerance optns) (degree_aform X)) X"

definition aform_inf_inner where "aform_inf_inner optns X n = Inf_aform' (precision optns) (inner_aform X n)"
definition aform_sup_inner where "aform_sup_inner optns X n = Sup_aform' (precision optns) (inner_aform X n)"

definition "support_aform optns b X = (let ia = inner_aform X b in fst X \<bullet> b + tdev' (precision optns) (snd ia))"

definition "width_aform optns X =
  (let t = tdev' (precision optns) (snd X) in listsum' (precision optns) (map (\<lambda>i. t \<bullet> i) Basis_list))"

locale aform_approximate_sets0 =
  approximate_sets0
    aform_of_ivl msum_aform' Affine Joints
    Inf_aform Sup_aform
    "uncurry_options add_aform_componentwise"
    "uncurry_options scaleQ_aform_componentwise"
    "uncurry_options scaleR_aform_ivl"
    split_aform_largest
    support_aform
    aform_inf_inner
    aform_sup_inner
    inter_aform_plane
    reduce_aform
    independent_aform
    width_aform

interpretation aform_approximate_sets0 .

locale aform_approximate_ivp0 =
  approximate_ivp0
    aform_of_ivl msum_aform' Affine Joints
    Inf_aform Sup_aform
    "uncurry_options add_aform_componentwise"
    "uncurry_options scaleQ_aform_componentwise"
    "uncurry_options scaleR_aform_ivl"
    split_aform_largest
    support_aform
    aform_inf_inner
    aform_sup_inner
    inter_aform_plane
    reduce_aform
    independent_aform
    width_aform

interpretation aform_approximate_ivp0 ode ode_d ode_d2 for ode ode_d ode_d2 .

definition print_rectangle
  where
  "print_rectangle m i t0 X =
    (let
      _ = print (STR ''rect: '');
      _ = print (int_to_string i);
      _ = print (STR '': '');
      _ = print_real t0
    in
      if i mod m = 0 then
        let
          R = Radius X;
          _ = print (STR '': '');
          _ = print_eucl (fst X - R);
          _ = print (STR '' - '');
          _ = print_eucl (fst X + R);
          _ = print (STR ''; devs: '');
          _ = print (int_to_string (length (list_of_pdevs (snd X))));
          _ = print (STR ''; width: '');
          _ = print_real (infnorm R);
          _ = print (STR ''; tdev: '');
          _ = print_eucl R;
          _ = print (STR ''; maxdev: '');
          _ = print_eucl (snd (max_pdev (snd X)));
          _ = println (STR '''')
        in ()
      else println (STR ''''))"

definition print_aform::"'a::executable_euclidean_space aform \<Rightarrow> unit"
  where
  "print_aform X =
    (let
      _ = print (STR ''aform('');
      _ = print (int_to_string (length (Basis_list::'a list)));
      _ = print (STR ''): '');
      _ = print_eucl (fst X);
      _ = print (STR '' -- '');
      _ = map (\<lambda>(i, x). print_eucl x) (list_of_pdevs (snd X));
      _ = println (STR '''')
    in ())"

definition "ivls_of_aforms p ress = map (\<lambda>(t0, CX, t1, X).
    (t0, eucl_truncate_down p (Inf_aform CX), eucl_truncate_up p (Sup_aform CX), t1,
      eucl_truncate_down p (Inf_aform X), eucl_truncate_up p (Sup_aform X))) ress"

primrec summarize_ivls where
  "summarize_ivls [] = None"
| "summarize_ivls (x#xs) = (case summarize_ivls xs of
      None \<Rightarrow> Some x
    | Some (t0', cl', cu', t1', xl', xu') \<Rightarrow>
      case x of (t0, cl, cu, t1, xl, xu) \<Rightarrow>
      if t0 = t1' then
        Some (min t0 t0', inf cl cl', sup cu cu', max t1 t1',
          if t1 \<le> t1' then xl' else xl, if t1 \<le> t1' then xu' else xu)
      else None)"

fun "set_res_of_ivl_res"
where "set_res_of_ivl_res (t0, CXl, CXu, t1 ,Xl, Xu) = (t0, {CXl .. CXu}, t1, {Xl .. Xu})"

fun parts::"nat\<Rightarrow>'a list\<Rightarrow>'a list list"
where
  "parts n [] = []"
| "parts 0 xs = [xs]"
| "parts n xs = take n xs # parts n (drop n xs)"

definition summarize_enclosure
where "summarize_enclosure p m xs =
  map the (filter (-Option.is_none) (map summarize_ivls (parts m (ivls_of_aforms p xs))))"

definition "ivls_result p m = (apsnd (summarize_enclosure p m o snd)) o snd"

definition "default_optns =
    \<lparr>
    precision = 53,
    tolerance = FloatR 1 (- 8),
    stepsize  = FloatR 1 (- 8),
    iterations = 40,
    halve_stepsizes = 10,
    widening_mod = 40,
    rk2_param = FloatR 1 (- 1),
    max_tdev_thres = FloatR 1 100,
    presplit_summary_tolerance = FloatR 1 0,
    collect_granularity = FloatR 1 100,
    start_section = Sctn 0 0,
    checkpoint_sections = [Sctn 0 0],
    stop_sections = [Sctn 0 0],
    safe_appr = (\<lambda>_. True),
    snap_sections = 1,
    min_section_distance = 1/2,
    next_section_distance_factor = 50,
    next_section_turn_distance_factor = 0,
    printing_fun = print_aform,
    irect_mod = 8,
    result_fun = ivls_result 23 1,
    max_intersection_step = FloatR 1 (- 8),
    adaptive_atol = FloatR 1 (-12),
    method_id = 2
  \<rparr>"

end

