# Theory Abstract_Reachability_Analysis_C1

```theory Abstract_Reachability_Analysis_C1
imports
Abstract_Reachability_Analysis
"../Refinement/Weak_Set"
"../Refinement/Refine_Parallel"
"../Refinement/Refine_Default"
"../Refinement/Refine_Phantom"
"../Refinement/Refine_ScaleR2"
begin

definition blinfun_of_list :: "real list ⇒ 'a ⇒⇩L 'a::executable_euclidean_space"
where "blinfun_of_list xs = blinfun_of_matrix (λi j. xs ! ((index Basis_list i) * DIM('a) + (index Basis_list j)))"

definition vec1_of_list :: "real list ⇒ 'n::{finite, one, plus} vec1"
where "vec1_of_list xs =
(vector (take CARD('n) xs), vector (map (λi. vector (nths xs {CARD('n)*i..CARD('n)*Suc i})) [1..<Suc (CARD('n))]))"

definition flow1_of_vec1 :: "'n vec1 ⇒ ('n rvec * ('n rvec ⇒⇩L 'n::finite rvec))"
where "flow1_of_vec1 xs = (fst xs, blinfun_of_vmatrix (snd xs))"

definition vec1_of_flow1 :: "('n::finite eucl1) ⇒ 'n vec1"
where "vec1_of_flow1 xs = (fst xs, matrix (snd xs))"

lemma vec1_of_flow1_flow1_of_vec1[simp]:
"vec1_of_flow1 (flow1_of_vec1 X) = X"
unfolding vec1_of_flow1_def flow1_of_vec1_def
by (transfer) (auto simp: matrix_of_matrix_vector_mul)

definition "flow1_of_list xs =
(eucl_of_list (take DIM('a::executable_euclidean_space) xs)::'a,
blinfun_of_list (take (DIM('a)*DIM('a)) (drop DIM('a) xs @
replicate (DIM('a)*DIM('a) - (length xs - DIM('a))) 0))::'a⇒⇩L'a)"

lemma blinfun_of_list_eq_blinfun_of_vmatrix:
assumes "length xs = CARD('n)*CARD('n::enum)"
shows "blinfun_of_list xs = blinfun_of_vmatrix (eucl_of_list xs::((real, 'n) vec, 'n) vec)"
using assms
apply (auto simp: blinfun_of_list_def)
apply (auto intro!: simp: blinfun_ext blinfun_of_vmatrix.rep_eq blinfun_of_matrix.rep_eq)
subgoal for i
apply (subst (2) eucl_of_list_list_of_eucl[symmetric, of i])
apply (subst eucl_of_list_matrix_vector_mult_eq_sum_nth_Basis_list)
by (auto simp: sum_Basis_sum_nth_Basis_list scaleR_sum_left intro!: sum.cong)
done

definition "ghost_rel = Pair () ` UNIV"
consts i_ghost::interface
lemmas [autoref_rel_intf] = REL_INTFI[of ghost_rel i_ghost]
lemma ghost_relI: "((), x) ∈ ghost_rel" by (auto simp: ghost_rel_def)

definition [refine_vcg_def]: "GSPEC x = SPEC x"
lemma [autoref_op_pat_def]: "GSPEC x ≡ Autoref_Tagging.OP (GSPEC x)" by auto

lemma GSPEC_impl[autoref_rules]:
assumes "SIDE_PRECOND (Ex P)"
shows "(RETURN (), GSPEC P) ∈ ⟨ghost_rel⟩nres_rel"
using assms by (auto simp: nres_rel_def ghost_rel_def GSPEC_def intro!: RETURN_SPEC_refine)

context approximate_sets_ode' begin

definition "c1_info_of_appr XD =
(case snd XD of None ⇒ eucl_of_list ` set_of_appr (fst XD) × UNIV
| Some DD ⇒ flow1_of_list ` set_of_appr (fst XD @ DD))"
definition "c1_info_of_apprs x = ⋃(set (map c1_info_of_appr x))"
definition "c1_info_of_appr' x = Affine_Code.the_default UNIV (map_option c1_info_of_apprs x)"
definition "c1_info_of_appre X = scaleR2 (fst (fst X)) (snd (fst X)) (c1_info_of_appr (snd X))"
definition "c1_info_of_apprse x = ⋃(set (map c1_info_of_appre x))"

definition [simp]: "op_image_flow1_of_vec1 = (`) flow1_of_vec1"

definition [simp]: "op_image_flow1_of_vec1_coll = (`) flow1_of_vec1"

definition [simp]: "op_image_fst = (`) fst"

definition [refine_vcg_def]:
"vec1rep CX = SPEC (λR. case R of None ⇒ True | Some X ⇒ X = vec1_of_flow1 ` CX)"

definition [simp]: "op_times_UNIV X = X × UNIV"

