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))::'aL'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_relnres_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_relscaleR2_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_relsctn_rel), Ainvar_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_reloption_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 
      (xPS. (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)) 
          ((xPS  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 optionunit" where "trace_set1e _ _ = ()"
definition trace_set1::"string'a set optionunit" 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, sinfo_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_relscaleR2_rel)):::appr1_relscaleR2_rel ×r clw_rel (appr1_relscaleR2_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_relscaleR2_rel);
        RETURN (mk_coll (a:::appr1_relscaleR2_rel)  mk_coll (b:::appr1_relscaleR2_rel)  Xs', Ys)
      }
    }) (X:::clw_rel (appr1_relscaleR2_rel), op_empty_coll:::clw_rel (appr1_relscaleR2_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_boolii_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_relnres_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:::Idphantom_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:::Idphantom_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_relnres_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_relnres_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_relnres_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_relinfo_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_relnres_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_relnres_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_relinfo_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_relnres_rel) guardsi X0:::
      clw_rel appr_rel ×r
       clw_rel (iplane_rel lvivl_rel::(_×'n rvec set) set, iinfo_rel appr1e_relinfo_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_relnres_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_relnres_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_relnres_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_relnres_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_relnres_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_relsctn_rel ×r
        clw_rel (isbelows_rel appr_rel) ×r clw_rel appr_rellist_wset_relnres_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_rellist_wset_relnres_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_relsctn_rel
            ×r clw_rel (isbelows_rel appr_rel) ×r clw_rel appr_rellist_wset_relnres_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_relnres_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_relscaleR2_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 (widthpig)
  }"

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_relnres_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_relsctn_rel ×r
      clw_rel (isbelows_rel appr_rel) ×r clw_rel appr_rellist_wset_relnres_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_relnres_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_relsctn_rel ×r
      clw_rel (isbelows_rel appr_rel) ×r clw_rel appr_rellist_wset_relnres_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 = (xset 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