Theory Continuous_Time_Markov_Chain

(* Author: Johannes Hölzl <hoelzl@in.tum.de> *)

section ‹Continuous-time Markov chains›

theory Continuous_Time_Markov_Chain
  imports Discrete_Time_Markov_Process Discrete_Time_Markov_Chain
begin

subsection ‹Trace Operations: relate @{typ "('a × real) stream"} and @{typ "real  'a"}

partial_function (tailrec) trace_at :: "'a  (real × 'a) stream  real  'a"
where
  "trace_at s ω j = (case ω of (t', s')##ω  if t'  j then trace_at s' ω j else s)"

lemma trace_at_simp[simp]: "trace_at s ((t', s')##ω) j = (if t'  j then trace_at s' ω j else s)"
  by (subst trace_at.simps) simp

lemma trace_at_eq:
  "trace_at s ω j = (case sfirst (λx. j < fst (shd x)) ω of   undefined | enat i  (s ## smap snd ω) !! i)"
proof (split enat.split; safe)
  assume "sfirst (λx. j < fst (shd x)) ω = "
  with sfirst_finite[of "λx. j < fst (shd x)" ω]
  have "alw (λx. fst (shd x)  j) ω"
    by (simp add: not_ev_iff not_less)
  then show "trace_at s ω j = undefined"
    by (induction arbitrary: s ω rule: trace_at.fixp_induct) (auto split: stream.split)
next
  show "sfirst (λx. j < fst (shd x)) ω = enat n  trace_at s ω j = (s ## smap snd ω) !! n" for n
  proof (induction n arbitrary: s ω)
    case 0 then show ?case
      by (subst trace_at.simps) (auto simp add: enat_0 sfirst_eq_0 split: stream.split)
  next
    case (Suc n) show ?case
      using sfirst.simps[of "λx. j < fst (shd x)" ω] Suc.prems Suc.IH[of "stl ω" "snd (shd ω)"]
      by (cases ω) (auto simp add: eSuc_enat[symmetric] split: stream.split if_split_asm)
  qed
qed

lemma trace_at_shift: "trace_at s (smap (λ(t, s'). (t + t', s')) ω) t = trace_at s ω (t - t')"
  by (induction arbitrary: s ω rule: trace_at.fixp_induct) (auto split: stream.split)

primcorec merge_at :: "(real × 'a) stream  real  (real × 'a) stream  (real × 'a) stream"
where
  "merge_at ω j ω' = (case ω of (t, s) ## ω  if t  j then (t, s)##merge_at ω j ω' else ω')"

lemma merge_at_simp[simp]: "merge_at (x##ω) j ω' = (if fst x  j then x##merge_at ω j ω' else ω')"
  by (cases x) (subst merge_at.code; simp)

subsection ‹Exponential Distribution›

definition exponential :: "real  real measure"
where
  "exponential l = density lborel (exponential_density l)"

lemma space_exponential: "space (exponential l) = UNIV"
  by (simp add: exponential_def)

lemma sets_exponential[measurable_cong]: "sets (exponential l) = sets borel"
  by (simp add: exponential_def)

lemma prob_space_exponential: "0 < l  prob_space (exponential l)"
  unfolding exponential_def by (intro prob_space_exponential_density)

lemma AE_exponential: "0 < l  AE x in exponential l. 0 < x"
  unfolding exponential_def using AE_lborel_singleton[of 0] by (auto simp add: AE_density exponential_density_def)

lemma emeasure_exponential_Ioi_cutoff:
  assumes "0 < l"
  shows "emeasure (exponential l) {x <..} = exp (- (max 0 x) * l)"
proof -
  interpret prob_space "exponential l"
    unfolding exponential_def using 0<l by (rule prob_space_exponential_density)
  have *: "prob {xa  space (exponential l). max 0 x < xa} = exp (- max 0 x * l)"
    apply (rule exponential_distributedD_gt[OF _ _ 0<l])
    apply (auto simp: exponential_def distributed_def)
    apply (subst (6) distr_id[symmetric])
    apply (subst (2) distr_cong)
    apply simp_all
    done
  have "emeasure (exponential l) {x <..} = emeasure (exponential l) {max 0 x <..}"
    using AE_exponential[OF 0<l] by (intro emeasure_eq_AE) auto
  also have " = exp (- (max 0 x) * l)"
    using * unfolding emeasure_eq_measure by (simp add: space_exponential greaterThan_def)
  finally show ?thesis .
qed

lemma emeasure_exponential_Ioi:
  "0 < l  0  x  emeasure (exponential l) {x <..} = exp (- x * l)"
  using emeasure_exponential_Ioi_cutoff[of l x] by simp

lemma exponential_eq_stretch:
  assumes "0 < l"
  shows "exponential l = distr (exponential 1) borel (λx. (1/l) * x)"