definition appr1_rel::"(('b list × 'b list option) ×
('a::executable_euclidean_space c1_info set)) set"
where appr1_rel_internal: "appr1_rel = {((xs, None), X × UNIV) |xs X. (xs, X) ∈ appr_rel} ∪
{((xs, Some ys), X::('a c1_info) set) |xs ys X. X = flow1_of_list ` set_of_appr (xs @ ys) ∧
length xs = DIM('a::executable_euclidean_space) ∧
length ys = DIM('a) * DIM('a)}"

abbreviation "appr1e_rel ≡ ⟨appr1_rel⟩scaleR2_rel"

text ‹TODO: remove ‹...:::relation› from this file›
definition "solve_poincare_plane (n::'n::enum rvec) (CX::'n eucl1 set) = do {
X ← mk_safe (((`) fst CX::'n rvec set));
F ← ode_set (X);
nonzero_component (ST ''solve_poincare_map: not nonzero!'') F n;
let i = index Basis_list n;
ASSERT (i < length solve_poincare_slp);
vCX ← vec1rep CX;
case vCX of None ⇒
do {
RETURN (op_times_UNIV (X))
}
| Some vCX ⇒
do {
(R::'n vec1 set) ← approx_slp_appr (map fold_const_fa (solve_poincare_fas i)) (solve_poincare_slp ! i) (list_of_eucl ` vCX);
let R = (op_image_flow1_of_vec1 (R:::appr_rel):::appr1_rel);
X ← mk_safe (op_image_fst R);
F ← ode_set (X);
nonzero_component (ST ''solve_poincare_map: not nonzero!'') (F) n;
RETURN (R::'n eucl1 set)
}
}"

definition embed1::"'n::enum rvec ⇒ ('n rvec * (real^'n::enum^'n::enum))" where [simp]: "embed1 x = (x, 0)"

definition "choose_step1 (X::'n::enum eucl1 set) h = do {
lX ← vec1rep X;
case lX of None ⇒
do {
sX ← mk_safe (fst ` X);
(h, err, CX', X') ← choose_step sX h;
⌦‹err ← width_spec (err:::appr_rel);›
RETURN (h, err × UNIV, (CX') × UNIV, (X') × UNIV)
}
| Some vX ⇒
do {
sX ← var.mk_safe vX;
(h, err, CX', X') ← var.choose_step sX h;
let CX' = flow1_of_vec1 ` ((CX':::appr_rel):::appr_rel);
let X' = flow1_of_vec1 ` ((X':::appr_rel):::appr_rel);
let err' = flow1_of_vec1 ` (err:::appr_rel);
⌦‹err ← width_spec (fst ` (flow1_of_vec1 ` (err:::appr_rel)):::appr_rel);›
RETURN (h, err', CX', X'::'n eucl1 set)
}
}"

definition "plane1_of x = plane_of x × UNIV"
definition "below1_halfspaces x = below_halfspaces x × UNIV"
definition "sbelow1_halfspaces x = sbelow_halfspaces x × UNIV"
abbreviation "plane1_invar_rel ≡ λA. ⟨(⟨lv_rel⟩sctn_rel), A⟩invar_rel plane1_of "

definition "c1_info_invar N XD ⟷ length (fst XD) = N ∧ (case snd XD of Some DD ⇒ length DD = (length (fst XD))⇧2 | None ⇒ True)"

definition op_image_zerofst :: "('a × 'c) set ⇒ ('a::zero × 'c) set"
where [simp]: "op_image_zerofst ≡ λX. (λx. (0, snd x)) ` X"

definition op_image_zerofst_vec :: "('n::enum vec1) set ⇒ ('n::enum vec1) set"
where [simp]: "op_image_zerofst_vec ≡ λX. (λx. (0, snd x)) ` X"

definition [simp]: "op_image_embed1 X = embed1 ` X"

definition "inter_sctn1_spec (X::'n::enum eucl1 set) (sctn::'n rvec sctn) = do {
R ← inter_sctn_spec (fst ` X) sctn;
vX ← vec1rep X;
case vX of None ⇒
do {
let R1 = R × UNIV;
RETURN (R1, R1)
}
| Some X ⇒
do {
let sctn = ((Sctn (embed1 (normal sctn)) (pstn sctn)::'n vec1 sctn));
R1 ← inter_sctn_spec X sctn;
let R2 = op_image_zerofst_vec X + op_image_embed1 R;
RETURN ((flow1_of_vec1 ` R1), (flow1_of_vec1 ` R2))
}
}"

definition "op_image_fst_coll_nres XS = do {
XSs ← sets_of_coll XS;
FORWEAK XSs (RETURN op_empty_coll) (λX.
RETURN (mk_coll (op_image_fst X:::appr_rel))) (λA B. RETURN (B ∪ A))
}"

lemma op_image_fst_coll_nres_spec[le, refine_vcg]: "op_image_fst_coll_nres X ≤ SPEC (λR. R = fst ` X)"
unfolding op_image_fst_coll_nres_def
by (refine_vcg FORWEAK_mono_rule[where I="λit s. fst ` ⋃it ⊆ s ∧ s ⊆ fst ` X"])
(auto, force+)

definition [simp]: "op_image_fst_coll = (`) fst"

definition "fst_safe_coll XS = do {
C ← op_image_fst_coll_nres XS;
mk_safe_coll (C:::clw_rel appr_rel)
}"

definition [refine_vcg_def]:
"vec1reps X = do {
XS ← sets_of_coll X;
FORWEAK XS (RETURN (Some op_empty_coll))
(λx. do {
x ← vec1rep x;
RETURN (map_option mk_coll x:::⟨clw_rel appr_rel⟩option_rel)
})
(λa b.
case (a, b) of (Some a, Some b) ⇒ RETURN (Some (b ∪ a))
| _ ⇒ RETURN None)
}"

definition "do_intersection_spec S guards ivl sctn X0 = (λ(PS, CXS).
poincare_mapsto {x ∈ ivl. x ∈ plane_of sctn} X0 S CXS PS ∧
CXS ∩ guards = {} ∧
CXS ∩ ivl ∩ plane_of sctn = {} ∧
fst ` PS ∩ guards = {} ∧
fst ` PS ⊆ {x ∈ ivl. x ∈ plane_of sctn} ∧
fst ` PS ∪ CXS ⊆ Csafe ∧
0 ∉ (λx. ode x ∙ normal sctn) ` fst ` PS ∧
(∀x∈PS. (∀⇩F x in at (fst x) within plane_of sctn. x ∈ ivl)))"

abbreviation "inter_sbelows X sctns ≡ mk_inter X (sbelow_halfspaces sctns)"

definition "list_of_appr1 X = fst X @ the_default [] (snd X)"

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

definition "nonzero_component_within ivl sctn PDP = do {
fPDP ← mk_safe (fst ` PDP);
F ← ode_set (fPDP);
nonzero_component (ST ''solve_poincare_map: not nonzero!'') F (normal sctn);
op_eventually_within_sctn (op_image_fst PDP:::appr_rel) sctn ivl
}"

