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 (norm2e 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