proof (intro measure_eqI)
  fix A assume "A  sets (exponential l)"
  then have [measurable]: "A  sets borel"
    by (simp add: sets_exponential)
  then have [measurable]: "(λx. x / l) -` A  sets borel"
    by (rule measurable_sets_borel[rotated]) simp
  have "emeasure (exponential l) A =
    (+x. ennreal l * (indicator (((*) (1/l) -` A)  {0 ..}) (l * x) * ennreal (exp (- (l * x)))) lborel)"
    using 0 < l
    by (auto simp: ac_simps emeasure_distr exponential_def emeasure_density exponential_density_def
                   ennreal_mult zero_le_mult_iff
             intro!: nn_integral_cong split: split_indicator)
  also have " = (+x. indicator (((*) (1/l) -` A)  {0 ..}) x * ennreal (exp (- x)) lborel)"
    using 0<l
    apply (subst nn_integral_stretch)
      apply (auto simp: nn_integral_cmult)
    apply (simp add: ennreal_mult[symmetric] mult.assoc[symmetric])
    done
  also have " = emeasure (distr (exponential 1) borel (λx. (1/l) * x)) A"
    by (auto simp add: emeasure_distr exponential_def emeasure_density exponential_density_def
        intro!: nn_integral_cong split: split_indicator)
  finally show "emeasure (exponential l) A = emeasure (distr (exponential 1) borel (λx. (1/l) * x)) A" .
qed (simp add: sets_exponential)

lemma uniform_measure_exponential:
  assumes "0 < l" "0  t"
  shows "uniform_measure (exponential l) {t <..} = distr (exponential l) borel ((+) t)" (is "?L = ?R")
proof (rule measure_eqI_lessThan)
  fix x
  have "0 < emeasure (exponential l) {t<..}"
    unfolding emeasure_exponential_Ioi[OF assms] by simp
  with assms show "?L {x<..} < "
    by (simp add: ennreal_divide_eq_top_iff less_top[symmetric] lessThan_Int_lessThan
      emeasure_exponential_Ioi)
  have *: "((+) t -` {x<..}  space (exponential l)) = {x - t <..}"
    by (auto simp: space_exponential)
  show "?L {x<..} = ?R {x<..}"
    using assms by (simp add: lessThan_Int_lessThan emeasure_exponential_Ioi divide_ennreal
      emeasure_distr * emeasure_exponential_Ioi_cutoff exp_diff[symmetric] field_simps split: split_max)
qed (auto simp: sets_exponential)

lemma emeasure_PiM_exponential_Ioi_finite:
  assumes "J  I" "finite J" "i. i  I  0 < R i" "0  x"
  shows "emeasure (ΠM iI. exponential (R i)) (prod_emb I (λi. exponential (R i)) J (ΠE jJ. {x<..})) = exp (- x * (iJ. R i))"
proof (subst emeasure_PiM_emb)
  from assms show "(iJ. emeasure (exponential (R i)) {x<..}) = ennreal (exp (- x * sum R J))"
    by (subst prod.cong[OF refl emeasure_exponential_Ioi])
       (auto simp add: prod_ennreal exp_sum sum_negf[symmetric] sum_distrib_left)
qed (insert assms, auto intro!: prob_space_exponential)

lemma emeasure_PiM_exponential_Ioi_sequence:
  assumes R: "summable R" "i. 0 < R i" "0  x"
  shows "emeasure (ΠM iUNIV. exponential (R i)) (Π iUNIV. {x<..}) = exp (- x * suminf R)"
proof -
  let ?R = "λi. exponential (R i)" let ?P = "ΠM iUNIV. ?R i"
  let ?N = "λn::nat. prod_emb UNIV ?R {..<n} (ΠE i{..<n}. {x<..})"
  interpret prob_space ?P
    by (intro prob_space_PiM prob_space_exponential R)
  have "(ΠM iUNIV. exponential (R i)) (n. ?N n) = (INF n. (ΠM iUNIV. exponential (R i)) (?N n))"
    by (intro INF_emeasure_decseq[symmetric] decseq_emb_PiE) (auto simp: incseq_def)
  also have " = (INF n. ennreal (exp (- x * (i<n. R i))))"
    using R by (intro INF_cong emeasure_PiM_exponential_Ioi_finite) auto
  also have " = ennreal (exp (- x * (SUP n. (i<n. R i))))"
    using R
    by (subst continuous_at_Sup_antimono[where f="λr. ennreal (exp (- x * r))"])
       (auto intro!: bdd_aboveI2[where M="i. R i"] sum_le_suminf summable_mult mult_left_mono
                     continuous_mult continuous_at_ennreal continuous_within_exp[THEN continuous_within_compose3] continuous_minus
             simp: less_imp_le antimono_def image_comp)
  also have " = ennreal (exp (- x * (i. R i)))"
    using R by (subst suminf_eq_SUP_real) (auto simp: less_imp_le)
  also have "(n. ?N n) = (Π iUNIV. {x<..})"
    by (fastforce simp: prod_emb_def Pi_iff PiE_iff space_exponential)
  finally show ?thesis
    using R by simp
qed

lemma emeasure_PiM_exponential_Ioi_countable:
  assumes R: "J  I" "countable J" "i. i  I  0 < R i" "0  x" and finite: "integrable (count_space J) R"
  shows "emeasure (ΠM iI. exponential (R i)) (prod_emb I (λi. exponential (R i)) J (ΠE jJ. {x<..})) =
    exp (- x * (LINT i|count_space J. R i))"
proof cases
  assume "finite J" with assms show ?thesis
    by (subst emeasure_PiM_exponential_Ioi_finite)
       (auto simp: lebesgue_integral_count_space_finite)
next
  assume "infinite J"
  let ?R = "λi. exponential (R i)" let ?P = "ΠM iI. ?R i"
  define f where "f = from_nat_into J"
  have J_eq: "J = range f" and f: "inj f" "f  UNIV  I"
    using from_nat_into_inj_infinite[of J] range_from_nat_into[of J] countable J infinite J J  I
    by (auto simp: inj_on_def f_def simp del: range_from_nat_into)
  have Bf: "bij_betw f UNIV J"
    unfolding J_eq using inj_on_imp_bij_betw[OF f(1)] .

  have summable_R: "summable (λi. R (f i))"
    using finite unfolding integrable_bij_count_space[OF Bf, symmetric] integrable_count_space_nat_iff
    by (rule summable_norm_cancel)

  have "emeasure (ΠM iUNIV. exponential (R (f i))) (Π iUNIV. {x<..}) = exp (- x * (i. R (f i)))"
    using finite assms unfolding J_eq by (intro emeasure_PiM_exponential_Ioi_sequence[OF summable_R]) auto
  also have "(ΠM iUNIV. exponential (R (f i))) = distr ?P (ΠM iUNIV. exponential (R (f i))) (λω. λiUNIV. ω (f i))"
    using R by (intro distr_PiM_reindex[symmetric, OF _ f] prob_space_exponential) auto
  also have " (Π iUNIV. {x<..}) = ?P ((λω. λiUNIV. ω (f i)) -` (Π iUNIV. {x<..})  space ?P)"
    using f(2) by (intro emeasure_distr infprod_in_sets) (auto simp: Pi_iff)
  also have "(λω. λiUNIV. ω (f i)) -` (Π iUNIV. {x<..})  space ?P = prod_emb I ?R J (ΠE jJ. {x<..})"
    by (auto simp: prod_emb_def space_PiM space_exponential Pi_iff J_eq)
  also have "(i. R (f i)) = (LINT i|count_space J. R i)"
    using finite
    by (subst integral_count_space_nat[symmetric])
       (auto simp: integrable_bij_count_space[OF Bf] integral_bij_count_space[OF Bf])
  finally show ?thesis .
qed

lemma AE_PiM_exponential_suminf_infty:
  fixes R :: "nat  real"
  assumes R: "n. 0 < R n" and finite: "(n. ennreal (1 / R n)) = top"
  shows "AE ω in ΠM nUNIV. exponential (R n). (n. ereal (ω n)) = "
proof -
  let ?P = "ΠM nUNIV. exponential (R n)"
  interpret prob_space "exponential (R n)" for n
    by (intro prob_space_exponential R)
  interpret product_prob_space "λn. exponential (R n)" UNIV
    proof qed

  have AE_pos: "AE ω in ?P. i. 0 < ω i"
    unfolding AE_all_countable by (intro AE_PiM_component allI prob_space_exponential R AE_exponential) simp

  have indep: "indep_vars (λi. borel) (λi x. x i) UNIV"
    using PiM_component
    apply (subst P.indep_vars_iff_distr_eq_PiM)
     apply (auto simp: restrict_UNIV distr_id2)
    apply (subst distr_id2)
     apply (intro sets_PiM_cong)
      apply (auto simp: sets_exponential cong: distr_cong)
    done

  have [simp]: "0  x + x * R i  0  x" for x i
    using zero_le_mult_iff[of x "1 + R i"] R[of i] by (simp add: field_simps)

  have "(+ω. eexp (n. - ereal (ω n)) ?P) = (+ω. (INF n. i<n. eexp (- ereal (ω i))) ?P)"
  proof (intro nn_integral_cong_AE, use AE_pos in eventually_elim)
    fix ω :: "nat  real" assume ω: "i. 0 < ω i"
    show "eexp (n. - ereal (ω n)) = (n. i<n. eexp (- ereal (ω i)))"
    proof (rule LIMSEQ_unique[OF _ LIMSEQ_INF])
      show "(λi. i<i. eexp (- ereal (ω i)))  eexp (n. - ereal (ω n))"
        using ω by (intro eexp_suminf summable_minus_ereal summable_ereal_pos) (auto intro: less_imp_le)
      show "decseq (λn. i<n. eexp (- ereal (ω i)))"
        using ω by (auto simp: decseq_def intro!: prod_mono3 intro: less_imp_le)
    qed
  qed
  also have " = (INF n. (+ω. (i<n. eexp (- ereal (ω i))) ?P))"
  proof (intro nn_integral_monotone_convergence_INF_AE')
    show "AE ω in ?P. (i<Suc n. eexp (- ereal (ω i)))  (i<n. eexp (- ereal (ω i)))" for n
      using AE_pos
    proof eventually_elim
      case (elim ω)
      show ?case
        by (rule prod_mono3) (auto simp: elim le_less)
    qed
  qed (auto simp: less_top[symmetric])
  also have " = (INF n. (i<n. (+ω. eexp (- ereal (ω i)) ?P)))"
  proof (intro INF_cong refl indep_vars_nn_integral)
    show "indep_vars (λ_. borel) (λi ω. eexp (- ereal (ω i))) {..<n}" for n
    proof (rule indep_vars_compose2[of _ _ _ "λi x. eexp(- ereal x)"])
      show "indep_vars (λi. borel) (λi x. x i) {..<n}"
        by (rule indep_vars_subset[OF indep]) auto
    qed auto
  qed auto
  also have " = (INF n. (i<n. R i * (+x. indicator {0 ..} ((1 + R i) * x) * ennreal (exp (- ((1 + R i) * x))) lborel)))"
    by (subst product_nn_integral_component)
       (auto simp: field_simps exponential_def nn_integral_density ennreal_mult'[symmetric] ennreal_mult''[symmetric]
                   exponential_density_def exp_diff exp_minus nn_integral_cmult[symmetric]
             intro!: INF_cong prod.cong nn_integral_cong split: split_indicator)
  also have " = (INF n. (i<n. ennreal (R i / (1 + R i))))"
  proof (intro INF_cong prod.cong refl)
    show "R i * (+ x. indicator {0..} ((1 + R i) * x) * ennreal (exp (- ((1 + R i) * x))) lborel) =
      ennreal (R i / (1 + R i))" for i
      using nn_intergal_power_times_exp_Ici[of 0] 0 < R i
      by (subst nn_integral_stretch[where c="1 + R i"])
         (auto simp: mult.assoc[symmetric] ennreal_mult''[symmetric] less_imp_le mult.commute)
  qed
  also have " = (INF n. ennreal (i<n. R i / (1 + R i)))"
    using R by (intro INF_cong refl prod_ennreal divide_nonneg_nonneg) (auto simp: less_imp_le)
  also have " = (INF n. ennreal (inverse (i<n. (1 + R i) / R i)))"
    by (subst prod_inversef[symmetric]) simp_all
  also have " = (INF n. inverse (ennreal (i<n. (1 + R i) / R i)))"
    using R by (subst inverse_ennreal) (auto intro!: prod_pos divide_pos_pos simp: add_pos_pos)
  also have " = inverse (SUP n. ennreal (i<n. (1 + R i) / R i))"
    by (subst continuous_at_Sup_antimono [where f = inverse])
      (auto simp: antimono_def image_comp intro!: continuous_on_imp_continuous_within[OF continuous_on_inverse_ennreal'])
  also have "(SUP n. ennreal (i<n. (1 + R i) / R i)) = top"
  proof (cases "SUP n. ennreal (i<n. (1 + R i) / R i)")
    case (real r)
    have "(λn. ennreal (i<n. (1 + R i) / R i))  r"
      using R unfolding real(2)[symmetric]
      by (intro LIMSEQ_SUP monoI ennreal_leI prod_mono2) (auto intro!: divide_nonneg_nonneg add_nonneg_nonneg intro: less_imp_le)
    then have "(λn. (i<n. (1 + R i) / R i))  r"
      by (rule tendsto_ennrealD)
         (use R real in auto intro!: always_eventually prod_nonneg divide_nonneg_nonneg add_nonneg_nonneg intro: less_imp_le)
    moreover have "(1 + R i) / R i = 1 / R i + 1" for i
      using 0 < R i by (auto simp: field_simps)
    ultimately have "convergent (λn. i<n. 1 / R i + 1)"
      by (auto simp: convergent_def)
    then have "summable (λi. 1 / R i)"
      using R by (subst summable_iff_convergent_prod) (auto intro: less_imp_le)
    moreover have "0  1 / R i" for i
      using R by (auto simp: less_imp_le)
    ultimately show ?thesis
      using finite ennreal_suminf_neq_top[of "λi. 1 / R i"] by blast
  qed
  finally have "(+ω. eexp (n. - ereal (ω n)) ?P) = 0"
    by simp
  then have "AE ω in ?P. eexp (n. - ereal (ω n)) = 0"
    by (subst (asm) nn_integral_0_iff_AE) auto
  then show ?thesis
    using AE_pos
  proof eventually_elim
    show "(i. 0 < ω i)  eexp (n. - ereal (ω n)) = 0  (n. ereal (ω n)) = " for ω
      apply (auto simp del: uminus_ereal.simps simp add: uminus_ereal.simps[symmetric]
                  intro!: summable_iff_suminf_neq_top intro: less_imp_le)
      apply (subst (asm) suminf_minus_ereal)
      apply (auto intro!: summable_ereal_pos intro: less_imp_le)
      done
  qed
qed

subsection ‹Transition Rates›

locale transition_rates =
  fixes R :: "'a  'a  real"
  assumes R_nonneg[simp]: "x y. 0  R x y"
  assumes R_diagonal_0[simp]: "x. R x x = 0"
  assumes finite_weight: "x. (+y. R x y count_space UNIV) < "
  assumes positive_weight: "x. 0 < (+y. R x y count_space UNIV)"
begin

abbreviation S :: "(real × 'a) measure"
where "S  (borel M count_space UNIV)"

abbreviation T :: "(real × 'a) stream measure"
where "T  stream_space S"

abbreviation I :: "'a  'a set"
where "I x  {y. 0 < R x y}"

lemma I_countable: "countable (I x)"
proof -
  let ?P = "point_measure UNIV (R x)"
  interpret finite_measure ?P
  proof
    show "emeasure ?P (space ?P)  "
      using finite_weight
      by (simp add: emeasure_density point_measure_def less_top)
  qed
  from countable_support emeasure_point_measure_finite2[of "{_}" UNIV "R x"]
  show ?thesis
    by (simp add: emeasure_eq_measure less_le)
qed

definition escape_rate :: "'a  real" where
  "escape_rate x = y. R x y count_space UNIV"

lemma ennreal_escape_rate: "ennreal (escape_rate x) = (+y. R x y count_space UNIV)"
  using finite_weight[of x] unfolding escape_rate_def
  by (intro nn_integral_eq_integral[symmetric]) (auto simp: integrable_iff_bounded)

lemma escape_rate_pos: "0 < escape_rate x"
  using positive_weight unfolding ennreal_escape_rate[symmetric] by simp

lemma nonneg_escape_rate[simp]: "0  escape_rate x"
  using escape_rate_pos[THEN less_imp_le] .

lemma prob_space_exponential_escape_rate: "prob_space (exponential (escape_rate x))"
  using escape_rate_pos by (rule prob_space_exponential)

lemma measurable_escape_rate[measurable]: "escape_rate  count_space UNIV M borel"
  by auto

lemma measurable_exponential_escape_rate[measurable]: "(λx. exponential (escape_rate x))  count_space UNIV M prob_algebra borel"
  by (auto simp: space_prob_algebra sets_exponential prob_space_exponential_escape_rate)

interpretation pmf_as_function .

lift_definition J :: "'a  'a pmf" is "λx y.  R x y / escape_rate x"
proof safe
  show "0  R x y / escape_rate x" for x y
    by (auto intro!: integral_nonneg_AE divide_nonneg_nonneg R_nonneg simp: escape_rate_def)
  show "(+y. R x y / escape_rate x count_space UNIV) = 1" for x
    using escape_rate_pos[of x]
    by (auto simp add: divide_ennreal[symmetric] nn_integral_divide ennreal_escape_rate[symmetric] intro!: ennreal_divide_self)
qed

lemma set_pmf_J: "set_pmf (J x) = I x"
  using escape_rate_pos[of x] by (auto simp: set_pmf_iff J.rep_eq less_le)

interpretation exp_esc: pair_prob_space "distr (exponential (escape_rate x)) borel ((+) t)" "J x" for x
proof -
  interpret prob_space "distr (exponential (escape_rate x)) borel ((+) t)"
    by (intro prob_space.prob_space_distr prob_space_exponential_escape_rate) simp
  show "pair_prob_space (distr (exponential (escape_rate x)) borel ((+) t)) (measure_pmf (J x))"
    by standard
qed

subsection ‹Continuous-time Kernel›

definition K :: "(real × 'a)  (real × 'a) measure" where
  "K = (λ(t, x). (distr (exponential (escape_rate x)) borel ((+) t)) M J x)"

interpretation K: discrete_Markov_process "borel M count_space UNIV" K
proof
  show "K  borel M count_space UNIV M prob_algebra (borel M count_space UNIV)"
    unfolding K_def
    apply measurable
    apply (rule measurable_snd[THEN measurable_compose])
    apply (auto simp: space_prob_algebra prob_space_measure_pmf)
    done
qed

interpretation DTMC: MC_syntax J .

lemma in_space_S[simp]: "x  space S"
  by (simp add: space_pair_measure)

lemma in_space_T[simp]: "x  space T"
  by (simp add: space_pair_measure space_stream_space)

lemma in_space_lim_stream: "ω  space (K.lim_stream x)"
  unfolding K.space_lim_stream space_stream_space[symmetric] by simp

lemma prob_space_K_lim: "prob_space (K.lim_stream x)"
  using K.lim_stream[THEN measurable_space] by (simp add: space_prob_algebra)

definition select_first :: "'a  ('a  real)  'a  bool"
where "select_first x p y = (y  I x  (y'I x - {y}. p y < p y'))"

lemma select_firstD1: "select_first x p y  y  I x"
  by (simp add: select_first_def)

lemma select_first_unique:
  assumes y: "select_first x p y1" " select_first x p y2" shows "y1 = y2"
proof -
  have "y1  y2  p y1 < p y2" "y1  y2  p y2 < p y1"
    using y by (auto simp: select_first_def)
  then show "y1 = y2"
    by (rule_tac ccontr) auto
qed

lemma The_select_first[simp]: "select_first x p y  The (select_first x p) = y"
  by (intro the_equality select_first_unique)

lemma select_first_INF:
  "select_first x p y  (INF xI x. p x) = p y"
  by (intro antisym cINF_greatest cINF_lower bdd_belowI2[where m="p y"])
     (auto simp: select_first_def le_less)

lemma measurable_select_first[measurable]:
  "(λp. select_first x p y)  (ΠM yI x. borel) M count_space UNIV"
  using I_countable unfolding select_first_def by (intro measurable_pred_countable pred_intros_conj1') measurable

lemma measurable_THE_select_first[measurable]:
  "(λp. The (select_first x p))  (ΠM yI x. borel) M count_space UNIV"
  by (rule measurable_THE) (auto intro: select_first_unique I_countable dest: select_firstD1)

lemma sets_S_eq: "sets S = sigma_sets UNIV { {t ..} × A | t A. A  - I x  (sI x. A = {s}) }"
proof (subst sets_pair_eq)
  let ?CI = "λa::real. {a ..}" let ?Ea = "range ?CI"

  show "?Ea  Pow (space borel)" "sets borel = sigma_sets (space borel) ?Ea"
    unfolding borel_Ici by auto
  show "?CI`Rats  ?Ea" "(iRats. ?CI i) = space borel"
    using Rats_dense_in_real[of "x - 1" "x" for x] by (auto intro: less_imp_le)

  let ?Eb = "Pow (- I x)  (λs. {s}) ` I x"
  have "b  sigma_sets UNIV (Pow (- I x)  (λs. {s}) ` I x)" for b
  proof -
    have "b = (b - I x)  (xb  I x. {x})"
      by auto
    also have "  sigma UNIV (Pow (- I x)  (λs. {s}) ` I x)"
      using I_countable by (intro sets.Un sets.countable_UN') auto
    finally show ?thesis
      by simp
  qed
  then show "sets (count_space UNIV) = sigma_sets (space (count_space UNIV)) ?Eb"
    by auto
  show "countable ({- I x}  (sI x. {{s}}))"
    using I_countable by auto
  show "sets (sigma (space borel × space (count_space UNIV)) {a × b |a b. a  ?Ea  b  ?Eb}) =
    sigma_sets UNIV {{t ..} × A |t A. A  - I x  (sI x. A = {s})}"
    apply simp
    apply (intro arg_cong[where f="sigma_sets _"])
    apply auto
    done
qed (auto intro: countable_rat)

subsection ‹Kernel equals Parallel Choice›

abbreviation PAR :: "'a  ('a  real) measure"
where
  "PAR x  (ΠM yI x. exponential (R x y))"

lemma PAR_least:
  assumes y: "y  I x"
  shows "PAR x {pspace (PAR x). t  p y  select_first x p y} =
     emeasure (exponential (escape_rate x)) {t ..} * ennreal (pmf (J x) y)"
proof -
  let ?E = "λy. exponential (R x y)" let ?P' = "ΠM yI x - {y}. ?E y"
  interpret P': prob_space ?P'
    by (intro prob_space_PiM prob_space_exponential) simp
  have *: "PAR x = (ΠM yinsert y (I x - {y}). ?E y)"
    using y by (intro PiM_cong) auto
  have "0 < R x y"
    using y by simp
  have **: "(λ(x, X). X(y := x))  exponential (R x y) M PiM (I x - {y}) (λi. exponential (R x i)) M PAR x"
    using y
    apply (subst measurable_cong_sets[OF sets_pair_measure_cong[OF sets_exponential sets_PiM_cong[OF refl sets_exponential]] sets_PiM_cong[OF refl sets_exponential]])
    apply measurable
    apply (rule measurable_fun_upd[where J="I x - {y}"])
    apply auto
    done
  have "PAR x {pspace (PAR x). t  p y  (y'I x-{y}. p y < p y')} =
    (+ty. indicator {t..} ty * ?P' {pspace ?P'. y'I x-{y}. ty < p y'} ?E y)"
    unfolding * using y  I x
    apply (subst distr_pair_PiM_eq_PiM[symmetric])
    apply (auto intro!: prob_space_exponential simp: emeasure_distr insert_absorb)
    apply (subst emeasure_distr[OF **])
    subgoal
      using I_countable by (auto simp: pred_def[symmetric])
    apply (subst P'.emeasure_pair_measure_alt)
    subgoal
      using I_countable[of x]
      apply (intro measurable_sets[OF **])
      apply (auto simp: pred_def[symmetric])
      done
    apply (auto intro!: nn_integral_cong arg_cong2[where f=emeasure] split: split_indicator if_split_asm
      simp: space_exponential space_PiM space_pair_measure PiE_iff extensional_def)
    done
  also have " = (+ty. indicator {t..} ty * ennreal (exp (- ty * (escape_rate x - R x y))) ?E y)"
    apply (intro nn_integral_cong_AE)
    using AE_exponential[OF 0 < R x y]
  proof eventually_elim
    fix ty :: real assume "0 < ty"
    have "escape_rate x =
      (+y'. R x y' * indicator {y} y' count_space UNIV) + (+y'. R x y' * indicator (I x - {y}) y' count_space UNIV)"
      unfolding ennreal_escape_rate by (subst nn_integral_add[symmetric]) (auto simp: less_le split: split_indicator intro!: nn_integral_cong)
    also have " = R x y + (+y'. R x y' count_space (I x - {y}))"
      by (auto simp add: nn_integral_count_space_indicator less_le simp del: nn_integral_indicator_singleton
               intro!: arg_cong2[where f="(+)"] nn_integral_cong split: split_indicator)
    finally have "(+y'. R x y' count_space (I x - {y})) = escape_rate x - R x y  R x y  escape_rate x"
      using escape_rate_pos[THEN less_imp_le]
      by (cases "(+y'. R x y' count_space (I x - {y}))")
         (auto simp: add_top ennreal_plus[symmetric] simp del: ennreal_plus)
    then have "integrable (count_space (I x - {y})) (R x)" "(LINT y'|count_space (I x - {y}). R x y') = escape_rate x - R x y"
      by (auto simp: nn_integral_eq_integrable)
    then have "?P' (prod_emb (I x-{y}) ?E (I x-{y}) (ΠE j(I x-{y}). {ty<..})) = exp (- ty * (escape_rate x - R x y))"
      using I_countable 0 < ty by (subst emeasure_PiM_exponential_Ioi_countable) auto
    also have "prod_emb (I x-{y}) ?E (I x-{y}) (ΠE j(I x-{y}). {ty<..}) = {pspace ?P'. y'I x-{y}. ty < p y'}"
      by (simp add: set_eq_iff prod_emb_def space_PiM space_exponential ac_simps Pi_iff)
    finally show "indicator {t..} ty * ?P' {pspace ?P'. y'I x-{y}. ty < p y'} =
      indicator {t..} ty * ennreal (exp (- ty * (escape_rate x - R x y)))"
      by simp
  qed
  also have " = (+ty. ennreal (R x y) * (ennreal (exp (- ty * escape_rate x)) * indicator {max 0 t..} ty) lborel)"
    by (auto simp add: exponential_def exponential_density_def nn_integral_density ennreal_mult[symmetric] exp_add[symmetric] field_simps
                intro!: nn_integral_cong split: split_indicator)
  also have " = (R x y / escape_rate x) * emeasure (exponential (escape_rate x)) {max 0 t..}"
    using escape_rate_pos[of x]
    by (auto simp: exponential_def exponential_density_def emeasure_density nn_integral_cmult[symmetric] ennreal_mult[symmetric]
             split: split_indicator intro!: nn_integral_cong )
  also have " = pmf (J x) y * emeasure (exponential (escape_rate x)) {t..}"
    using AE_exponential[OF escape_rate_pos[of x]]
    by (intro arg_cong2[where f="(*)"] emeasure_eq_AE) (auto simp: J.rep_eq )
  finally show ?thesis
    using assms by (simp add: mult_ac select_first_def)
qed

lemma AE_PAR_least: "AE p in PAR x. yI x. select_first x p y"
proof -
  have D: "disjoint_family_on (λy. {p  space (PAR x). select_first x p y}) (I x)"
    by (auto simp: disjoint_family_on_def dest: select_first_unique)
  have "PAR x {pspace (PAR x). yI x. select_first x p y} =
    PAR x (yI x. {pspace (PAR x). select_first x p y})"
    by (auto intro!: arg_cong2[where f=emeasure])
  also have " = (+y. PAR x {pspace (PAR x). select_first x p y} count_space (I x))"
    using I_countable by (intro emeasure_UN_countable D) auto
  also have " = (+y. PAR x {pspace (PAR x). 0  p y  select_first x p y} count_space (I x))"
  proof (intro nn_integral_cong emeasure_eq_AE, goal_cases)
    case (1 y) with AE_PiM_component[of "I x" "λy. exponential (R x y)" y "(<) 0"] AE_exponential[of "R x y"] show ?case
      by (auto simp: prob_space_exponential)
  qed (insert I_countable, auto)
  also have " = (+y. emeasure (exponential (escape_rate x)) {0 ..} * ennreal (pmf (J x) y) count_space (I x))"
    by (auto simp add: PAR_least intro!: nn_integral_cong)
  also have " = (+y. emeasure (exponential (escape_rate x)) {0 ..} J x)"
    by (auto simp: nn_integral_measure_pmf nn_integral_count_space_indicator ac_simps pmf_eq_0_set_pmf set_pmf_J
             simp del: nn_integral_const intro!: nn_integral_cong split: split_indicator)
  also have " = 1"
    using AE_exponential[of "escape_rate x"]
    by (auto intro!: prob_space.emeasure_eq_1_AE prob_space_exponential simp: escape_rate_pos less_imp_le)
  finally show ?thesis
    using I_countable
    by (subst prob_space.AE_iff_emeasure_eq_1 prob_space_PiM prob_space_exponential)
       (auto intro!: prob_space_PiM prob_space_exponential simp del: Set.bex_simps(6))
qed

lemma K_alt: "K (t, x) = distr (ΠM yI x. exponential (R x y)) S (λp. (t + (INF yI x. p y), The (select_first x p)))" (is "_ = ?R")
proof (rule measure_eqI_generator_eq_countable)
  let ?E = "{ {t ..} × A | (t::real) A. A  - I x  (sI x. A = {s}) }"
  show "Int_stable ?E"
    apply (auto simp: Int_stable_def)
    subgoal for t1 A1 t2 A2
      by (intro exI[of _ "max t1 t2"] exI[of _ "A1  A2"]) auto
    subgoal for t1 t2 y1 y2
      by (intro exI[of _ "max t1 t2"] exI[of _ "{y1}  {y2}"]) auto
    done
  show "sets (K (t, x)) = sigma_sets UNIV ?E"
    unfolding K.sets_K[OF in_space_S] by (subst sets_S_eq) rule
  show "sets ?R = sigma_sets UNIV ?E"
    using sets_S_eq by simp
  show "countable ((λ(t, A). {t ..} × A) ` ( × ({- I x}  (λs. {s}) ` I x)))"
    by (intro countable_image countable_SIGMA countable_rat countable_Un I_countable) auto

   have *: "(+) t -` {t'..}  space (exponential (escape_rate x)) = {t' - t..}" for t'
     by (auto simp: space_exponential)
  { fix X assume "X  ?E"
    then consider
        t' s where "s  I x" "X = {t' ..} × {s}"
      | t' A where "A  - I x" "X = {t' ..} × A"
      by auto
    then show "K (t, x) X = ?R X"
    proof cases
      case 1
      have "AE p in PAR x. (t' - t  p s  select_first x p s) =
              (t'  t + (xI x. p x)  The (select_first x p) = s)"
        using AE_PAR_least by eventually_elim (auto dest: select_first_unique simp: select_first_INF)
      with 1 I_countable show ?thesis
        by (auto simp add: K_def measure_pmf.emeasure_pair_measure_Times emeasure_distr emeasure_pmf_single *
          PAR_least[symmetric] intro!: emeasure_eq_AE)
    next
      case 2
      moreover
      then have "emeasure (measure_pmf (J x)) A = 0"
        by (subst AE_iff_measurable[symmetric, where P="λx. x  A"])
           (auto simp: AE_measure_pmf_iff set_pmf_J subset_eq)
      moreover
      have "PAR x ((λp. (t + (p ` (I x)), The (select_first x p))) -` ({t'..} × A)  space (PAR x)) = 0"
        using A  - I x AE_PAR_least[of x] I_countable
        by (subst AE_iff_measurable[symmetric, where P="λp. (t + (p ` (I x)), The (select_first x p))  {t'..} × A"])
           (auto simp del: all_simps(5) simp add: imp_ex imp_conjL subset_eq)
      ultimately show ?thesis
        using I_countable
        by (simp add: K_def measure_pmf.emeasure_pair_measure_Times emeasure_distr *)
    qed }

  interpret prob_space "K ts" for ts
    by (rule K.prob_space_K) simp
  show "emeasure (K (t, x)) a  " for a
    using emeasure_finite by simp
qed (insert  Rats_dense_in_real[of "x - 1" x for x], auto, blast intro: less_imp_le)

lemma AE_K: "AE y in K x. fst x < fst y  snd y  J (snd x)"
  unfolding K_def split_beta
  apply (subst exp_esc.AE_pair_iff[symmetric])
  apply measurable
  apply (simp_all add: AE_distr_iff AE_measure_pmf_iff exponential_def AE_density exponential_density_def cong del: AE_cong)
  using AE_lborel_singleton[of 0]
  apply eventually_elim
  apply simp
  done

lemma AE_lim_stream:
  "AE ω in K.lim_stream x. i. snd ((x ## ω) !! i)  DTMC.acc``{snd x}  snd (ω !! i)  J (snd ((x ## ω) !! i))  fst ((x ## ω) !! i) < fst (ω !! i)"
  (is "AE ω in K.lim_stream x. i. ?P ω i")
  unfolding AE_all_countable
proof
  let ?F = "λi x ω. fst ((x ## ω) !! i)" and ?S = "λi x ω. snd ((x ## ω) !! i)"
  fix i show "AE ω in K.lim_stream x. ?P ω i"
  proof (induction i arbitrary: x)
    case 0 with AE_K[of x] show ?case
      by (subst K.AE_lim_stream) (auto simp add: space_pair_measure cong del: AE_cong)
  next
    case (Suc i)
    show ?case
    proof (subst K.AE_lim_stream, goal_cases)
      case 2 show ?case
        using DTMC.countable_reachable
        by (intro measurable_compose_countable_restrict[where f="?S (Suc i) x"])
           (simp_all del: Image_singleton_iff)
    next
      case 3 show ?case
        apply (simp del: AE_conj_iff cong del: AE_cong)
        using AE_K[of x]
        apply eventually_elim
        subgoal premises K_prems for y
          using Suc
          by eventually_elim (insert K_prems, auto intro: converse_rtrancl_into_rtrancl)
        done
    qed (simp add: space_pair_measure)
  qed
qed

lemma measurable_merge_at[measurable]: "(λ(ω, ω'). merge_at ω j ω')  (T M T) M T"
proof (rule measurable_stream_space2)
  define F where "F x n = (case x of (ω::(real × 'a) stream, ω')  merge_at ω j ω') !! n" for x n
  fix n
  have "(λx. F x n)  stream_space S M stream_space S M S"
  proof (induction n)
    case 0 then show ?case
      by (simp add: F_def split_beta' stream.case_eq_if)
  next
    case (Suc n)
    from Suc[measurable]
    have eq: "F x (Suc n) = (case fst x of (t, s) ## ω  if t  j then F (ω, snd x) n else snd x !! Suc n)" for x
      by (auto simp: F_def split: prod.split stream.split)
    show ?case
      unfolding eq stream.case_eq_if by measurable
  qed
  then show "(λx. (case x of (ω, ω')  merge_at ω j ω') !! n)  stream_space S M stream_space S M S"
    unfolding F_def by auto
qed

lemma measurable_trace_at[measurable]: "(λ(s, ω). trace_at s ω j)  (count_space UNIV M T) M count_space UNIV"
  unfolding trace_at_eq by measurable

lemma measurable_trace_at': "(λ((s, j), ω). trace_at s ω j)  ((count_space UNIV M borel) M T) M count_space UNIV"
  unfolding trace_at_eq split_beta' by measurable

lemma K_time_split:
  assumes "t  j" and [measurable]: "f  S M borel"
  shows "(+x. f x * indicator {j <..} (fst x) K (t, s)) = (+x. f x K (j, s)) * exponential (escape_rate s) {j - t <..}"
proof -
  have "(+ y. + x. f (t + x, y) * indicator {j<..} (t + x) exponential (escape_rate s) J s) =
    (+ y. + x. f (t + x, y) * indicator {j - t<..} x exponential (escape_rate s) J s)"
    by (intro nn_integral_cong) (auto split: split_indicator)
  also have " = (+ y. + x. f (t + x, y) uniform_measure (exponential (escape_rate s)) {j-t <..} J s) *
      emeasure (exponential (escape_rate s)) {j - t <..}"
    using t  j escape_rate_pos
    by (subst nn_integral_uniform_measure)
       (auto simp: nn_integral_divide ennreal_divide_times emeasure_exponential_Ioi)
  also have " = (+ y. + x. f (j + x, y) exponential (escape_rate s) J s) *
      emeasure (exponential (escape_rate s)) {j - t <..}"
    using t  j escape_rate_pos by (simp add: uniform_measure_exponential nn_integral_distr)
  finally show ?thesis
    by (simp add: K_def exp_esc.nn_integral_snd[symmetric] nn_integral_distr)
qed

lemma K_in_space[simp]: "K x  space (prob_algebra S)"
  by (rule measurable_space [OF K.K]) simp

lemma L_in_space[simp]: "K.lim_stream x  space (prob_algebra T)"
  by (rule measurable_space [OF K.lim_stream]) simp

subsection ‹Markov Chain Property›

lemma lim_time_split:
  "t  j  K.lim_stream (t, s) = do { ω  K.lim_stream (t, s) ; ω'  K.lim_stream (j, trace_at s ω j) ; return T (merge_at ω j ω')}"
    (is "_  _ = ?DO t s")
proof (coinduction arbitrary: t s rule: K.lim_stream_eq_coinduct)
  case step let ?L = K.lim_stream

  note measurable_compose[OF measurable_prob_algebraD measurable_emeasure_subprob_algebra, measurable (raw)]

  define B' where "B' = (λ(t', s). if t'  j then ?DO t' s else ?L (t', s))"
  show ?case
  proof (intro bexI conjI AE_I2)
    show [measurable]: "B'  S M prob_algebra T"
      unfolding B'_def by measurable
    show "(t s. y = (t, s)  B' y = ?DO t s  t  j)  ?L y = B' y" for y
      by (cases y; cases "fst y  j") (auto simp: B'_def)
    let ?C = "λx. do { ω  ?L x; ω'  ?L (j, trace_at s (x##ω) j); return T (merge_at (x##ω) j ω') }"
    have "?DO t s = do { x  K (t, s); ?C x }"
      apply (subst K.lim_stream_eq[OF in_space_S])
      apply (subst bind_assoc[OF measurable_prob_algebraD measurable_prob_algebraD])
      apply (subst measurable_cong_sets[OF K.sets_K[OF in_space_S] refl])
      apply measurable
      apply (subst bind_assoc[OF measurable_prob_algebraD measurable_prob_algebraD])
      apply measurable
      apply (subst bind_cong[OF refl bind_cong[OF refl bind_return[OF measurable_prob_algebraD]]])
      apply measurable
      done
    also have " = K (t, s)  (λy. B' y  (λω. return T (y ## ω)))" (is "?DO' = ?R")
    proof (rule measure_eqI)
      have "sets ?DO' = sets T"
        by (intro sets_bind'[OF K_in_space]) measurable
      moreover have "sets ?R = sets T"
        by (intro sets_bind'[OF K_in_space]) measurable
      ultimately show "sets ?DO' = sets ?R"
        by simp
      fix A assume "A  sets ?DO'"
      then have A[measurable]: "A  T"
        unfolding sets ?DO' = sets T .
      have "?DO' A = (+x. ?C x A K (t, s))"
        by (subst emeasure_bind_prob_algebra[OF K_in_space]) measurable
      also have " = (+x. ?C x A * indicator {.. j} (fst x) K (t, s)) +
        (+x. ?C x A * indicator {j <..} (fst x) K (t, s))"
        by (subst nn_integral_add[symmetric]) (auto intro!: nn_integral_cong split: split_indicator)
      also have "(+x. ?C x A * indicator {.. j} (fst x) K (t, s)) =
        (+y. emeasure (B' y  (λω. return T (y ## ω))) A * indicator {.. j} (fst y) K (t, s))"
      proof (intro nn_integral_cong ennreal_mult_right_cong refl arg_cong2[where f=emeasure])
        fix x :: "real × 'a" assume "indicator {..j} (fst x)  (0::ennreal)"
        then have "fst x  j"
          by (auto split: split_indicator_asm)
        then show "?C x = (B' x  (λω. return T (x ## ω)))"
          apply (cases x)
          apply (simp add: B'_def)
          apply (subst bind_assoc[OF measurable_prob_algebraD measurable_prob_algebraD])
          apply measurable
          apply (subst bind_assoc[OF measurable_prob_algebraD measurable_prob_algebraD])
          apply measurable
          apply (subst bind_return)
          apply measurable
          done
      qed
      also have "(+x. ?C x A * indicator {j <..} (fst x) K (t, s)) =
        (+y. emeasure (B' y  (λω. return T (y ## ω))) A * indicator {j <..} (fst y) K (t, s))"
      proof -
        have *: "(+) t -` {j<..} = {j - t <..}"
          by auto

        have "(+x. ?C x A * indicator {j <..} (fst x) K (t, s)) =
          (+x. ?L (j, s) A * indicator {j <..} (fst x) K (t, s))"
          by (intro nn_integral_cong ennreal_mult_right_cong refl arg_cong2[where f=emeasure])
             (auto simp: K.sets_lim_stream bind_return'' bind_const' prob_space_K_lim prob_space_imp_subprob_space split: split_indicator_asm)
        also have " = ?L (j, s) A * exponential (escape_rate s) {j - t <..}"
          by (subst nn_integral_cmult) (simp_all add: K_def exp_esc.nn_integral_snd[symmetric] emeasure_distr space_exponential *)
        also have " = (+x. emeasure (?L x  (λω. return T (x ## ω))) A K (j, s)) * exponential (escape_rate s) {j - t <..}"
          by (subst K.lim_stream_eq) (auto simp: emeasure_bind_prob_algebra[OF K_in_space _ A])
        also have " = (+y. emeasure (?L y  (λω. return T (y ## ω))) A * indicator {j <..} (fst y) K (t, s))"
          using t  j by (rule K_time_split[symmetric]) measurable
        also have " = (+y. emeasure (B' y  (λω. return T (y ## ω))) A * indicator {j <..} (fst y) K (t, s))"
          by (intro nn_integral_cong ennreal_mult_right_cong refl arg_cong2[where f=emeasure])
             (auto simp add: B'_def split: split_indicator_asm)
        finally show ?thesis .
      qed
      also have "(+y. emeasure (B' y  (λω. return T (y ## ω))) A * indicator {.. j} (fst y) K (t, s)) +
        (+y. emeasure (B' y  (λω. return T (y ## ω))) A * indicator {j <..} (fst y) K (t, s)) =
        (+y. emeasure (B' y  (λω. return T (y ## ω))) A K (t, s))"
        by (subst nn_integral_add[symmetric]) (auto intro!: nn_integral_cong split: split_indicator)
      also have " = emeasure (K (t, s)  (λy. B' y  (λω. return T (y ## ω)))) A"
        by (rule emeasure_bind_prob_algebra[symmetric, OF K_in_space _ A]) auto
      finally show "?DO' A = emeasure (K (t, s)  (λy. B' y  (λω. return T (y ## ω)))) A" .
    qed
    finally show "?DO t s = K (t, s)  (λy. B' y  (λω. return T (y ## ω)))" .
  qed
qed (simp add: space_pair_measure)

lemma K_eq: "K (t, s) = distr (exponential (escape_rate s) M J s) S (λ(t', s). (t + t', s))"
proof -
  have "distr (exponential (escape_rate s)) borel ((+) t) M distr (J s) (J s) (λx. x) =
    distr (exponential (escape_rate s) M J s) (borel M J s) (λ(x, y). (t + x, y))"
  proof (intro pair_measure_distr)
    interpret prob_space "distr (measure_pmf (J s)) (measure_pmf (J s)) (λx. x)"
      by (intro measure_pmf.prob_space_distr) simp
    show "sigma_finite_measure (distr (measure_pmf (J s)) (measure_pmf (J s)) (λx. x))"
      by unfold_locales
  qed auto
  also have " = distr (exponential (escape_rate s) M J s) S (λ(x, y). (t + x, y))"
    by (intro distr_cong refl sets_pair_measure_cong) simp
  finally show ?thesis
    by (simp add: K_def)
qed

lemma K_shift: "K (t + t', s) = distr (K (t, s)) S (λ(t, s). (t + t', s))"
  unfolding K_eq by (subst distr_distr) (auto simp: comp_def split_beta' ac_simps)

lemma K_not_empty: "space (K x)  {}"
  by (simp add: K_def space_pair_measure split: prod.split)

lemma lim_stream_not_empty: "space (K.lim_stream x)  {}"
  by (simp add: K.space_lim_stream space_pair_measure split: prod.split)

lemma lim_shift: ― ‹Generalize to bijective function on @{const K.lim_stream} invariant on @{const K}
  "K.lim_stream (t + t', s) = distr (K.lim_stream (t, s)) T (smap (λ(t, s). (t + t', s)))"
  (is "_ = ?D t s")
proof (coinduction arbitrary: t s rule: K.lim_stream_eq_coinduct)
  case step then show ?case
  proof (intro bexI[of _ "λ(t, s). ?D (t - t') s"] conjI)
    show "?D t s = K (t + t', s)  (λy. (case y of (t, s)  ?D (t - t') s)  (λω. return T (y ## ω)))"
      apply (subst K.lim_stream_eq[OF in_space_S])
      apply (subst K_shift)
      apply (subst distr_bind[OF measurable_prob_algebraD K_not_empty])
      apply (measurable; fail)
      apply (measurable; fail)
      apply (subst bind_distr[OF _ measurable_prob_algebraD K_not_empty])
      apply (measurable; fail)
      apply (measurable; fail)
      apply (intro bind_cong refl)
      apply (subst distr_bind[OF measurable_prob_algebraD lim_stream_not_empty])
      apply (measurable; fail)
      apply (measurable; fail)
      apply (simp add: distr_return split_beta)
      apply (subst bind_distr[OF _ measurable_prob_algebraD lim_stream_not_empty])
      apply (measurable; fail)
      apply (measurable; fail)
      apply (simp add: split_beta')
      done
  qed (auto cong: conj_cong intro!: exI[of _ "_ - t'"])
qed simp

lemma lim_0: "K.lim_stream (t, s) = distr (K.lim_stream (0, s)) T (smap (λ(t', s). (t' + t, s)))"
  using lim_shift[of 0 t s] by simp

subsection ‹Explosion time›

definition explosion :: "(real × 'a) stream  ereal"
  where "explosion ω = (SUP i. ereal (fst (ω !! i)))"

lemma ball_less_Suc_eq: "(i<Suc n. P i)  (P 0  (i<n. P (Suc i)))"
  using less_Suc_eq_0_disj by auto

lemma lim_stream_timediff_eq_exponential_1:
  "distr (K.lim_stream ts) (PiM UNIV (λ_. borel))
    (λω i. escape_rate (snd ((ts##ω) !! i)) * (fst (ω !! i) - fst ((ts##ω) !! i))) =
    PiM UNIV (λ_. exponential 1)"
  (is "?D = ?P")
proof (rule measure_eqI_PiM_sequence)
  show "sets ?D = sets (PiM UNIV (λ_. borel))" "sets ?P = sets (PiM UNIV (λ_. borel))"
    by (auto intro!: sets_PiM_cong simp: sets_exponential)
  have [measurable]: "ts  space S"
    by auto
  { interpret prob_space ?D
      by (intro prob_space.prob_space_distr K.prob_space_lim_stream measurable_abs_UNIV) auto
    show "finite_measure ?D"
      by unfold_locales }

  interpret E: prob_space "exponential 1"
    by (rule prob_space_exponential) simp
  interpret P: product_prob_space "λ_. exponential 1" UNIV
    by unfold_locales

  let "distr _ _ (?f ts)" = ?D

  fix A :: "nat  real set" and n :: nat assume A[measurable]: "i. A i  sets borel"
  define n' where "n' = Suc n"
  have "emeasure ?D (prod_emb UNIV (λ_. borel) {..n} (PiE {..n} A)) =
    emeasure (K.lim_stream ts) {ωspace (stream_space S). i<n'. ?f ts ω i  A i}"
    apply (subst emeasure_distr)
      apply (auto intro!: measurable_abs_UNIV arg_cong[where f="emeasure _"])
      apply (auto simp: prod_emb_def K.space_lim_stream space_pair_measure n'_def)
    done
  also have " = (i<n'. emeasure (exponential 1) (A i))"
    using A
  proof (induction n' arbitrary: A ts)
    case 0 then show ?case
      using prob_space.emeasure_space_1[OF prob_space_K_lim]
      by (simp add: K.space_lim_stream space_pair_measure)
  next
    case (Suc n A ts)
    from Suc.prems[measurable]
    have [measurable]: "ts  space S"
      by auto

    have "emeasure (K.lim_stream ts) {ω  space (stream_space S). i<Suc n. ?f ts ω i  A i} =
      (+ts'. indicator (A 0) (escape_rate (snd ts) * (fst ts' - fst ts)) *
        emeasure (K.lim_stream ts') {ω  space (stream_space S). i<n. ?f ts' ω i  A (Suc i)} K ts)"
      apply (subst K.emeasure_lim_stream)
      apply simp
       apply measurable
      apply (auto intro!: nn_integral_cong arg_cong2[where f=emeasure] split: split_indicator
        simp: ball_less_Suc_eq)
      done
    also have " = (+ts'. indicator (A 0) (escape_rate (snd ts) * (fst ts' - fst ts)) K ts) *
      (i<n. emeasure (exponential 1) (A (Suc i)))"
      by (subst Suc.IH) (simp_all add: nn_integral_multc)
    also have "(+ts'. indicator (A 0) (escape_rate (snd ts) * (fst ts' - fst ts)) K ts) =
      (+t. indicator (A 0) (escape_rate (snd ts) * t) exponential (escape_rate (snd ts)))"
      by (simp add: K_def exp_esc.nn_integral_snd[symmetric] nn_integral_distr split: prod.split)
    also have " = emeasure (exponential 1) (A 0)"
      using escape_rate_pos[of "snd ts"]
      by (subst exponential_eq_stretch) (simp_all add: nn_integral_distr)
    also have "emeasure (exponential 1) (A 0) * (i<n. emeasure (exponential 1) (A (Suc i))) =
      (i<Suc n. emeasure (exponential 1) (A i))"
      by (rule prod.lessThan_Suc_shift[symmetric])
    finally show ?case .
  qed
  also have " = emeasure ?P (prod_emb UNIV (λ_. borel) {..<n'} (PiE {..<n'} A))"
    using P.emeasure_PiM_emb[of "{..<n'}" A] by (simp add: prod_emb_def space_exponential)
  finally show "emeasure ?D (prod_emb UNIV (λ_. borel) {..n} (PiE {..n} A)) =
    emeasure ?P (prod_emb UNIV (λ_. borel) {..n} (PiE {..n} A))"
    by (simp add: n'_def lessThan_Suc_atMost)
qed

lemma AE_explosion_infty:
  assumes bdd: "bdd_above (range escape_rate)"
  shows "AE ω in K.lim_stream x. explosion ω = "
proof -
  have "escape_rate undefined  (SUP x. escape_rate x)"
    using bdd by (intro cSUP_upper) auto
  then have SUP_escape_pos: "0 < (SUP x. escape_rate x)"
    using escape_rate_pos[of undefined] by simp
  then have SUP_escape_nonneg: "0  (SUP x. escape_rate x)"
    by (rule less_imp_le)

  have [measurable]: "x  space S" by auto
  have "(i. 1::ennreal) = top"
    by (rule sums_unique[symmetric]) (auto simp: sums_def of_nat_tendsto_top_ennreal)
  then have "AE ω in (PiM UNIV (λ_. exponential 1)). (i. ereal (ω i)) = "
    by (intro AE_PiM_exponential_suminf_infty) auto
  then have "AE ω in K.lim_stream x.
    (i. ereal (escape_rate (snd ((x##ω) !! i)) * (fst (ω !! i) - fst ((x##ω) !! i)))) = "
    apply (subst (asm) lim_stream_timediff_eq_exponential_1[symmetric, of x])
    apply (subst (asm) AE_distr_iff)
    apply (auto intro!: measurable_abs_UNIV)
    done
  then show ?thesis
    using AE_lim_stream
  proof eventually_elim
    case (elim ω)
    then have le: "fst ((x##ω) !! n)  fst ((x ## ω) !! m)" if "n  m" for n m
      by (intro lift_Suc_mono_le[OF _ n  m, of "λi. fst ((x ## ω) !! i)"]) (auto intro: less_imp_le)
    have [simp]: "fst x  fst ((x##ω) !! i)" "fst ((x##ω) !! i)  fst (ω !! i)" for i
      using le[of "i" "Suc i"] le[of 0 i] by auto

    have "(i. ereal (escape_rate (snd ((x ## ω) !! i)) * (fst (ω !! i) - fst ((x ## ω) !! i)))) =
      (SUP n. i<n. ereal (escape_rate (snd ((x ## ω) !! i)) * (fst (ω !! i) - fst ((x ## ω) !! i))))"
      by (intro suminf_ereal_eq_SUP) (auto intro!: mult_nonneg_nonneg)
    also have "  (SUP n. (SUP x. escape_rate x) * (ereal (fst ((x ## ω) !! n)) - ereal (fst x)))"
    proof (intro SUP_least SUP_upper2)
      fix n
      have "(i<n. ereal (escape_rate (snd ((x ## ω) !! i)) * (fst (ω !! i) - fst ((x ## ω) !! i)))) 
        (i<n. ereal ((SUP i. escape_rate i) * (fst (ω !! i) - fst ((x ## ω) !! i))))"
        using elim bdd by (intro sum_mono) (auto intro!: cSUP_upper)
      also have " = (SUP i. escape_rate i) * (i<n. fst ((x ## ω) !! Suc i) - fst ((x ## ω) !! i))"
        using elim bdd by (subst sum_ereal) (auto simp: sum_distrib_left)
      also have " = (SUP i. escape_rate i) * (fst ((x ## ω) !! n) - fst x)"
        by (subst sum_lessThan_telescope) simp
      finally show "(i<n. ereal (escape_rate (snd ((x ## ω) !! i)) * (fst (ω !! i) - fst ((x ## ω) !! i))))
          (SUP x. escape_rate x) * (