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) * (ereal (fst ((x ## ω) !! n)) - ereal (fst x))"
        by simp
    qed simp
    also have " = (SUP x. escape_rate x) * ((SUP n. ereal (fst ((x ## ω) !! n))) - ereal (fst x))"
      using elim SUP_escape_nonneg by (subst SUP_ereal_mult_left) (auto simp: SUP_ereal_minus_left[symmetric])
    also have "(SUP n. ereal (fst ((x ## ω) !! n))) = explosion ω"
      unfolding explosion_def
      apply (intro SUP_eq)
      subgoal for i by (intro bexI[of _ i]) auto
      subgoal for i by (intro bexI[of _ "Suc i"]) auto
      done
    finally show "explosion ω = "
      using elim SUP_escape_pos by (cases "explosion ω") (auto split: if_splits)
  qed
qed

subsection ‹Transition probability $p_t$›

context
begin

declare [[inductive_internals = true]]

inductive trace_in :: "'a set  real  'a  (real × 'a) stream  bool" for S t
where
  "t < t'  s  S  trace_in S t s ((t', s')##ω)"
| "t  t'  trace_in S t s' ω  trace_in S t s ((t', s')##ω)"

end

lemma trace_in_simps[simp]:
  "trace_in ss t s (x##ω) = (if t < fst x then s  ss else trace_in ss t (snd x) ω)"
  by (cases x) (subst trace_in.simps; auto)

lemma trace_in_eq_lfp:
  "trace_in ss t = lfp (λF s. λ(t', s')##ω  if t < t' then s  ss else F s' ω)"
  unfolding trace_in_def by (intro arg_cong[where f=lfp] ext) (auto split: stream.splits)

lemma trace_in_shiftD: "trace_in ss t s ω  trace_in ss (t + t') s (smap (λ(t, s'). (t + t', s')) ω)"
  by (induction rule: trace_in.induct) auto

lemma trace_in_shift[simp]: "trace_in ss t s (smap (λ(t, s'). (t + t', s')) ω)  trace_in ss (t - t') s ω"
  using trace_in_shiftD[of ss t s "smap (λ(t, s'). (t + t', s')) ω" "- t'"]
    trace_in_shiftD[of ss "t - t'" s ω t']
  by (auto simp add: stream.map_comp prod.case_eq_if)

lemma measurable_trace_in':
  "Measurable.pred (borel M count_space UNIV M T) (λ(t, s, ω). trace_in ss t s ω)"
    (is "?M (λ(t, s, ω). trace_in ss t s ω)")
proof -
  let ?F = "λF. λ(t, s, (t', s')##ω)  if t < t' then s  ss else F (t, s', ω)"
  have [measurable]: "Measurable.pred (count_space UNIV) (λx. x  ss)"
    by simp
  have "trace_in ss = (λt s ω. lfp ?F (t, s, ω))"
    unfolding trace_in_def
    apply (subst lfp_arg)
    apply (subst lfp_rolling[where g="λF t s ω. F (t, s, ω)"])
    subgoal by (auto simp: mono_def le_fun_def split: stream.splits)
    subgoal by (auto simp: mono_def le_fun_def split: stream.splits)
    subgoal
      by (intro arg_cong[where f=lfp])
         (auto simp: mono_def le_fun_def split_beta' not_less fun_eq_iff split: stream.splits intro!: arg_cong[where f=lfp])
    done
  then have eq: "(λ(t, s, ω). trace_in ss t s ω) = lfp ?F"
    by simp
  have "sup_continuous ?F"
    by (auto simp: sup_continuous_def fun_eq_iff split: stream.splits)
  then show ?thesis
    unfolding eq
  proof (rule measurable_lfp)
    fix F assume "?M F" then show "?M (?F F)"
      by measurable
  qed
qed

lemma measurable_trace_in[measurable (raw)]:
  assumes [measurable]: "f  M M borel" "g  M M count_space UNIV" "h  M M T"
  shows "Measurable.pred M (λx. trace_in ss (f x) (g x) (h x))"
  using measurable_compose[of "λx. (f x, g x, h x)" M, OF _ measurable_trace_in'[of ss]] by simp

definition p :: "'a  'a  real  real"
where "p s s' t = 𝒫(ω in K.lim_stream (0, s). trace_in {s'} t s ω)"

lemma p[measurable]: "(λ(s, t). p s s' t)  (count_space UNIV M borel) M borel"
proof -
  have *: "(SIGMA x:space (count_space UNIV M borel). {ω  streams (space S). trace_in {s'} (snd x) (fst x) ω}) =
    {xspace ((count_space UNIV M borel) M T). trace_in {s'} (snd (fst x)) (fst (fst x)) (snd x)}"
    by (auto simp: space_pair_measure)

  note measurable_trace_at'[measurable]
  show ?thesis
    unfolding p_def[abs_def] split_beta'
    by (rule measure_measurable_prob_algebra2[where N=T])
       (auto simp: K.space_lim_stream * pred_def[symmetric]
                intro!: pred_count_space_const1 measurable_trace_at'[unfolded split_beta'])
qed

lemma p_nonpos: assumes "t  0" shows "p s s' t = of_bool (s = s')"
proof -
  have "AE ω in K.lim_stream (0, s). trace_in {s'} t s ω = (s = s')"
  proof (subst K.AE_lim_stream)
    show "AE y in K (0, s). AE ω in K.lim_stream y. trace_in {s'} t s (y ## ω) = (s = s')"
      using AE_K
    proof eventually_elim
      fix y :: "real × 'a" assume "fst (0, s) < fst y  snd y  set_pmf (J (snd (0, s)))"
      with t0 show "AE ω in K.lim_stream y. trace_in {s'} t s (y ## ω) = (s = s')"
        by (cases y) auto
    qed
  qed auto
  then have "p s s' t = 𝒫(ω in K.lim_stream (0, s). s = s')"
    unfolding p_def by (intro prob_space.prob_eq_AE K.prob_space_lim_stream) auto
  then show ?thesis
    using prob_space.prob_space[OF K.prob_space_lim_stream] by simp
qed

lemma p_0: "p s s' 0 = of_bool (s = s')"
  using p_nonpos[of 0] by simp

lemma in_sets_T[measurable (raw)]: "Measurable.pred T P  {ω. P ω}  sets T"
  unfolding pred_def by simp

lemma distr_id': "sets M = sets N  distr M N (λx. x) = M"
  by (subst distr_cong[of M M N M _ "λx. x"] ) simp_all

lemma p_nonneg[simp]: "0  p s s' t"
  by (simp add: p_def)

lemma p_le_1[simp]: "p s s' t  1"
  unfolding p_def by (intro prob_space.prob_le_1 K.prob_space_lim_stream) simp

lemma p_eq:
  assumes "0  t"
  shows "p s s'' t = (of_bool (s = s'') + (LINT u:{0..t}|lborel. escape_rate s * exp (escape_rate s * u) * (LINT s'|J s. p s' s'' u))) / exp (t * escape_rate s)"
proof -
  have *: "(+) 0 = (λx::real. x)"
    by auto
  interpret L: prob_space "K.lim_stream x" for x
    by (rule K.prob_space_lim_stream) simp
  interpret E: prob_space "exponential (escape_rate s)" for s
    by (intro escape_rate_pos prob_space_exponential)
  have "p s s'' t = emeasure (K.lim_stream (0, s)) {ωspace T. trace_in {s''} t s ω}"
    by (simp add: p_def L.emeasure_eq_measure K.space_lim_stream space_stream_space del: in_space_T)
  also have " = (+y. emeasure (K.lim_stream y) {ωspace T. trace_in {s''} t s (y##ω) } K (0, s))"
    apply (subst K.lim_stream_eq[OF in_space_S])
    apply (subst emeasure_bind_prob_algebra[OF K_in_space])
    apply (measurable; fail)
    apply (measurable; fail)
    apply (subst bind_return_distr'[OF lim_stream_not_empty])
    apply (measurable; fail)
    apply (simp add: emeasure_distr)
    done
  also have " = (+y. indicator {t <..} (fst y) * of_bool (s = s'') + indicator {0<..t} (fst y) * p (snd y) s'' (t - fst y) K (0, s))"
    apply (intro nn_integral_cong_AE)
    using AE_K
    apply eventually_elim
    subgoal for y
      using L.emeasure_space_1
      apply (cases y)
      apply (auto split: split_indicator simp del: in_space_T)
      subgoal for t' s2
        unfolding p_def L.emeasure_eq_measure[symmetric] K.space_lim_stream space_stream_space[symmetric]
        by (subst lim_0) (simp add: emeasure_distr)
      subgoal
        by (auto split: split_indicator cong: rev_conj_cong simp add: K.space_lim_stream space_stream_space simp del: in_space_T)
      done
    done
  also have " = (+u. +s'. indicator {t <..} u * of_bool (s = s'') +
    indicator {0<..t} u * p s' s'' (t - u) J s exponential (escape_rate s))"
    unfolding K_def
    by (simp add: K_def measure_pmf.nn_integral_fst[symmetric] * distr_id' sets_exponential)
  also have " = ennreal (exp (- t * escape_rate s) * of_bool (s = s'')) +
      (+u. indicator {0<..t} u * +s'. p s' s'' (t - u) J s exponential (escape_rate s))"
    using 0t by (simp add: nn_integral_add nn_integral_cmult ennreal_indicator ennreal_mult emeasure_exponential_Ioi escape_rate_pos)
  also have "(+u. indicator {0<..t} u * +s'. p s' s'' (t - u) J s exponential (escape_rate s)) =
      (+u. indicator {0<..t} u *R (LINT s'|J s. p s' s'' (t - u)) exponential (escape_rate s))"
    by (simp add: measure_pmf.integrable_const_bound[of _ 1] nn_integral_eq_integral ennreal_mult ennreal_indicator)
  also have " = (LINT u:{0<..t}|exponential (escape_rate s). (LINT s'|J s. p s' s'' (t - u)))"
    unfolding set_lebesgue_integral_def
    by (intro nn_integral_eq_integral E.integrable_const_bound[of _ 1] AE_I2)
       (auto intro!: mult_le_one measure_pmf.integral_le_const measure_pmf.integrable_const_bound[of _ 1])
  also have " = (LINT u:{0<..t}|lborel. escape_rate s * exp (- escape_rate s * u) * (LINT s'|J s. p s' s'' (t - u)))"
    unfolding exponential_def set_lebesgue_integral_def
    by (subst integral_density)
       (auto simp: ac_simps exponential_density_def fun_eq_iff split: split_indicator
             simp del: integral_mult_right integral_mult_right_zero intro!: arg_cong2[where f="integralL"])
  also have " = (LINT u:{0..t}|lborel. escape_rate s * exp (- escape_rate s * (t - u)) * (LINT s'|J s. p s' s'' u))"
    using AE_lborel_singleton[of 0] AE_lborel_singleton[of t] unfolding set_lebesgue_integral_def
    by (subst lborel_integral_real_affine[where t=t and c="-1"])
       (auto intro!: integral_cong_AE split: split_indicator)
  also have " = exp (- t * escape_rate s) * escape_rate s * (LINT u:{0..t}|lborel. exp (escape_rate s * u) * (LINT s'|J s. p s' s'' u))"
    by (simp add: field_simps exp_diff exp_minus)
  finally show "p s s'' t = (of_bool (s = s'') + (LBINT u:{0..t}. escape_rate s * exp (escape_rate s * u) * (LINT s'|J s. p s' s'' u))) / exp (t * escape_rate s)"
    unfolding set_lebesgue_integral_def
    by (simp del: ennreal_plus add: ennreal_plus[symmetric] exp_minus field_simps)
qed

lemma continuous_on_p: "continuous_on A (p s s')"
proof -
  interpret E: prob_space "exponential (escape_rate s'')" for s''
    by (intro escape_rate_pos prob_space_exponential)
  have "continuous_on {..0} (p s s')"
    by (simp add: p_nonpos continuous_on_const cong: continuous_on_cong_simp)
  moreover have "continuous_on {0..} (p s s')"
  proof (subst continuous_on_cong[OF refl p_eq])
    let ?I = "λt. escape_rate s * exp (escape_rate s * t) * (LINT s''|J s. p s'' s' t)"
    show "continuous_on {0..} (λt. (of_bool (s = s') + (LBINT u:{0..t}. ?I u)) / exp (t * escape_rate s))"
    proof (intro continuous_intros continuous_on_LBINT[THEN continuous_on_subset])
      fix t :: real assume t: "0  t"
      then have "0  x  x  t  exp (x * escape_rate s) * (LINT s''|J s. p s'' s' x)  exp (t * escape_rate s) * 1" for x
        by (intro mult_mono) (auto intro!: mult_mono measure_pmf.integral_le_const measure_pmf.integrable_const_bound[of _ 1])
      with t show "set_integrable lborel {0..t} ?I"
        using escape_rate_pos[of s] unfolding set_integrable_def
        by (intro integrableI_bounded_set_indicator[where B="escape_rate s * exp (escape_rate s * t)"])
           (auto simp: field_simps)
    qed auto
  qed simp
  ultimately have "continuous_on ({0..}  {..0}) (p s s')"
    by (intro continuous_on_closed_Un) auto
  also have "{0..}  {..0::real} = UNIV" by auto
  finally show ?thesis
    by (rule continuous_on_subset) simp
qed

lemma p_vector_derivative: ― ‹Backward equation›
  assumes "0  t"
  shows "(p s s' has_vector_derivative (LINT s''|count_space UNIV. R s s'' * p s'' s' t) - escape_rate s * p s s' t)
    (at t within {0..})"
    (is "(_ has_vector_derivative ?A) _")
proof -
  let ?I = "λt. escape_rate s * exp (escape_rate s * t) * (LINT s''|J s. p s'' s' t)"
  let ?p = "λt. (of_bool (s = s') + integral {0..t} ?I) * exp (t *R - escape_rate s)"

  { fix t :: real assume "0  t"
    have "p s s' t = (of_bool (s = s') + (LBINT u:{0..t}. ?I u)) * exp (- t * escape_rate s)"
      using p_eq[OF 0  t, of s s'] by (simp add: exp_minus field_simps)
    also have "(LBINT u:{0..t}. ?I u) = integral {0..t} ?I"
    proof (intro set_borel_integral_eq_integral)
      have "0  x  x  t  exp (x * escape_rate s) * (LINT s''|J s. p s'' s' x)  exp (t * escape_rate s) * 1" for x
        by (intro mult_mono) (auto intro!: mult_mono measure_pmf.integral_le_const measure_pmf.integrable_const_bound[of _ 1])
      with 0t show "set_integrable lborel {0..t} ?I"
        using escape_rate_pos[of s] unfolding set_integrable_def
        by (intro integrableI_bounded_set_indicator[where B="escape_rate s * exp (escape_rate s * t)"])
           (auto simp: field_simps)
    qed
    finally have "p s s' t = ?p t"
      by simp }
  note p_eq = this

  have at_eq: "at t within {0..} = at t within {0 .. t + 1}"
    by (intro at_within_nhd[where S="{..< t+1}"]) auto

  have c_I: "continuous_on {0..t + 1} ?I"
    by (intro continuous_intros continuous_on_LINT_pmf[where B=1] continuous_on_p) simp

  show ?thesis
  proof (subst has_vector_derivative_cong_ev)
    show "F u in nhds t. u  {0..}  p s s' u = ?p u" "p s s' t = ?p t"
      using 0t by (simp_all add: p_eq)
    have "(?p has_vector_derivative escape_rate s * ((LINT s''|J s. p s'' s' t) - p s s' t)) (at t within {0..})"
      unfolding at_eq
      apply (intro refl derivative_eq_intros)
      apply rule
      apply (rule integral_has_vector_derivative[OF c_I])
      apply (simp add: 0  t)
      apply rule
      apply (rule exp_scaleR_has_vector_derivative_right)
      apply (simp add: field_simps exp_minus p_eq 0t split del: split_of_bool)
      done
    also have "escape_rate s * ((LINT s''|J s. p s'' s' t) - p s s' t) =
        (LINT s''|count_space UNIV. R s s'' * p s'' s' t) - escape_rate s * p s s' t"
      using escape_rate_pos[of s]
      by (simp add: measure_pmf_eq_density integral_density J.rep_eq field_simps)
    finally show "(?p has_vector_derivative  ?A) (at t within {0..})" .
  qed
qed

coinductive wf_times :: "real  (real × 'a) stream  bool"
where
  "t < t'  wf_times t' ω  wf_times t ((t', s') ## ω)"

lemma wf_times_simp[simp]: "wf_times t (x ## ω)  t < fst x  wf_times (fst x) ω"
  by (cases x) (subst wf_times.simps; simp)

lemma trace_in_merge_at:
  assumes ω': "wf_times t' ω'"
  shows "trace_in ss t x (merge_at ω t' ω') 
    (if t < t' then trace_in ss t x ω else y. trace_in {y} t' x ω  trace_in ss t y ω')"
    (is "?merge  ?cases")
proof safe
  assume ?merge from this ω' show ?cases
  proof (induction ω"merge_at ω t' ω'" arbitrary: ω ω')
    case (1 j s' y ω'') then show ?case
      by (cases ω) (auto split: if_splits)
  next
    case (2 j x ω' s' ω ω'') then show ?case
      by (cases ω) (auto split: if_splits)
  qed
next
  assume ?cases then show ?merge
  proof (split if_split_asm)
    assume "trace_in ss t x ω" "t < t'" from this ω' show ?thesis
    proof induction
      case 1 then show ?case
        by (cases ω') auto
    qed auto
  next
    assume "y. trace_in {y} t' x ω  trace_in ss t y ω'" "¬ t < t'"
    then obtain y where "trace_in {y} t' x ω" "trace_in ss t y ω'" "t'  t"
      by auto
    from this ω' show ?thesis
      by induction auto
  qed
qed

lemma AE_lim_wf_times: "AE ω in K.lim_stream (t, s). wf_times t ω"
  using AE_lim_stream
proof eventually_elim
  fix ω assume *: "i. snd (((t, s) ## ω) !! i)  DTMC.acc `` {snd (t, s)} 
             snd (ω !! i)  J (snd (((t, s) ## ω) !! i)) 
             fst (((t, s) ## ω) !! i) < fst (ω !! i)"
  have "(t ## smap fst ω) !! i < fst (ω !! i)" for i
    using *[THEN spec, of i] by (cases i) auto
  then show "wf_times t ω"
  proof (coinduction arbitrary: t ω)
    case wf_times from this[THEN spec, of 0] this[THEN spec, of "Suc i" for i]  show ?case
      by (cases ω) auto
  qed
qed

lemma wf_times_shiftD: "wf_times t' (smap (λ(t', y). (t' + t, y)) ω)  wf_times (t' - t) ω"
  apply (coinduction arbitrary: t' t ω)
  subgoal for t' t ω
    apply (cases ω; cases "shd ω")
    apply auto
    subgoal for ω' j x
      by (rule exI[of _ "j + t"]) auto
    done
  done

lemma wf_times_shift[simp]: "wf_times t' (smap (λ(t', y). (t' + t, y)) ω) = wf_times (t' - t) ω"
  using wf_times_shiftD[of "t' - t" "-t" "smap (λ(t', y). (t' + t, y)) ω"]
  by (auto simp: stream.map_comp stream.case_eq_if prod.case_eq_if wf_times_shiftD)

lemma trace_in_unique: "trace_in {y1} t x ω  trace_in {y2} t x ω  y1 = y2"
  by (induction rule: trace_in.induct) auto

lemma trace_at_eq: "trace_in {z} t x ω  trace_at x ω t = z"
  by (induction rule: trace_in.induct) auto

lemma AE_lim_acc: "AE ω in K.lim_stream (t, x). t z. trace_in {z} t x ω  (x, z)  DTMC.acc"
  using AE_lim_stream
proof (eventually_elim, safe)
  fix t' z ω assume *: "i. snd (((t, x) ## ω) !! i)  DTMC.acc `` {snd (t, x)} 
    snd (ω !! i)  J (snd (((t, x) ## ω) !! i))  fst (((t, x) ## ω) !! i) < fst (ω !! i)"
    and t: "trace_in {z} t' x ω"
  define X where "X = DTMC.acc `` {x}"
  have "(x ## smap snd ω) !! i  X" for i
    using *[THEN spec, of i] by (cases i) (auto simp: X_def)
  from t this have "z  X"
  proof induction
    case (1 j y x ω) with "1.prems"[of 0] show ?case
      by simp
  next
    case (2 j y ω x) with "2.prems"[of "Suc i" for i] show ?case
      by simp
  qed
  then show "(x, z)  DTMC.acc"
    by (simp add: X_def)
qed

lemma p_add:
  assumes "0  t" "0  t'"
  shows "p x y (t + t') = (LINT z|count_space (DTMC.acc``{x}). p x z t * p z y t')"
proof -
  interpret L: prob_space "K.lim_stream xy" for xy
    by (rule K.prob_space_lim_stream) simp
  interpret A: sigma_finite_measure "count_space (DTMC.acc``{x})"
    by (intro sigma_finite_measure_count_space_countable DTMC.countable_acc) simp
  interpret LA: pair_sigma_finite "count_space (DTMC.acc``{x})" "K.lim_stream xy" for xy
    by unfold_locales

  have "p x y (t + t') = (+ ω. +ω'. indicator {ωspace T. trace_in {y} (t + t') x ω} (merge_at ω t ω')
    K.lim_stream (t, trace_at x ω t) K.lim_stream (0, x))"
    unfolding p_def L.emeasure_eq_measure[symmetric]
    apply (subst lim_time_split[OF 0  t])
    apply (subst emeasure_bind[OF lim_stream_not_empty measurable_prob_algebraD])
    apply (measurable; fail)
    apply (measurable; fail)
    apply (intro nn_integral_cong)
    apply (subst emeasure_bind[OF lim_stream_not_empty measurable_prob_algebraD])
    apply (measurable; fail)
    apply (measurable; fail)
    apply (simp add: in_space_lim_stream)
    done
  also have " = (+ ω. +ω'. indicator {ωspace T. trace_in {y} (t + t') x ω} (merge_at ω t (smap (λ(t'', s). (t'' + t, s)) ω'))
    K.lim_stream (0, trace_at x ω t) K.lim_stream (0, x))"
    unfolding lim_0[of t] by (subst nn_integral_distr) (measurable; fail)+
  also have " = (+ ω. +ω'. of_bool (zDTMC.acc``{x}. trace_in {z} t x ω  trace_in {y} t' z ω')
    K.lim_stream (0, trace_at x ω t) K.lim_stream (0, x))"
    apply (rule nn_integral_cong_AE)
    using AE_lim_wf_times AE_lim_acc
    apply eventually_elim
    subgoal premises ω for ω
      apply (rule nn_integral_cong_AE)
      using AE_lim_wf_times AE_lim_acc
      apply eventually_elim
      using ω assms
      apply (auto simp add: trace_in_merge_at indicator_eq_1_iff)
      done
    done
  also have " = (+ ω. +ω'. +z. of_bool (trace_in {z} t x ω  trace_in {y} t' z ω')
    count_space (DTMC.acc``{x}) K.lim_stream (0, trace_at x ω t) K.lim_stream (0, x))"
    by (intro nn_integral_cong of_bool_Bex_eq_nn_integral) (auto dest: trace_in_unique)
  also have " = (+ ω. +z. +ω'. of_bool (trace_in {z} t x ω  trace_in {y} t' z ω')
    K.lim_stream (0, trace_at x ω t) count_space (DTMC.acc``{x}) K.lim_stream (0, x))"
    apply (subst LA.Fubini')
    apply (subst measurable_split_conv)
    apply (rule measurable_compose_countable'[OF _ measurable_fst])
    apply (auto simp: DTMC.countable_acc)
    done
  also have " = (+z. + ω. of_bool (trace_in {z} t x ω) * +ω'. of_bool (trace_in {y} t' z ω')
    K.lim_stream (0, z) K.lim_stream (0, x) count_space (DTMC.acc``{x}))"
    apply (subst LA.Fubini')
    apply (subst measurable_split_conv)
    apply (rule measurable_compose_countable'[OF _ measurable_fst])
    apply (rule nn_integral_measurable_subprob_algebra2)
    apply (measurable; fail)
    apply (rule measurable_prob_algebraD)
    apply (auto simp: DTMC.countable_acc trace_at_eq intro!: nn_integral_cong)
    done
  also have " = (+z. (+ ω. of_bool (trace_in {z} t x ω)K.lim_stream (0, x)) *
      (+ω'. of_bool (trace_in {y} t' z ω') K.lim_stream (0, z)) count_space (DTMC.acc``{x}))"
     by (auto intro!: nn_integral_cong simp: nn_integral_multc)
  also have " = (+z. ennreal (p x z t) * ennreal (p z y t') count_space (DTMC.acc``{x}))"
    unfolding p_def L.emeasure_eq_measure[symmetric]
    by (auto intro!: nn_integral_cong arg_cong2[where f="(*)"]
             simp: nn_integral_indicator[symmetric] simp del: nn_integral_indicator )
  finally have "(+z. p x z t * p z y t' count_space (DTMC.acc``{x})) = p x y (t + t')"
    by (simp add: ennreal_mult)
  then show ?thesis
    by (subst (asm) nn_integral_eq_integrable) auto
qed

end

end