definition "do_intersection_invar guards GUARDS ivl sctn X ≡ λ(X', T, PS, PS2, CXS, intersects, inside).
(inside ⟶
(fst ` X ∩ GUARDS = {} ∧
fst ` X ⊆ sbelow_halfspace sctn ∧
ivl ⊆ plane_of sctn ∧
fst ` X ⊆ CXS ∧
fst ` PS ∪ fst ` PS2 ∪ CXS ∪ fst ` X' ⊆ Csafe ∧
T ⊆ nonneg_reals ∧
(¬intersects ⟶ (fst ` X' ∩ plane_of sctn = {} ∧ T = pos_reals)) ∧
CXS ⊆ (sbelow_halfspace sctn - guards) ∧
X' ⊆ (- guards) × UNIV ∧
fst ` (PS ∪ PS2) ∩ guards = {} ∧
(0 ∉ (λx. ode x ∙ normal sctn) ` fst ` (PS ∪ PS2)) ∧
((∀x∈PS ∪ PS2. (∀⇩F x in at (fst x) within plane_of sctn. x ∈ ivl))) ∧
(∃A B. X = A ∪ B ∧
flowsto A T (CXS × UNIV) (X' ∩ (sbelow_halfspace sctn) × UNIV) ∧
poincare_mapsto {x ∈ ivl. x ∙ normal sctn = pstn sctn} B UNIV CXS PS ∧
poincare_mapsto {x ∈ ivl. x ∙ normal sctn = pstn sctn} B UNIV CXS PS2)))"

definition "list_of_appr1e X = fst (snd X) @ the_default [] (snd (snd X)) @
(let (l, u) = fst X;
roer = (λx. if x = - ∞ then FloatR 1 (-88) else if x = ∞ then FloatR 1 88 else real_of_ereal x)
in
appr_of_ivl ops [roer l] [roer u]
)"

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

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

definition "split_spec_param1 X = do {
(vX) ← vec1rep X;
let D = CARD('n);
case vX of None ⇒ do {
(a, b) ← split_spec_param D (fst ` X::'n::finite rvec set);
RETURN (a × UNIV, b × UNIV)
}
| Some X ⇒ do {
(a, b) ← split_spec_param D X;
RETURN (op_image_flow1_of_vec1 a, op_image_flow1_of_vec1 b)
}
}"

abbreviation iinfo_rel :: "('c × 'a set) set ⇒ ((real × 'c) × 'a::real_normed_vector set) set"
where "iinfo_rel ≡ λs. ⟨rnv_rel, s⟩info_rel"

definition "split_spec_param1e X = do {
((l, u), Y) ← scaleR2_rep X;
(a, b) ← split_spec_param1 Y;
a ← scaleRe_ivl_spec l u a;
b ← scaleRe_ivl_spec l u b;
RETURN (a, b)
}"

definition "reduce_spec1 C X = do {
vX ← vec1rep X;
case vX of None ⇒ do {
X ← reduce_spec C (fst ` X);
RETURN (X × UNIV)
}
| Some vX ⇒ do {
vX ← reduce_spec C vX;
RETURN (flow1_of_vec1 ` vX)
}
}"
definition "reduce_spec1e C X = do {
((l, u), X) ← scaleR2_rep X;
X ← reduce_spec1 C X;
scaleRe_ivl_spec l u X
}"

definition [refine_vcg_def]: "pre_split_reduce_spec (ro::unit) = SPEC (λx::unit. True)"

definition split_under_threshold::"_ ⇒ real ⇒ 'n::enum eucl1 set ⇒ 'n eucl1 set nres" where
"split_under_threshold ro th X = do {
(_, Ys) ← WHILE⇗λ(Xs, Ys). X ⊆ Xs ∪ Ys⇖ (λ(Xs, Ys). ¬ op_coll_is_empty Xs) (λ(Xs, Ys). do {
(X, Xs') ← (split_spec_coll (Xs:::clw_rel (⟨appr1_rel⟩scaleR2_rel)):::⟨⟨appr1_rel⟩scaleR2_rel ×⇩r clw_rel (⟨appr1_rel⟩scaleR2_rel)⟩nres_rel);
w ← width_spec (op_image_fste X:::appr_rel);
if w ≤ th then RETURN (Xs', mk_coll X ∪ Ys)
else do {
ra ← pre_split_reduce_spec ro;
X ← reduce_spec1e ra X;
(a, b) ← split_spec_param1e (X:::⟨appr1_rel⟩scaleR2_rel);
RETURN (mk_coll (a:::⟨appr1_rel⟩scaleR2_rel) ∪ mk_coll (b:::⟨appr1_rel⟩scaleR2_rel) ∪ Xs', Ys)
}
}) (X:::clw_rel (⟨appr1_rel⟩scaleR2_rel), op_empty_coll:::clw_rel (⟨appr1_rel⟩scaleR2_rel));
RETURN Ys
}"

definition "choose_step1e X h = do {
((l, u), X) ← scaleR2_rep X;
(h', error, CY, Y) ← choose_step1 X h;
Y ← scaleRe_ivl_spec l u Y;
RETURN (h', error, fst ` CY, Y)
}"

definition "step_split ro (X::'n::enum eucl1 set) =
do {
ra ← pre_split_reduce_spec ro;
X ← reduce_spec1e ra X;
(a, b) ← split_spec_param1e (X:::appr1e_rel);
_ ← mk_safe (op_image_fste a);
_ ← mk_safe (op_image_fste b);
width_X ← width_spec (op_image_fste X:::appr_rel);
wa ← width_spec (op_image_fste a:::appr_rel);
wb ← width_spec (op_image_fste b:::appr_rel);
let _ = trace_set (ST ''splitting: '' @ show (lfloat10 width_X) @ ST ''-->'' @ show (lfloat10 wa) @
ST '', ''  @ show (lfloat10 wb)) (None::'n::enum eucl1 set option);
RETURN (mk_coll a ∪ mk_coll b)
}"

definition "width_spec_appr1 X = do {
vX ← vec1rep X;
case vX of None ⇒ width_spec (fst ` X:::appr_rel)
| Some vX ⇒ width_spec (vX:::appr_rel)
}"

definition "tolerate_error1 Y E = tolerate_error (fst ` Y) (fst ` E)"

definition "step_adapt_time (X::'n::enum eucl1 set) h =
do {
let X0 = fst ` X;
_ ← mk_safe (X0:::appr_rel);
(h', error, CY, Y) ← choose_step1e X h;
(te, e) ← tolerate_error1 Y error;
let _ = trace_set1 (ST ''discrete time step: stepsize = '' @ show (lfloat10 h)) (None::'n eucl1 set option);
let _ = trace_set1 (ST ''discrete time step: stepsize' = '' @ show (lfloat10 h')) (None::'n eucl1 set option);
let _ = trace_set1 (ST ''error estimation: '' @ show (lfloat10 e)) (None::'n eucl1 set option);
if ¬ te
then do {
let _ = trace_set (ST ''revoking step'') (None::'n eucl1 set option);
RETURN (0, fst ` X, X, 3 * h' / 2 / 2)
} else do {
let _ = trace_set1 (ST ''error OK, step_adapt_time stepped '') (None::'n eucl1 set option);
let _ = trace_set (ST ''interpolated step:'') (Some (CY));
let _ = print_set True CY;
let _ = trace_set1e (ST ''discrete step:'') (Some (Y));
let _ = print_set1e False Y;
rtol ← adaptive_rtol_spec;
method_id ← method_spec;
prec ← precision_spec;
case approx prec (adapt_stepsize_fa rtol method_id e h') []
of Some ivl_h'' ⇒
let h'' = lower ivl_h'';
_ = trace_set1 (ST ''increase step: stepsize = '' @ show (lfloat10 h'')) (None::'n eucl1 set option)
in RETURN (h', CY, Y, 15/2/2/2/2 * h'')
| None ⇒
let _ = trace_set1 (ST ''increase time step (failure): stepsize = '' @ show (lfloat10 h')) (None::'n eucl1 set option)
in RETURN (h', CY, Y, h' * 5 / 2 / 2)
}
}"

definition [simp]: "eq_spec x y = SPEC (λr. r ⟶ x = y)"
lemma [autoref_itype]: "eq_spec ::⇩i A →⇩i A →⇩i ⟨i_bool⟩⇩ii_nres" by simp

definition "select_with_inter ci a = do {
CIs ← (sets_of_coll ci);
As ← sets_of_coll a;
FORWEAK CIs (RETURN op_empty_coll)
(λci. do {
(c, I) ← (get_inter ci);
FORWEAK As (RETURN op_empty_coll)
(λa. do {
b ← eq_spec a c;
if b then RETURN (mk_coll ci)
else RETURN (op_empty_coll)
})
(λCIS CIS'. RETURN (CIS' ∪ CIS))
})
(λCIS CIS'. RETURN (CIS' ∪ CIS))
}"

abbreviation "fst_safe_colle XS ≡ (mk_safe_coll (op_image_fst_colle XS:::clw_rel appr_rel):::⟨clw_rel appr_rel⟩nres_rel)"

definition "do_intersection_body GUARDS ivl sctn h ≡ λ(X, T, PDPS, PDPS2, CXS, _, _).
do {
(_, _, CX', X') ← choose_step1 (X:::appr1_rel) (h:::rnv_rel);
let _ = trace_set1 (ST ''interpolated step during intersection:'') (Some (CX'));
let _ = print_set1 True (CX');
let _ = trace_set1 (ST ''step during intersection:'') (Some (X'));
let _ = print_set1 False (X');
CHECKs (ST ''unnormal intersection'') (abs (normal sctn) ∈ set Basis_list);
CPDP ← solve_poincare_plane (abs (normal sctn)) CX';
let _ = trace_set1 (ST ''CPDP: '') (Some CPDP);
let _ = print_set1 False (CPDP);
(PDP, PDP2) ← inter_sctn1_spec CPDP sctn;
b1 ← disjoints_spec (mk_coll (op_image_fst X')) (GUARDS);
b2 ← disjoints_spec (mk_coll (op_image_fst CX')) (GUARDS);
b3 ← disjoints_spec (mk_coll (op_image_fst PDP)) (GUARDS);
b4 ← disjoints_spec (mk_coll (op_image_fst PDP2)) (GUARDS);
CHECKs (ST ''do_intersection: hitting several planes :('') (b1 ∧ b2 ∧ b3 ∧ b4);
intersects ← op_intersects (op_image_fst X') sctn;
CX's ← mk_safe (op_image_fst CX');
c1 ← nonzero_component_within ivl sctn PDP;
c2 ← nonzero_component_within ivl sctn PDP2;
RETURN (X', pos_reals:::⟨Id⟩phantom_rel, mk_coll PDP ∪ PDPS,
mk_coll PDP2 ∪ PDPS2,
mk_coll (inter_sbelows (CX's:::appr_rel) {sctn}) ∪ CXS, intersects, c1 ∧ c2)
}"

definition "do_intersection guards ivl sctn (X::'n::enum eucl1 set) (h::real) =
do {
ASSUME (closed ivl);
sp ← subset_spec_plane ivl sctn;
sX ← mk_safe (op_image_fst (X:::appr1_rel));
GUARDS ← unintersect_coll guards;
a ← sbelow_sctn (op_image_fst X) sctn;
b ← disjoints_spec (mk_coll (op_image_fst X)) GUARDS;
let inside = sp ∧ a ∧ b; ― ‹this is a bit of a hack: if the ‹ivl› is not subset of the plane,›
― ‹then do not do intersections›
(X, T, PDPS, PDPS2, CXS, intersects, inside) ←
WHILE⇗do_intersection_invar guards GUARDS ivl sctn X⇖
(λ(X, T, PDPS, PDPS2, CXS, intersects, inside). intersects ∧ inside)
(do_intersection_body GUARDS ivl sctn h)
(X, nonneg_reals:::⟨Id⟩phantom_rel, op_empty_coll:::clw_rel appr1_rel::'n eucl1 set,
op_empty_coll:::clw_rel appr1_rel::'n eucl1 set,
mk_coll (inter_sbelows (sX:::appr_rel) {sctn}), True, inside);
a ← above_sctn (op_image_fst X) sctn;
b ← subset_spec_coll (op_image_fst_coll PDPS) ivl;
b2 ← subset_spec_coll (op_image_fst_coll PDPS2) ivl;
RETURN (inside ∧ b ∧ b2 ∧ a, PDPS, PDPS2, CXS)
}"

definition "resolve_step roptns X h = do {
width_X ← width_spec (op_image_fste X:::appr_rel);
mtdt ← max_tdev_thres_spec roptns;
if ¬ width_X ≤ mtdt
then do {
Y ← step_split roptns X;
RETURN (h, fst ` Y, Y, h)
}
else do {
(h0, CY, Y, h') ← step_adapt_time (X::'n::enum eucl1 set) h;
RETURN (h0, mk_coll (fst ` Y) ∪ mk_coll CY, mk_coll Y, h')
}
}"

definition "pre_intersection_step ro X h = do {
mis ← max_intersection_step_spec ro;
if h > mis
then RETURN (with_infos (h/2) (mk_coll X), mk_coll (fst ` X), op_empty_coll:::clw_rel (iinfo_rel appr1e_rel))
else do {
width_X ← width_spec (op_image_fste X:::appr_rel);
pig ← pre_inter_granularity_spec ro;
if width_X ≤ pig then
RETURN (with_infos h (op_empty_coll:::clw_rel appr1e_rel), mk_coll (fst ` X),
with_infos (5 * h / 2 / 2) (mk_coll X))
else do {
X' ← step_split ro X;
RETURN (with_infos h X', fst ` X', op_empty_coll:::clw_rel (iinfo_rel appr1e_rel))
}
}
}"

definition "reach_cont ro (guardsi::'n::enum rvec set) XS0 =
do {
startstep ← start_stepsize_spec;
(_, XS0') ← scaleR2_rep_coll XS0;
sXS0 ← fst_safe_coll XS0';
let fX0 = op_image_fst_colle XS0;
guards ← (unintersect_coll (guardsi:::clw_rel (iplane_rel lvivl_rel)):::⟨clw_rel lvivl_rel⟩nres_rel);
d ← disjoints_spec fX0 (guards);
CHECKs (ST ''reach_cont: starting from guarded set'') d;
(_, CXS, GS) ←
WHILE⇗(λ(XS, CXS, GS).
flowsto XS0 {0..} (CXS × UNIV) (XS ∪ GS) ∧
(XS ∪ GS ∪ CXS × UNIV) ⊆ (Csafe - guards) × UNIV ∧
XS0 ∪ GS ⊆ CXS × UNIV)⇖
(λ(XS, CXS, GS). ¬ op_coll_is_empty XS) (λ(XS, CXS, GS).
do {
(hX, XS') ← (split_spec_exact XS:::⟨iinfo_rel (appr1e_rel) ×⇩r clw_rel (iinfo_rel (appr1e_rel))⟩nres_rel);
(h::real, X) ← get_info hX;
let _ = trace_set1e (ST ''next step in resolve_sctns using'') (Some X);
cXS::nat ← card_info XS;
cGS::nat ← card_info GS;
let _ = trace_set1 (ST ''XS: '' @ show cXS) (None::'n eucl1 set option);
let _ = trace_set1 (ST ''GS: '' @ show cGS) (None::'n eucl1 set option);
(h0, fCX', X', h') ← resolve_step ro X h;
sfCX' ← (mk_safe_coll (fCX':::clw_rel appr_rel):::⟨clw_rel appr_rel⟩nres_rel);
let fX' = (fst ` X');
fXS ← ivls_of_sets (fCX' ∪ fX');
IS ← inter_overappr guards fXS;
let d = op_coll_is_empty IS;
if d then RETURN (with_infos h' X' ∪ XS':::clw_rel (iinfo_rel appr1e_rel), sfCX' ∪ CXS, GS)
else do {
(hX', fCX', hG') ← pre_intersection_step ro X h0;
sfCX' ← (mk_safe_coll (fCX':::clw_rel appr_rel):::⟨clw_rel appr_rel⟩nres_rel);
_ ← fst_safe_colle (uninfo hX');
_ ← fst_safe_colle (uninfo hG');
fGs ← ivls_of_sets (op_image_fst_colle (uninfo hG') ∪ fCX' ∪ op_image_fst_colle (uninfo hX'));
d ← disjoints_spec (sets_of_ivls guards) fGs;
CHECKs (ST ''reach_cont: pre_intersection_step should not change disjointness condition!'') d;
iguards ← select_with_inter guardsi IS;
iG' ← with_coll_infos iguards hG';
RETURN (hX' ∪ XS', sfCX' ∪ CXS, iG' ∪ GS)
}
})
(with_infos startstep (XS0:::clw_rel appr1e_rel):::clw_rel (iinfo_rel appr1e_rel),
sXS0:::clw_rel appr_rel, op_empty_coll:::clw_rel (⟨iplane_rel (lvivl_rel::(_×'n rvec set)set), iinfo_rel appr1e_rel⟩info_rel));
RETURN (CXS, GS)
}"

definition "reach_cont_par roptns guards XS = do {
XS ← sets_of_coll XS;
PARS ← PAR_IMAGE (λX (CX, G).
G ∪ (CX × UNIV) ⊆ (Csafe - guards) × UNIV ∧
X ∪ G ⊆ CX × UNIV ∧ flowsto X {0..} (CX × UNIV) G)
(λX. reach_cont roptns guards (mk_coll X)) XS;
RETURN (⋃(fst ` snd ` PARS), ⋃(snd ` snd ` PARS))
}"

definition "subset_iplane_coll x ics = do {
X ← unintersect x;
ics ← sets_of_coll ics;
FORWEAK ics (RETURN False) (λic. do {
(i, c) ← get_inter ic;
sctn ← get_plane c;
b1 ← subset_spec_plane X sctn;
RETURN (b1 ∧ op_subset_ivl X i)
}) (λb c. RETURN (b ∨ c))
}"

definition "subsets_iplane_coll xs ics = FORWEAK xs (RETURN True) (λx. subset_iplane_coll x ics) (λa b. RETURN (a ∧ b))"

definition "stable_set p = {x. {0..} ⊆ existence_ivl0 x ∧ (flow0 x ⤏ p) at_top}"

definition symstart_coll::"('n::enum eucl1 set ⇒ ('n rvec set × 'n eucl1 set)nres) ⇒
'n eucl1 set ⇒ ('n rvec set × 'n eucl1 set)nres"
where
"symstart_coll symstart X0 = do {
_ ← (fst_safe_colle (X0:::clw_rel appr1e_rel):::⟨clw_rel appr_rel⟩nres_rel);
X0s ← sets_of_coll X0;
(CX1, X1) ← FORWEAK X0s (RETURN (op_empty_coll, op_empty_coll)) (symstart)
(λ(CX, X) (CX', X'). RETURN (CX' ∪ CX, X' ∪ X));
RETURN (CX1, X1)
}"

definition reach_cont_symstart ::
"_ ⇒ _ ⇒ 'n::enum rvec set
⇒ 'n eucl1 set ⇒ ('n rvec set × 'n eucl1 set) nres"
where "reach_cont_symstart ro symstart (guards::'n rvec set) X0 = do {
let fX0 = op_image_fst_colle X0;
GUARDS ← unintersect_coll guards;
d ← disjoints_spec fX0 GUARDS;
CHECKs (ST ''reach_cont_symstart: starting from guarded set'') d;
(CY, Y0) ← symstart_coll symstart X0;
sCY ← (mk_safe_coll (op_image_fst_colle X0 ∪ CY:::clw_rel appr_rel):::⟨clw_rel appr_rel⟩nres_rel);
b ← disjoints_spec (op_image_fst_colle Y0 ∪ CY) GUARDS;
CHECKs ''reach_cont_symstart with a stupid guard'' b;
(CX, GS) ← (reach_cont_par ro guards Y0:::⟨clw_rel appr_rel ×⇩r clw_rel (⟨iplane_rel lvivl_rel::(_ × 'n rvec set) set, iinfo_rel appr1e_rel⟩info_rel)⟩nres_rel);
let CZ = sCY ∪ CX;
RETURN (CZ, GS)
}"

context includes autoref_syntax begin― ‹TODO: should not be annotating relations here›
definition reach_conts ::
"_ ⇒ _ ⇒ _ ⇒ 'n::enum rvec set
⇒ 'n eucl1 set ⇒ ('n rvec set × ('n rvec set × 'n eucl1 set) set × ('n eucl1 set ⇒ 'n eucl1 set)) nres"
where "reach_conts ro symstart trap (guardsi::'n rvec set) X0 = do {
(CX, GS) ← (reach_cont_symstart ro (symstart:::appr1e_rel → ⟨clw_rel appr_rel ×⇩r clw_rel appr1e_rel⟩nres_rel) guardsi X0:::
⟨clw_rel appr_rel ×⇩r
clw_rel (⟨iplane_rel lvivl_rel::(_×'n rvec set) set, iinfo_rel appr1e_rel⟩info_rel)⟩nres_rel);
(IGSs:: ('n rvec set × 'n eucl1 set) set) ← explicit_info_set GS;
let GSs = snd ` IGSs;
ASSUME (finite GSs);
CHECKs '''' (GSs ≠ {});
ASSERT       (∃f. X0 = ⋃(f ` GSs) ∧ (∀G ∈ GSs. flowsto (f G - trap × UNIV) {0..} (CX × UNIV) (G)));
X0f ← GSPEC (λf. X0 = ⋃(f ` GSs) ∧ (∀G ∈ GSs. flowsto (f G - trap × UNIV) {0..} (CX × UNIV) (G)));
let K = (fst ` IGSs);
b ← subsets_iplane_coll K guardsi;
CHECKs (ST ''reach_conts: subsets_iplane_coll'') b;
RETURN (CX, IGSs:::⟨iplane_rel lvivl_rel ×⇩r clw_rel (iinfo_rel appr1e_rel)⟩list_wset_rel, X0f)
}"
end

definition [refine_vcg_def]: "get_sctns X = SPEC (λR. X = below_halfspaces R)"

definition "leaves_halfspace S X = do {
sctns ← get_sctns S;
sctnss ← op_set_to_list sctns;
(case sctnss of
[] ⇒ RETURN None
| [sctn] ⇒
do {
(Xi, Xs) ← ivl_rep_of_set_coll X;
ASSERT (Xi ≤ Xs);
b ← subset_spec_plane ({Xi .. Xs}:::lvivl_rel) sctn;
CHECKs (ST ''leaves_halfspace: not subset of plane'') b;
F ← ode_set ({Xi .. Xs}:::appr_rel);
sF ← Sup_inner F (normal sctn);
CHECKs (ST ''leaves_halfspace: not down from plane'') (sF < 0);
RETURN (Some sctn)
}
| _ ⇒ do {CHECKs (ST ''leaves_halfspace: not a good halfspace'') False; SUCCEED})
}"

definition "poincare_start_on guards sctn X0S =
do {
X0SS ← sets_of_coll X0S;
(FORWEAK X0SS (RETURN (op_empty_coll:::clw_rel appr1e_rel, op_empty_coll:::clw_rel appr_rel)) (λX0. do {
mk_safe (fst ` X0);
startstep ← start_stepsize_spec;
(h, err, CX1, X1) ← choose_step1e X0 (startstep);
let _ = trace_set (ST ''interpolated start step:'') (Some CX1);
let _ = print_set True CX1;
let _ = trace_set1e (ST ''discrete start step:'') (Some X1);
let _ = print_set1e False X1;
let fX1 = op_image_fste X1;
c0 ← below_sctn (op_image_fste X0) (sctn);
c1 ← sbelow_sctn (fX1) (sctn);
c2 ← disjoints_spec (mk_coll (fX1)) guards;
c3 ← disjoints_spec (mk_coll (CX1)) guards;
mk_safe (fX1);
mk_safe (CX1);
D ← (ode_set (CX1):::⟨appr_rel⟩nres_rel);
d ← Sup_inner D (normal sctn);
let _ = trace_set (ST ''poincare_start_on: D '') (Some D);
CHECKs (ST ''poincare_start_on: is away and really moves away'') (c0 ∧ c1 ∧ c2 ∧ c3 ∧ d < 0);
RETURN (mk_coll X1:::clw_rel appr1e_rel, (mk_coll CX1):::clw_rel appr_rel)
})
(λ(X1, CX1) (X1S, CX1S). RETURN (op_union_coll X1 X1S:::clw_rel appr1e_rel, op_union_coll CX1 CX1S:::clw_rel appr_rel)))
}"

definition [simp]: "isets_of_iivls x = x"

abbreviation "inter_plane A B ≡ mk_inter A (plane_of B)"

definition "do_intersection_core guards ivl sctn hX =
do {
(h, eX) ← get_info hX;
((l, u), X) ← scaleR2_rep eX;
(b, PS1, PS2, CXS) ← do_intersection (guards:::clw_rel (iplane_rel lvivl_rel)) ivl sctn X h;
if b then do {
PS1e ← scaleRe_ivl_coll_spec l u PS1;
PS2e ← scaleRe_ivl_coll_spec l u PS2;
RETURN (PS1e, PS2e, CXS, op_empty_coll)
}
else RETURN (op_empty_coll, op_empty_coll, op_empty_coll, mk_coll eX)
}"

definition "do_intersection_coll guards ivl sctn X =
do {
Xs ← sets_of_coll X;
CHECKs ''nonempty inter nonfinite: '' (Xs ≠ {});
RS ← PAR_IMAGE (λX. λ(P1, P2, CX, X0s).
do_intersection_spec UNIV guards ivl sctn (X - X0s) (P1, CX) ∧
do_intersection_spec UNIV guards ivl sctn (X - X0s) (P2, CX)
∧ fst ` (X - X0s) ⊆ CX
∧ X0s ⊆ X) (do_intersection_core guards ivl sctn) Xs;
ASSUME (finite RS);
RETURN ((⋃(X, (P1, P2, CX, X0s))∈RS. P1),
(⋃(X, (P1, P2, CX, X0s))∈RS. P2),
(⋃(X, (P1, P2, CX, X0s))∈RS. CX),
(⋃(X, (P1, P2, CX, X0s))∈RS. X0s))
}"

definition "op_enlarge_ivl_sctn ivl sctn d = do {
(l, u) ← ivl_rep ivl;
CHECKs ''op_enlarge_ivl_sctn: trying to shrink'' (d ≥ 0);
CHECKs ''op_enlarge_ivl_sctn: empty ivl'' (l ≤ u);
CHECKs ''op_enlarge_ivl_sctn: not in Basis'' (abs (normal sctn) ∈ set Basis_list);
let dOne = sum_list (map (λi. d *⇩R i) Basis_list) - d *⇩R abs (normal sctn);
ASSERT (l - dOne ≤ u + dOne);
RETURN (op_atLeastAtMost_ivl (l - dOne) (u + dOne))
}"

definition "guardset guards = Union (case_prod (∩) ` (λ(x, y). (x, plane_of y)) ` guards)"

definition "resolve_ivlplanes (guards::'n::enum rvec set)
(ivlplanes::'n::enum rvec set)
XS
=
FORWEAK XS (RETURN ({}))
(λ(ivlplane, X).
do {
(ivl, plane) ← (get_inter (ivlplane));
ASSUME (closed ivl);
sctn ← (get_plane plane);
b ← subset_iplane_coll ivlplane ivlplanes;
CHECKs ''reach_conts: subsets_iplane_coll'' b;
(PS1, PS2, CXS, RS) ← (do_intersection_coll guards ivl sctn X);
RETURN {(uninfo X, PS1, PS2, RS, ivl, sctn, CXS)}
})
(λ(PS) (PS'). RETURN (PS' ∪ PS))"

context includes autoref_syntax begin
definition "poincare_onto ro ― ‹options›
symstart trap ― ‹symbolic start and trap›
(guards::'n::enum rvec set) ― ‹avoiding guards›
(ivlplanes::'n::enum rvec set) ― ‹target sections›
(XS0::'n eucl1 set)
(CXS0::'n rvec set)
=
do {
(CXS, XS, X0s) ← (reach_conts ro (symstart:::appr1e_rel → ⟨clw_rel appr_rel ×⇩r clw_rel appr1e_rel⟩nres_rel) trap (ivlplanes ∪ guards) XS0:::
⟨clw_rel appr_rel ×⇩r ⟨iplane_rel lvivl_rel ×⇩r clw_rel (iinfo_rel appr1e_rel)⟩list_wset_rel ×⇩r
ghost_rel⟩nres_rel);
PS ← resolve_ivlplanes guards ivlplanes XS;
_ ← mk_safe_coll CXS0;
RETURN ((λ(X, P1, P2, R, ivl, sctn, CX). (X, P1, P2, R, ivl, sctn, CX, CXS ∪ CXS0)) ` PS)
}"

definition "empty_remainders PS =
FORWEAK PS (RETURN True) (λ(X, P1, P2, R, ivl, sctn, CX, CXS). do { e ← isEmpty_spec R; RETURN e})
(λa b. RETURN (a ∧ b))"

definition [simp]: "empty_trap = {}"

definition empty_symstart::"((real, 'a) vec × (real, 'a) vec ⇒⇩L (real, 'a) vec) set
⇒ ((real, 'a) vec set × ((real, 'a) vec × (real, 'a) vec ⇒⇩L (real, 'a) vec) set) nres"
where [simp]: "empty_symstart ≡ λX. RETURN (op_empty_coll, mk_coll X)"

definition "poincare_onto_empty ro ― ‹options›
(guards::'n::enum rvec set) ― ‹avoiding guards›
(ivlplanes::'n::enum rvec set) ― ‹target sections›
(XS0::'n eucl1 set) =
poincare_onto ro (OP empty_symstart:::appr1e_rel → ⟨clw_rel appr_rel ×⇩r clw_rel appr1e_rel⟩nres_rel)
empty_trap guards ivlplanes XS0"

definition "poincare_onto2 ro ― ‹options›
symstart trap ― ‹symbolic start and trap›
(guards::'n::enum rvec set) ― ‹avoiding guards›
(ivlplanes::'n::enum rvec set) ― ‹target sections›
(XS0::'n eucl1 set) =
do {
(PS) ← (poincare_onto ro (symstart:::appr1e_rel → ⟨clw_rel appr_rel ×⇩r clw_rel appr1e_rel⟩nres_rel)
trap guards ivlplanes XS0 op_empty_coll:::
⟨⟨clw_rel appr1e_rel ×⇩r clw_rel appr1e_rel ×⇩r clw_rel appr1e_rel ×⇩r clw_rel appr1e_rel ×⇩r lvivl_rel ×⇩r ⟨lv_rel⟩sctn_rel ×⇩r
clw_rel (isbelows_rel appr_rel) ×⇩r clw_rel appr_rel⟩list_wset_rel⟩nres_rel);
(PS2) ← FORWEAK PS (RETURN ({})) (λ(X, P1, P2, R, ivl, sctn, CX, CXS).
if op_coll_is_empty R then RETURN ({})
else do {
ivlplaness ← (sets_of_coll ivlplanes:::⟨⟨iplane_rel lvivl_rel⟩list_wset_rel⟩nres_rel);
ivlplaness' ← op_set_ndelete (mk_inter ivl (plane_of sctn)) ivlplaness;
let ivlplanes' = (⋃(mk_coll ` ivlplaness':::⟨clw_rel (iplane_rel lvivl_rel)⟩list_wset_rel));
PS' ← (poincare_onto_empty ro (guards) ivlplanes' R CXS:::
⟨⟨clw_rel appr1e_rel ×⇩r clw_rel appr1e_rel ×⇩r clw_rel appr1e_rel ×⇩r clw_rel appr1e_rel ×⇩r lvivl_rel ×⇩r ⟨lv_rel⟩sctn_rel
×⇩r clw_rel (isbelows_rel appr_rel) ×⇩r clw_rel appr_rel⟩list_wset_rel⟩nres_rel);
b ← empty_remainders PS';
CHECKs (ST ''poincare_onto2: empty remainders!'') b;
ASSUME (finite PS');
RETURN PS'
}) (λPS PS'. RETURN (PS' ∪ PS));
RETURN (Pair True ` PS2 ∪ Pair False ` PS)
}"

definition "width_spec_ivl M x = do {
(i, s) ← ivl_rep x;
RETURN (∑(i, s)←zip (take M (list_of_eucl i)) (take M (list_of_eucl s)). abs (s - i))
}"

definition partition_ivl::"_ ⇒ 'a::executable_euclidean_space set ⇒ 'a::executable_euclidean_space set nres"
where
"partition_ivl roptns xs = (if op_coll_is_empty xs then RETURN (op_empty_coll:::clw_rel lvivl_rel) else do {
(i, s) ← ivl_rep_of_set_coll (sets_of_ivls (xs:::clw_rel lvivl_rel):::clw_rel appr_rel);
ASSERT (i ≤ s);
let r = (op_atLeastAtMost_ivl i s);
(rs, ps) ←
WHILE⇗(λ(rs, ps). (xs) ⊆ rs ∪ ps)⇖ (λ(rs, ps). ¬ op_coll_is_empty (rs:::clw_rel lvivl_rel))
(λ(rs, ps).
do {
(r, rs') ← (split_spec_exact rs:::⟨lvivl_rel ×⇩r clw_rel lvivl_rel⟩nres_rel);
(ri, rs) ← ivl_rep r;
CHECKs (ST ''partition_ivl with strange ivl'') (ri ≤ rs);
width ← width_spec ({ri .. rs}:::appr_rel);
pig ← post_inter_granularity_spec roptns;
if width ≤ pig then
RETURN (rs', mk_coll r ∪ ps)
else do {
(a, b) ← split_spec_ivl (DIM('a)) r;
let isa = (op_inter_ivl_coll (xs:::clw_rel lvivl_rel) (a:::lvivl_rel));
let isb = (op_inter_ivl_coll(xs:::clw_rel lvivl_rel) (b:::lvivl_rel));
ra' ← (if op_coll_is_empty isa then RETURN op_empty_coll else do {
(i', s') ← ivl_rep_of_set_coll (sets_of_ivls isa);
RETURN (mk_coll (({i' .. s'}:::lvivl_rel) ∩ a))
});
rb' ← (if op_coll_is_empty isb then RETURN op_empty_coll else do {
(i', s') ← ivl_rep_of_set_coll (sets_of_ivls isb);
RETURN (mk_coll (({i' .. s'}:::lvivl_rel) ∩ b))
});
RETURN (ra' ∪ rb' ∪ rs', ps)
}
}) (mk_coll r:::clw_rel lvivl_rel, op_empty_coll :::clw_rel lvivl_rel);
RETURN ps
})"

definition
"vec1repse X = do {
XS ← sets_of_coll X;
FORWEAK XS (RETURN (Some op_empty_coll))
(λx. do {
((l, u), x) ← scaleR2_rep x;
xo ← vec1rep x;
case xo of None ⇒ RETURN None
| Some x ⇒ do {
xe ← scaleRe_ivl_spec l u x;
RETURN (Some (mk_coll xe))
}
})
(λa b.
case (a, b) of (Some a, Some b) ⇒ RETURN (Some (b ∪ a))
| _ ⇒ RETURN None)
}"

abbreviation "appre_rel ≡ ⟨appr_rel⟩scaleR2_rel"

definition "scaleR2_rep1 (Y::('a::executable_euclidean_space×_) set) = do {
let D = DIM('a);
((l, u), X) ← scaleR2_rep Y;
(i, s) ← op_ivl_rep_of_set X;
let mig = inf (abs i) (abs s);
CHECKs (ST ''scaleR2_rep1: strange'') (i ≤ s);
(N::real set) ← approx_slp_appr [floatarith.Inverse (norm2⇩e D)] (norm2_slp D) (list_of_eucl ` ({mig .. mig}:::appr_rel));
(sl, su) ← op_ivl_rep_of_set (N:::appr_rel);
let scale = (rnv_of_lv sl + rnv_of_lv su)/2;
CHECKs (ST ''scaleR2_rep1: scale 0'') (scale > 0);
CHECKs (ST ''scaleR2_rep1: l 0'') (l > 0);
CHECKs (ST ''scaleR2_rep1: u 0'') (u > 0);
precision ← precision_spec;
let scalel = real_divl (precision) 1 scale;
let scaleu = real_divr (precision) 1 scale;
CHECKs (ST ''scaleR2_rep1: scalel 0'') (scalel > 0);
CHECKs (ST ''scaleR2_rep1: scaleu 0'') (scaleu > 0);
(i, s) ← op_ivl_rep_of_set X;
let (i0, i1) = split_lv_rel i;
let (s0, s1) = split_lv_rel s;
scaleRe_ivl_spec (scalel * l) (scaleu * u)
(op_atLeastAtMost_ivl (Pair_lv_rel i0 (scale *⇩R i1)) (Pair_lv_rel s0 (scale *⇩R s1)))
}"

definition "ivlse_of_setse X = do {
Xs ← sets_of_coll X;
FORWEAK Xs (RETURN op_empty_coll) (λX. do {
I ← scaleR2_rep1 X;
I ← reduces_ivle I;
RETURN (mk_coll I)
}) (λX' X. RETURN (X' ∪ X))
}"

definition [simp]: "op_image_flow1_of_vec1_colle ≡ op_image_flow1_of_vec1_coll"

definition "okay_granularity ro r = do {
(ri, rs) ← ivl_rep r;
CHECKs (ST ''partition_ivle with strange ivl'') (ri ≤ rs);
width ← width_spec ({ri .. rs}:::appr_rel);
pig ← post_inter_granularity_spec ro;
RETURN (width≤pig)
}"

definition partition_set::"unit ⇒ 'n::enum eucl1 set ⇒ 'n eucl1 set nres"
where
"partition_set ro xs =
(if op_coll_is_empty xs then RETURN (op_empty_coll:::clw_rel appr1e_rel) else do {
ASSERT (xs ≠ {});
pcg ← pre_collect_granularity_spec ro;
xs ← split_under_threshold ro pcg
(xs:::clw_rel appr1e_rel);
vxs ← vec1repse xs;
case vxs of None ⇒ do {
xs ← ivls_of_sets (op_image_fst_colle xs);
ps ← partition_ivl ro xs;
scaleRe_ivl_coll_spec 1 1 (sets_of_ivls ps × UNIV:::clw_rel appr1_rel)
}
| Some xs ⇒ do {
xs ← ivlse_of_setse xs;
ps ← (OP partition_ivle) \$ ro \$ xs;
ps ← setse_of_ivlse ps;
RETURN (op_image_flow1_of_vec1_colle ps)
}
})"

definition partition_sets::"unit ⇒
(bool × 'a::enum eucl1 set × 'a::enum eucl1 set × 'a::enum eucl1 set × 'a::enum eucl1 set × 'a rvec set × 'a rvec sctn × 'a rvec set × 'a rvec set) set ⇒
'a eucl1 set nres"
where
"partition_sets ro xs =
FORWEAK xs (RETURN op_empty_coll) (λ(b, X, PS1, PS2, R, ivl', sctn', CX, CXS). do {
PS ← partition_set ro PS1;
RETURN PS
})
(λa b. RETURN (b ∪ a))"

definition "ivlsctn_to_set xs = (⋃(ivl, sctn)∈set xs. ivl ∩ plane_of sctn)"

definition [refine_vcg_def]: "singleton_spec X = SPEC (λx. X = {x})"

primrec poincare_onto_series where
"poincare_onto_series interrupt trap [] XS0 ivl sctn ro = do {
let guard0 = mk_coll (mk_inter ivl (plane_of sctn));
ASSUME (closed guard0);
XS1 ← (poincare_onto2 (ro:::Id)
(interrupt:::appr1e_rel → ⟨clw_rel appr_rel ×⇩r clw_rel appr1e_rel⟩nres_rel) trap
(op_empty_coll:::clw_rel (iplane_rel lvivl_rel)) guard0 XS0:::
⟨⟨bool_rel ×⇩r clw_rel appr1e_rel ×⇩r clw_rel appr1e_rel ×⇩r clw_rel appr1e_rel ×⇩r clw_rel appr1e_rel ×⇩r lvivl_rel ×⇩r ⟨lv_rel⟩sctn_rel ×⇩r
clw_rel (isbelows_rel appr_rel) ×⇩r clw_rel appr_rel⟩list_wset_rel⟩nres_rel);
(b, X, PS1, PS2, R, ivl', sctn', CX, CXS) ← singleton_spec XS1;
CHECKs (ST ''poincare_onto_series: last return!'') (ivl' = ivl ∧ sctn' = sctn);
RETURN PS2
}"
| "poincare_onto_series interrupt trap ((guardro)#guards) XS0 ivl sctn ro0 = (case guardro of (guard, ro) ⇒
do {
ASSUME (closed ivl);
let guard0 = mk_coll (mk_inter ivl (plane_of sctn));
ASSUME (closed guard0);
ASSUME (∀(guard, ro) ∈ set (guardro#guards). closed guard);
let guardset = (⋃(guard, ro)∈set ((guard0, ro0)#guards). guard);
XS1 ← (poincare_onto2 ro (interrupt:::appr1e_rel → ⟨clw_rel appr_rel ×⇩r clw_rel appr1e_rel⟩nres_rel) trap (guardset:::clw_rel (iplane_rel lvivl_rel)) guard XS0 :::
⟨⟨bool_rel ×⇩r clw_rel appr1e_rel ×⇩r clw_rel appr1e_rel ×⇩r clw_rel appr1e_rel ×⇩r clw_rel appr1e_rel ×⇩r lvivl_rel ×⇩r ⟨lv_rel⟩sctn_rel ×⇩r
clw_rel (isbelows_rel appr_rel) ×⇩r clw_rel appr_rel⟩list_wset_rel⟩nres_rel);
ASSUME (∀(b, X, PS1, PS1, R, ivl, sctn, CX, CXS) ∈ XS1. closed ivl);
XS2 ← partition_sets ro XS1;
_ ← fst_safe_colle XS2;
XS3 ← poincare_onto_series interrupt trap guards XS2 ivl sctn (ro0:::Id);
RETURN XS3
})"

definition "poincare_onto_from interrupt trap
S                      ― ‹leaving this (half)space in the beginning›
(guards)                    ― ‹avoiding guards›
(ivl::'n rvec set)          ― ‹onto ‹ivl››
sctn                        ― ‹which is part of ‹sctn››
ro
(XS0::'n::enum eucl1 set) =
do {
ASSUME (closed ivl);
let guardset = (⋃(ivlsctn, ro)∈set (guards). ivlsctn:::clw_rel (iplane_rel lvivl_rel));
lsctn ← leaves_halfspace S (op_image_fst_colle XS0);
XS0 ← (case lsctn of
None ⇒ RETURN XS0
| Some lsctn =>
do {
CHECKs (ST ''poincare_onto_from: section only makes sense if start section = end section'')
(lsctn = sctn ∨ normal lsctn = - normal sctn ∧ pstn lsctn = - pstn sctn);
guards ← unintersect_coll guardset;
b ← subset_spec_coll (op_image_fst_colle XS0) ivl;
CHECKs (ST ''poincare_onto_from: section only makes sense if we start from there'') b;
(XS0, _) ← poincare_start_on guards lsctn (XS0);
RETURN XS0
}
);
PS ← poincare_onto_series interrupt trap guards XS0 ivl sctn ro;
RETURN PS
}"

definition "subset_spec1 R P dP = do {
R1 ← vec1rep R;
dP ← default_rep UNIV dP;
case (R1, dP) of (_, None) ⇒
op_subset (fst ` R) P
| (Some RdR, Some dP) ⇒
op_subset RdR (P × dP)
| (None, Some _) ⇒ RETURN False
}"

definition "subset_spec1_coll R P dP = do {
XS ← sets_of_coll R;
WEAK_ALL (λx. x ⊆ flow1_of_vec1 ` (P × dP)) XS (λX. subset_spec1 X P dP)
}"

definition "one_step_until_time X0 (ph::bool) (t1::real) =
do {
CHECKs ''one_step_until_time optns'' (0 ≤ t1);
startstep ← start_stepsize_spec;
rk2param ← rk2_param_spec;
let fX0 = fst ` X0;
mk_safe (fX0);
(t, _, X, CX) ← WHILE (λ(t, _, _, _). t < t1) (λ(t, h, X, CXs). do {
let _ = trace_set1e (ST ''choose step from:'') (Some X);
(h0, CX, X, h') ← step_adapt_time X (min h (t1 - t));
CHECKs (ST ''one_step negative step'') (h0 ≥ 0 ∧ h' > 0 ∧ h0 ≤ min h (t1 - t));
let _ = trace_set (ST ''interpolated step:'') (Some CX);
let _ = print_set True CX;
let _ = trace_set1e (ST ''step:'') (Some X);
let _ = print_set1e False X;
let fCX = CX;
mk_safe fCX;
let fX = fst ` X;
mk_safe fX;
RETURN (t + h0, h', X, mk_phantom (mk_coll CX) ∪ CXs)
}) (0::real, startstep, X0, op_union_phantom (mk_phantom (mk_coll fX0)) (op_empty_phantom ph));
RETURN (X, CX)
}"

definition "ivl_of_eucl_coll CY =
do {
(i, s) ← ivl_rep_of_set_coll CY;
ASSERT (i ≤ s);
RETURN (({i .. s}×UNIV):::appr1_rel)
}"

definition "one_step_until_time_ivl X0 (ph::bool) (t1::real) (t2::real) =
do {
(X, CX) ←  one_step_until_time X0 ph t1;
CHECKs (ST ''one_step_until_time_ivl empty time interval'') (0 ≤ t1 ∧ t1 ≤ t2);
(if t2 = t1 then RETURN (X, CX)
else do {
(Y, CYp) ← one_step_until_time X False (t2 - t1);
CY ← get_phantom CYp;
R ← ivl_of_eucl_coll CY;
mk_safe (fst ` R);
R ← scaleRe_ivl_spec 1 1 R;
RETURN (R, CYp ∪ CX)
})
}"

definition "c1_info_invare n X = (let l = (fst (fst X)); u = (snd (fst X))
in (c1_info_invar n (snd X)) ∧ (l < u ∨ -∞ < l ∧ l ≤ u ∧ u < ∞))"

definition "c0_info_of_appr X = eucl_of_list ` set_of_appr X"
definition "c0_info_of_apprs X = (⋃x∈set X. c0_info_of_appr x)"

definition "c0_info_of_appr' X = the_default UNIV (map_option c0_info_of_apprs X)"

lemma the_default_eq: "the_default a x = (case x of None ⇒ a | Some b ⇒ b)"
by (auto split: option.splits)

definition "poincare_onto_from_in_ivl interrupt trap
S                      ― ‹leaving this (half)space in the beginning›
(guards)                    ― ‹avoiding guards›
(ivl::'n rvec set)          ― ‹onto ‹ivl››
sctn                        ― ‹which is part of ‹sctn››
ro
(XS0::'n::enum eucl1 set)
P dP =
do {
RS ← poincare_onto_from interrupt trap S guards ivl sctn ro XS0;
((l, u), R) ← scaleR2_rep_coll RS;
CHECKs (ST ''poincare_onto_from_in: there should not be scaleR2'') (l = 1 ∧ u = 1);
(l, u) ← ivl_rep P;
CHECKs (ST ''poincare_onto_from_in: strange interval'') (l ≤ u);
_ ← mk_safe {l .. u};
subset_spec1_coll R P dP
}"

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

definition "lvivl'_invar n x =
(case x of None ⇒ True | Some (l, u) ⇒ length l = length u ∧ length u = n)"

definition "one_step_until_time_ivl_in_ivl X0 (t1::real) (t2::real) R dR =
do {
(X, CX) ← one_step_until_time_ivl X0 True t1 t2;
((l, u), X) ← scaleR2_rep X;
CHECKs (ST ''one_step_until_time_ivl_in_ivl: there should not be scaleR2'') (l = 1 ∧ u = 1);
(l, u) ← ivl_rep R;
CHECKs (ST ''one_step_until_time_ivl_in_ivl: strange interval'') (l ≤ u);
let _ = trace_set1 (ST ''final step to:'') (Some X);
let _ = trace_set (ST ''contained in?'') (Some {l .. u});
let _ = print_set1 False X;
let _ = print_set False {l .. u};
subset_spec1 X R dR
}"

definition "poincare_onto_in_ivl
(guards)                    ― ‹avoiding guards›
(ivl::'n rvec set)          ― ‹onto ‹ivl››
sctn                        ― ‹which is part of ‹sctn››
ro
(XS0::'n::enum eucl1 set)
P dP =
do {
RS ← poincare_onto_series empty_symstart empty_trap guards XS0 ivl sctn ro;
((l, u), R) ← scaleR2_rep_coll RS;
CHECKs (ST ''poincare_onto_in_ivl: there should not be scaleR2'') (l = 1 ∧ u = 1);
(l, u) ← ivl_rep P;
CHECKs (ST ''poincare_onto_in_ivl: strange interval'') (l ≤ u);
(lR, uR) ← ivl_rep_of_set_coll (op_image_fst_coll R);
CHECKs (ST ''poincare_onto_in_ivl: strange interval2'') (lR ≤ uR);
let _ = trace_set (ST ''final step to:'') (Some {lR .. uR});
let _ = trace_set (ST ''contained in?'') (Some {l .. u});
_ ← mk_safe {l .. u};
subset_spec1_coll R P dP
}"

definition "poincare_maps_onto Σ X0 X1 ⟷ poincare_mapsto Σ X0 UNIV (Csafe - Σ) X1"

end

end

end```