Theory Classifying_Markov_Chain_States

section ‹Classifying Markov Chain States›

theory Classifying_Markov_Chain_States
  imports
    "HOL-Computational_Algebra.Group_Closure"
    Discrete_Time_Markov_Chain
begin

lemma eventually_mult_Gcd:
  fixes S :: "nat set"
  assumes S: "s t. s  S  t  S  s + t  S"
  assumes s: "s  S" "s > 0"
  shows "eventually (λm. m * Gcd S  S) sequentially"
proof -
  define T where "T = insert 0 (int ` S)"
  with s S have "int s  T" "0  T" and T: "r  T  t  T  r + t  T" for r t
    by (auto simp del: of_nat_add simp add: of_nat_add [symmetric])
  have "Gcd T  group_closure T"
    by (rule Gcd_in_group_closure)
  also have "group_closure T = {s - t | s t. s  T  t  T}"
  proof (auto intro: group_closure.base group_closure.diff)
    fix x assume "x  group_closure T"
    then show "s t. x = s - t  s  T  t  T"
    proof induction
      case (base x) with 0  T show ?case
        apply (rule_tac x=x in exI)
        apply (rule_tac x=0 in exI)
        apply auto
        done
    next
      case (diff x y)
      then obtain a b c d where
        "a  T" "b  T" "x = a - b"
        "c  T" "d  T" "y = c - d"
        by auto
      then show ?case
        apply (rule_tac x="a + d" in exI)
        apply (rule_tac x="b + c" in exI)
        apply (auto intro: T)
        done
    qed
  qed
  finally obtain s' t' :: int
    where "s'  T" "t'  T" "Gcd T = s' - t'"
    by blast
  moreover define s and t where "s = nat s'" and "t = nat t'"
  moreover have "int (Gcd S) = - int t  S  {0}  t = 0"
    by auto (metis Gcd_dvd_nat dvd_0_right dvd_antisym nat_int nat_zminus_int) 
  ultimately have 
    st: "s = 0  s  S" "t = 0  t  S" and Gcd_S: "Gcd S = s - t"
    using T_def by safe simp_all
  with s
  have "t < s"
    by (rule_tac ccontr) auto

  { fix s n have "0 < n  s  S  n * s  S"
    proof (induct n)
      case (Suc n) then show ?case
        by (cases n) (auto intro: S)
    qed simp }
  note cmult_S = this

  show ?thesis
    unfolding eventually_sequentially
  proof cases
    assume "s = 0  t = 0"
    with st Gcd_S s have *: "Gcd S  S"
      by (auto simp: int_eq_iff)
    then show "N. nN. n * Gcd S  S" by (auto intro!: exI[of _ 1] cmult_S)
  next
    assume "¬ (s = 0  t = 0)"
    with st have "s  S" "t  S" "t  0" by auto
    then have "Gcd S dvd t" by auto
    then obtain a where a: "t = Gcd S * a" ..
    with t  0 have "0 < a" by auto

    show "N. nN. n * Gcd S  S"
    proof (safe intro!: exI[of _ "a * a"])
      fix n
      define m where "m = (n - a * a) div a"
      define r where "r = (n - a * a) mod a"
      with 0 < a have "r < a" by simp
      moreover define am where "am = a + m"
      ultimately have "r < am" by simp
      assume "a * a  n" then have n: "n = a * a + (m * a + r)"
        unfolding m_def r_def by simp
      have "n * Gcd S = am * t + r * Gcd S"
        unfolding n a by (simp add: field_simps am_def)
      also have " = r * s + (am - r) * t"
        unfolding Gcd S = s - t
        using t < s r < am by (simp add: field_simps diff_mult_distrib2)
      also have "  S"
        using s  S t  S r < am
        by (cases "r = 0") (auto intro!: cmult_S S)
      finally show "n * Gcd S  S" .
    qed
  qed
qed

context MC_syntax
begin

subsection ‹Expected number of visits›

definition "G s t = (+ω. scount (HLD {t}) (s ## ω) T s)"

lemma G_eq: "G s t = (+ω. emeasure (count_space UNIV) {i. (s ## ω) !! i = t} T s)"
  by (simp add: G_def scount_eq_emeasure HLD_iff)

definition "p s t n = 𝒫(ω in T s. (s ## ω) !! n = t)"

definition "gf_G s t z = (n. p s t n *R z ^ n)"

definition "convergence_G s t z  summable (λn. p s t n * norm z ^ n)"

lemma p_nonneg[simp]: "0  p x y n"
  by (simp add: p_def)

lemma p_le_1: "p x y n  1"
  by (simp add: p_def)

lemma p_x_x_0[simp]: "p x x 0 = 1"
  by (simp add: p_def T.prob_space del: space_T)

lemma p_0: "p x y 0 = (if x = y then 1 else 0)"
  by (simp add: p_def T.prob_space del: space_T)

lemma p_in_reachable: assumes "(x, y)  (SIGMA x:UNIV. K x)*" shows "p x y n = 0"
  unfolding p_def
proof (rule T.prob_eq_0_AE)
  from AE_T_reachable show "AE ω in T x. (x ## ω) !! n  y"
  proof eventually_elim
    fix ω assume "alw (HLD ((SIGMA ω:UNIV. K ω)* `` {x})) ω"
    then have "alw (HLD (- {y})) ω"
      using assms by (auto intro: alw_mono simp: HLD_iff)
    then show "(x ## ω) !! n  y"
      using assms by (cases n) (auto simp: alw_HLD_iff_streams streams_iff_snth)
  qed
qed

lemma p_Suc: "ennreal (p x y (Suc n)) = (+ w. p w y n K x)"
  unfolding p_def T.emeasure_eq_measure[symmetric] by (subst emeasure_Collect_T) simp_all

lemma p_Suc':
  "p x y (Suc n) = (x'. p x' y n K x)"
  using p_Suc[of x y n]
  by (subst (asm) nn_integral_eq_integral)
     (auto simp: p_le_1 intro!: measure_pmf.integrable_const_bound[where B=1])

lemma p_add: "p x y (n + m) = (+ w. p x w n * p w y m count_space UNIV)"
proof (induction n arbitrary: x)
  case 0
  have [simp]: "w. (if x = w then 1 else 0) * p w y m = ennreal (p x y m) * indicator {x} w"
    by auto
  show ?case
    by (simp add: p_0 one_ennreal_def[symmetric] max_def)
next
  case (Suc n)
  define X where "X = (SIGMA x:UNIV. K x)* `` K x"
  then have X: "countable X"
    by (blast intro: countable_Image countable_reachable countable_set_pmf)

  then interpret X: sigma_finite_measure "count_space X"
    by (rule sigma_finite_measure_count_space_countable)
  interpret XK: pair_sigma_finite "K x" "count_space X"
    by unfold_locales

  have "ennreal (p x y (Suc n + m)) = (+t. (+w. p t w n * p w y m count_space UNIV) K x)"
    by (simp add: p_Suc Suc)
  also have " = (+t. (+w. ennreal (p t w n * p w y m) * indicator X w count_space UNIV) K x)"
    by (auto intro!: nn_integral_cong_AE simp: AE_measure_pmf_iff AE_count_space Image_iff p_in_reachable X_def             split: split_indicator)
  also have " = (+t. (+w. p t w n * p w y m count_space X) K x)"
    by (subst nn_integral_restrict_space[symmetric]) (simp_all add: restrict_count_space)
  also have " = (+w. (+t. p t w n * p w y m K x) count_space X)"
    apply (rule XK.Fubini'[symmetric])
    unfolding measurable_split_conv
    apply (rule measurable_compose_countable'[OF _ measurable_snd X])
    apply (rule measurable_compose[OF measurable_fst])
    apply simp
    done
  also have " = (+w. (+t. ennreal (p t w n * p w y m) * indicator X w K x) count_space UNIV)"
    by (simp add: nn_integral_restrict_space[symmetric] restrict_count_space nn_integral_multc)
  also have " = (+w. (+t. ennreal (p t w n * p w y m) K x) count_space UNIV)"
    by (auto intro!: nn_integral_cong_AE simp: AE_measure_pmf_iff AE_count_space Image_iff p_in_reachable X_def             split: split_indicator)
  also have " = (+w. (+t. p t w n K x) * p w y m count_space UNIV)"
    by (simp add: nn_integral_multc[symmetric] ennreal_mult)
  finally show ?case
    by (simp add: ennreal_mult p_Suc)
qed

lemma prob_reachable_le:
  assumes [simp]: "m  n"
  shows "p x y m * p y w (n - m)  p x w n"
proof -
  have "p x y m * p y w (n - m) = (+y'. ennreal (p x y m * p y w (n - m)) * indicator {y} y' count_space UNIV)"
    by simp
  also have "  p x w (m + (n - m))"
    by (subst p_add)
       (auto intro!: nn_integral_mono split: split_indicator simp del: nn_integral_indicator_singleton)
  finally show ?thesis
    by simp
qed

lemma G_eq_suminf: "G x y = (i. ennreal (p x y i))"
proof -
  have *: "i ω. indicator {ω  space S. (x ## ω) !! i = y} ω = indicator {i. (x ## ω) !! i = y} i"
    by (auto simp: space_stream_space split: split_indicator)

  have "G x y = (+ ω. (i. indicator {ωspace (T x). (x ## ω) !! i = y} ω) T x)"
    unfolding G_eq by (simp add: nn_integral_count_space_nat[symmetric] *)
  also have " = (i. ennreal (p x y i))"
    by (simp add: T.emeasure_eq_measure[symmetric] p_def nn_integral_suminf)
  finally show ?thesis .
qed

lemma G_eq_real_suminf:
  "convergence_G x y (1::real)  G x y = ennreal (i. p x y i)"
  unfolding G_eq_suminf
  by (intro suminf_ennreal ennreal_suminf_neq_top p_nonneg)
     (auto simp: convergence_G_def p_def)

lemma convergence_norm_G:
  "convergence_G x y z  summable (λn. p x y n * norm z ^ n)"
  unfolding convergence_G_def .

lemma convergence_G:
  "convergence_G x y (z::'a::{banach, real_normed_div_algebra})  summable (λn. p x y n *R z ^ n)"
  unfolding convergence_G_def
  by (rule summable_norm_cancel) (simp add: abs_mult norm_power)

lemma convergence_G_less_1:
  fixes z :: "_ :: {banach, real_normed_field}"
  assumes z: "norm z < 1" shows "convergence_G x y z"
  unfolding convergence_G_def
proof (rule summable_comparison_test)
  have "n. p x y n * norm (z ^ n)  1 * norm (z ^ n)"
    by (intro mult_right_mono p_le_1) simp_all
  then show "N. nN. norm (p x y n * norm z ^ n)  norm z ^ n"
    by (simp add: norm_power)
qed (simp add: z summable_geometric)

lemma lim_gf_G: "((λz. ennreal (gf_G x y z))  G x y) (at_left (1::real))"
  unfolding gf_G_def G_eq_suminf real_scaleR_def
  by (intro power_series_tendsto_at_left p_nonneg p_le_1 summable_power_series)

subsection ‹Reachability probability›

definition "u x y n = 𝒫(ω in T x. ev_at (HLD {y}) n ω)"

definition "U s t = 𝒫(ω in T s. ev (HLD {t}) ω)"

definition "gf_U x y z = (n. u x y n *R z ^ Suc n)"

definition "f x y n = 𝒫(ω in T x. ev_at (HLD {y}) n (x ## ω))"

definition "F s t = 𝒫(ω in T s. ev (HLD {t}) (s ## ω))"

definition "gf_F x y z = (n. f x y n * z ^ n)"

lemma f_Suc: "x  y  f x y (Suc n) = u x y n"
  by (simp add: u_def f_def)

lemma f_Suc_eq: "f x x (Suc n) = 0"
  by (simp add: f_def)

lemma f_0: "f x y 0 = (if x = y then 1 else 0)"
  using T.prob_space by (simp add: f_def)

lemma shows u_nonneg: "0  u x y n" and u_le_1: "u x y n  1"
  by (simp_all add: u_def)

lemma shows f_nonneg: "0  f x y n" and f_le_1: "f x y n  1"
  by (simp_all add: f_def)

lemma U_nonneg[simp]: "0  U x y"
  by (simp add: U_def)

lemma U_le_1: "U s t  1"
  by (auto simp add: U_def intro!: antisym)

lemma U_cases: "U s s = 1  U s s < 1"
  by (auto simp add: U_def intro!: antisym)

lemma u_sums_U: "u x y sums U x y"
  unfolding u_def[abs_def] U_def ev_iff_ev_at by (intro T.prob_sums) (auto intro: ev_at_unique)

lemma gf_U_eq_U: "gf_U x y 1 = U x y"
  using u_sums_U[THEN sums_unique] by (simp add: gf_U_def U_def)

lemma f_sums_F: "f x y sums F x y"
  unfolding f_def[abs_def] F_def ev_iff_ev_at
  by (intro T.prob_sums) (auto intro: ev_at_unique)

lemma F_nonneg[simp]: "0  F x y"
  by (auto simp: F_def)

lemma F_le_1: "F x y  1"
  by (simp add: F_def)

lemma gf_F_eq_F: "gf_F x y 1 = F x y"
  using f_sums_F[THEN sums_unique] by (simp add: gf_F_def F_def)

lemma gf_F_le_1:
  fixes z :: real
  assumes z: "0  z" "z  1"
  shows "gf_F x y z  1"
proof -
  have "gf_F x y z  gf_F x y 1"
    using z unfolding gf_F_def
    by (intro suminf_le[OF _ summable_comparison_test[OF _ sums_summable[OF f_sums_F[of x y]]]] mult_left_mono allI f_nonneg)
       (simp_all add: power_le_one f_nonneg mult_right_le_one_le f_le_1 sums_summable[OF f_sums_F[of x y]])
  also have "  1"
    by (simp add: gf_F_eq_F F_def)
  finally show ?thesis .
qed

lemma u_le_p: "u x y n  p x y (Suc n)"
  unfolding u_def p_def by (auto intro!: T.finite_measure_mono dest: ev_at_HLD_imp_snth)

lemma f_le_p: "f x y n  p x y n"
  unfolding f_def p_def by (auto intro!: T.finite_measure_mono dest: ev_at_HLD_imp_snth)

lemma convergence_norm_U:
  fixes z :: "_ :: real_normed_div_algebra"
  assumes z: "convergence_G x y z"
  shows "summable (λn. u x y n * norm z ^ Suc n)"
  using summable_ignore_initial_segment[OF convergence_norm_G[OF z], of 1]
  by (rule summable_comparison_test[rotated])
     (auto simp add: u_nonneg abs_mult intro!: exI[of _ 0] mult_right_mono u_le_p)

lemma convergence_norm_F:
  fixes z :: "_ :: real_normed_div_algebra"
  assumes z: "convergence_G x y z"
  shows "summable (λn. f x y n * norm z ^ n)"
  using convergence_norm_G[OF z]
  by (rule summable_comparison_test[rotated])
     (auto simp add: f_nonneg abs_mult intro!: exI[of _ 0] mult_right_mono f_le_p)

lemma gf_G_nonneg:
  fixes z :: real
  shows "0  z  z < 1  0  gf_G x y z"
  unfolding gf_G_def
  by (intro suminf_nonneg convergence_G convergence_G_less_1) simp_all

lemma gf_F_nonneg:
  fixes z :: real
  shows "0  z  z < 1  0  gf_F x y z"
  unfolding gf_F_def
  using convergence_norm_F[OF convergence_G_less_1, of z x y]
  by (intro suminf_nonneg) (simp_all add: f_nonneg)

lemma convergence_U:
  fixes z :: "_ :: banach"
  shows "convergence_G x y z  summable (λn. u x y n * z ^ Suc n)"
  by (rule summable_norm_cancel)
     (auto simp add: abs_mult u_nonneg power_abs dest!: convergence_norm_U)

lemma p_eq_sum_p_u: "p x y (Suc n) = (in. p y y (n - i) * u x y i)"
proof -
  have "ω. ω !! n = y  (i. i  n  ev_at (HLD {y}) i ω)"
  proof (induction n)
    case (Suc n)
    then obtain i where "i  n" "ev_at (HLD {y}) i (stl ω)"
      by auto
    then show ?case
      by (auto intro!: exI[of _ "if HLD {y} ω then 0 else Suc i"])
  qed (simp add: HLD_iff)
  then have "p x y (Suc n) = (in. 𝒫(ω in T x. ev_at (HLD {y}) i ω  ω !! n = y))"
    unfolding p_def by (intro T.prob_sum) (auto intro: ev_at_unique)
  also have " = (in. p y y (n - i) * u x y i)"
  proof (intro sum.cong refl)
    fix i assume i: "i  {.. n}"
    then have "ω. (Suc i  n  ω !! (n - Suc i) = y)  ((y ## ω) !! (n - i) = y)"
      by (auto simp: Stream_snth diff_Suc split: nat.split)
    from i have "i  n" by auto
    then have "𝒫(ω in T x. ev_at (HLD {y}) i ω  ω !! n = y) =
      (ω'. 𝒫(ω in T y. (y ## ω) !! (n - i) = y) *
        indicator {ω'space (T x). ev_at (HLD {y}) i ω' } ω' T x)"
      by (subst prob_T_split[where n="Suc i"])
         (auto simp: ev_at_shift ev_at_HLD_single_imp_snth shift_snth diff_Suc
               split: split_indicator nat.split intro!: Bochner_Integration.integral_cong arg_cong2[where f=measure]
               simp del: stake.simps integral_mult_right_zero)
    then show "𝒫(ω in T x. ev_at (HLD {y}) i ω  ω !! n = y) = p y y (n - i) * u x y i"
      by (simp add: p_def u_def)
  qed
  finally show ?thesis .
qed

lemma p_eq_sum_p_f: "p x y n = (in. p y y (n - i) * f x y i)"
  by (cases n)
     (simp_all del: sum.atMost_Suc
               add: f_0 p_0 p_eq_sum_p_u atMost_Suc_eq_insert_0 zero_notin_Suc_image sum.reindex
                    f_Suc f_Suc_eq)

lemma gf_G_eq_gf_F:
  assumes z: "norm z < 1"
  shows "gf_G x y z = gf_F x y z * gf_G y y z"
proof -
  have "gf_G x y z = (n. in. p y y (n - i) * f x y i * z^n)"
    by (simp add: gf_G_def p_eq_sum_p_f[of x y] sum_distrib_right)
  also have " = (n. in. (f x y i * z^i) * (p y y (n - i) * z^(n - i)))"
    by (intro arg_cong[where f=suminf] sum.cong ext atLeast0AtMost[symmetric])
       (simp_all add: power_add[symmetric])
  also have " = (n. f x y n * z^n) * (n. p y y n * z^n)"
    using convergence_norm_F[OF convergence_G_less_1[OF z]] convergence_norm_G[OF convergence_G_less_1[OF z]]
    by (intro Cauchy_product[symmetric]) (auto simp: f_nonneg abs_mult power_abs)
  also have " = gf_F x y z * gf_G y y z"
    by (simp add: gf_F_def gf_G_def)
  finally show ?thesis .
qed

lemma gf_G_eq_gf_U:
  fixes z :: "'z :: {banach, real_normed_field}"
  assumes z: "convergence_G x x z"
  shows "gf_G x x z = 1 / (1 - gf_U x x z)" "gf_U x x z  1"
proof -
  { fix n
    have "p x x (Suc n) *R z^Suc n = (in. (p x x (n - i) * u x x i) *R z^Suc n)"
      unfolding scaleR_sum_left[symmetric] by (simp add: p_eq_sum_p_u)
    also have " = (in. (u x x i *R z^Suc i) * (p x x (n - i) *R z^(n - i)))"
      by (intro sum.cong refl) (simp add: field_simps power_diff cong: disj_cong)
    finally have "p x x (Suc n) *R z^(Suc n) = (in. (u x x i *R z^Suc i) * (p x x (n - i) *R z^(n - i)))"
      unfolding atLeast0AtMost . }
  note gfs_Suc_eq = this

  have "gf_G x x z = 1 + (n. p x x (Suc n) *R z^(Suc n))"
    unfolding gf_G_def
    by (subst suminf_split_initial_segment[OF convergence_G[OF z], of 1]) simp
  also have " = 1 + (n. in. (u x x i *R z^Suc i) * (p x x (n - i) *R z^(n - i)))"
    unfolding gfs_Suc_eq ..
  also have " = 1 + gf_U x x z * gf_G x x z"
    unfolding gf_U_def gf_G_def
    by (subst Cauchy_product)
       (auto simp: u_nonneg norm_power simp del: power_Suc
             intro!: z convergence_norm_G convergence_norm_U)
  finally show "gf_G x x z = 1 / (1 - gf_U x x z)" "gf_U x x z  1"
    apply -
    apply (cases "gf_U x x z = 1")
    apply (auto simp add: field_simps)
    done
qed

lemma gf_U: "(gf_U x y  U x y) (at_left 1)"
proof -
  have "((λz. ennreal (n. u x y n * z ^ n))  (n. ennreal (u x y n))) (at_left 1)"
    using u_le_1 u_nonneg by (intro power_series_tendsto_at_left summable_power_series)
  also have "(n. ennreal (u x y n)) = ennreal (suminf (u x y))"
    by (intro u_nonneg suminf_ennreal ennreal_suminf_neq_top sums_summable[OF u_sums_U])
  also have "suminf (u x y) = U x y"
    using u_sums_U by (rule sums_unique[symmetric])
  finally have "((λz. n. u x y n * z ^ n)  U x y) (at_left 1)"
    by (rule tendsto_ennrealD)
       (auto simp: u_nonneg u_le_1 intro!: suminf_nonneg summable_power_series eventually_at_left_1)
  then have "((λz. z * (n. u x y n * z ^ n))  1 * U x y) (at_left 1)"
    by (intro tendsto_intros) simp
  then have "((λz. n. u x y n * z ^ Suc n)  1 * U x y) (at_left 1)"
    apply (rule filterlim_cong[OF refl refl, THEN iffD1, rotated])
    apply (rule eventually_at_left_1)
    apply (subst suminf_mult[symmetric])
    apply (auto intro!: summable_power_series u_le_1 u_nonneg)
    apply (simp add: field_simps)
    done
  then show ?thesis
    by (simp add: gf_U_def[abs_def] U_def)
qed

lemma gf_U_le_1: assumes z: "0 < z" "z < 1" shows "gf_U x y z  (1::real)"
proof -
  note u = u_sums_U[of x y, THEN sums_summable]
  have "gf_U x y z  gf_U x y 1"
    using z
    unfolding gf_U_def real_scaleR_def
    by (intro suminf_le allI mult_mono power_mono summable_comparison_test_ev[OF _ u] always_eventually)
       (auto simp: u_nonneg intro!: mult_left_le mult_le_one power_le_one)
  also have "  1"
    unfolding gf_U_eq_U by (rule U_le_1)
  finally show ?thesis .
qed

lemma gf_F: "(gf_F x y  F x y) (at_left 1)"
proof -
  have "((λz. ennreal (n. f x y n * z ^ n))  (n. ennreal (f x y n))) (at_left 1)"
    using f_le_1 f_nonneg by (intro power_series_tendsto_at_left summable_power_series)
  also have "(n. ennreal (f x y n)) = ennreal (suminf (f x y))"
    by (intro f_nonneg suminf_ennreal ennreal_suminf_neq_top sums_summable[OF f_sums_F])
  also have "suminf (f x y) = F x y"
    using f_sums_F by (rule sums_unique[symmetric])
  finally have "((λz. n. f x y n * z ^ n)  F x y) (at_left 1)"
    by (rule tendsto_ennrealD)
       (auto simp: f_nonneg f_le_1 intro!: suminf_nonneg summable_power_series eventually_at_left_1)
  then show ?thesis
    by (simp add: gf_F_def[abs_def] F_def)
qed

lemma U_bounded: "0  U x y" "U x y  1"
  unfolding U_def by simp_all

subsection ‹Recurrent states›

definition recurrent :: "'s  bool" where
  "recurrent s  (AE ω in T s. ev (HLD {s}) ω)"

lemma recurrent_iff_U_eq_1: "recurrent s  U s s = 1"
    unfolding recurrent_def U_def by (subst T.prob_Collect_eq_1) simp_all

definition "H s t = 𝒫(ω in T s. alw (ev (HLD {t})) ω)"

lemma H_eq:
  "recurrent s  H s s = 1"
  "¬ recurrent s  H s s = 0"
  "H s t = U s t * H t t"
proof -
  define H' where "H' t n = {ωspace S. enat n  scount (HLD {t::'s}) ω}" for t n
  have [measurable]: "y n. H' y n  sets S"
    by (simp add: H'_def)
  let ?H' = "λs t n. measure (T s) (H' t n)"
  { fix x y :: 's and ω
    have "Suc 0  scount (HLD {y}) ω  ev (HLD {y}) ω"
      using scount_eq_0_iff[of "HLD {y}" ω]
      by (cases "scount (HLD {y}) ω" rule: enat_coexhaust)
         (auto simp: not_ev_iff[symmetric] eSuc_enat[symmetric] enat_0 HLD_iff[abs_def]) }
  then have H'_1: "x y. ?H' x y 1 = U x y"
    unfolding H'_def U_def by simp

  { fix n and x y :: 's
    let ?U = "(not (HLD {y}) suntil (HLD {y} aand nxt (λω. enat n  scount (HLD {y}) ω)))"
    { fix ω
      have "enat (Suc n)  scount (HLD {y}) ω  ?U ω"
      proof
        assume "enat (Suc n)  scount (HLD {y}) ω"
        with scount_eq_0_iff[of "HLD {y}" ω] have "ev (HLD {y}) ω" "enat (Suc n)  scount (HLD {y}) ω"
          by (auto simp add: not_ev_iff[symmetric] eSuc_enat[symmetric])
        then show "?U ω"
          by (induction rule: ev_induct_strong)
             (auto simp: scount_simps eSuc_enat[symmetric] intro: suntil.intros)
      next
        assume "?U ω" then show "enat (Suc n)  scount (HLD {y}) ω"
          by induction (auto simp: scount_simps  eSuc_enat[symmetric])
      qed }
    then have "emeasure (T x) (H' y (Suc n)) = emeasure (T x) {ωspace (T x). ?U ω}"
      by (simp add: H'_def)
    also have " = U x y * ?H' y y n"
      by (subst emeasure_suntil_HLD) (simp_all add: T.emeasure_eq_measure U_def H'_def ennreal_mult)
    finally have "?H' x y (Suc n) = U x y * ?H' y y n"
      by (simp add: T.emeasure_eq_measure) }
  note H'_Suc = this

  { fix m and x :: 's
    have "?H' x x (Suc m) = U x x^Suc m"
      using H'_1 H'_Suc by (induct m) auto }
  note H'_eq = this

  { fix x y
    have "?H' x y  measure (T x) (i. H' y i)"
    proof (rule T.finite_Lim_measure_decseq)
      show "range (H' y)  T.events x"
        by auto
    next
      show "decseq (H' y)"
        by (rule antimonoI) (simp add: subset_eq H'_def order_subst2)
    qed
    also have "(i. H' y i) = {ωspace (T x). alw (ev (HLD {y})) ω}"
      by (auto simp: H'_def scount_infinite_iff[symmetric]) (metis Suc_ile_eq enat.exhaust neq_iff)
    finally have "?H' x y  H x y"
      unfolding H_def . }
  note H'_lim = this

  from H'_lim[of s s, THEN LIMSEQ_Suc]
  have "(λn. U s s ^ Suc n)  H s s"
    by (simp add: H'_eq)
  then have lim_H: "(λn. U s s ^ n)  H s s"
    by (rule LIMSEQ_imp_Suc)

  have "U s s < 1  (λn. U s s ^ n)  0"
    by (rule LIMSEQ_realpow_zero) (simp_all add: U_def)
  with lim_H have "U s s < 1  H s s = 0"
    by (blast intro: LIMSEQ_unique)
  moreover have "U s s = 1  (λn. U s s ^ n)  1"
    by simp
  with lim_H have "U s s = 1  H s s = 1"
    by (blast intro: LIMSEQ_unique)
  moreover note recurrent_iff_U_eq_1 U_cases
  ultimately show "recurrent s  H s s = 1" "¬ recurrent s  H s s = 0"
    by (metis one_neq_zero)+

  from H'_lim[of s t, THEN LIMSEQ_Suc] H'_Suc[of s]
  have "(λn. U s t * ?H' t t n)  H s t"
    by simp
  moreover have "(λn. U s t * ?H' t t n)  U s t * H t t"
    by (intro tendsto_intros H'_lim)
  ultimately show "H s t = U s t * H t t"
    by (blast intro: LIMSEQ_unique)
qed

lemma recurrent_iff_G_infinite: "recurrent x  G x x = "
proof -
  have "((λz. ennreal (gf_G x x z))  G x x) (at_left 1)"
    by (rule lim_gf_G)
  then have G: "((λz. ennreal (1 / (1 - gf_U x x z)))  G x x) (at_left (1::real))"
    apply (rule filterlim_cong[OF refl refl, THEN iffD1, rotated])
    apply (rule eventually_at_left_1)
    apply (subst gf_G_eq_gf_U)
    apply (rule convergence_G_less_1)
    apply simp
    apply simp
    done

  { fix z :: real assume z: "0 < z" "z < 1"
    have 1: "summable (u x x)"
      using u_sums_U by (rule sums_summable)
    have "gf_U x x z  1"
      using gf_G_eq_gf_U[OF convergence_G_less_1[of z]] z by simp
    moreover
    have "gf_U x x z  U x x"
      unfolding gf_U_def gf_U_eq_U[symmetric]
      using z
      by (intro suminf_le)
         (auto simp add: 1 convergence_U convergence_G_less_1 u_nonneg simp del: power_Suc
               intro!: mult_right_le_one_le power_le_one)
    ultimately have "gf_U x x z < 1"
      using U_bounded[of x x] by simp }
  note strict = this

  { assume "U x x = 1"
    moreover have "((λxa. 1 - gf_U x x xa :: real)  1 - U x x) (at_left 1)"
      by (intro tendsto_intros gf_U)
    moreover have "eventually (λz. gf_U x x z < 1) (at_left (1::real))"
      by (auto intro!: eventually_at_left_1 strict simp: U x x = 1 gf_U_eq_U)
    ultimately have "((λz. ennreal (1 / (1 - gf_U x x z)))  top) (at_left 1)"
      unfolding ennreal_tendsto_top_eq_at_top
      by (intro LIM_at_top_divide[where a=1] tendsto_const zero_less_one)
         (auto simp: field_simps)
    with G have "G x x = top"
      by (rule tendsto_unique[rotated]) simp }
  moreover
  { assume "U x x < 1"
    then have "((λxa. ennreal (1 / (1 - gf_U x x xa)))  1 / (1 - U x x)) (at_left 1)"
      by (intro tendsto_intros gf_U tendsto_ennrealI) simp
    from tendsto_unique[OF _ G this] have "G x x  "
      by simp }
  ultimately show ?thesis
    using U_cases recurrent_iff_U_eq_1 by auto
qed

definition communicating :: "('s × 's) set" where
  "communicating = acc  acc¯"

definition essential_class :: "'s set  bool" where
  "essential_class C  C  UNIV // communicating  acc `` C  C"

lemma accI_U:
  assumes "0 < U x y" shows "(x, y)  acc"
proof (rule ccontr)
  assume *: "(x, y)  acc"

  { fix ω assume "ev (HLD {y}) ω" "alw (HLD (acc `` {x})) ω" from this * have False
      by induction (auto simp: HLD_iff) }
  with AE_T_reachable[of x] have "U x y = 0"
    unfolding U_def by (intro T.prob_eq_0_AE) auto
  with 0 < U x y show False by auto
qed

lemma accD_pos:
  assumes "(x, y)  acc"
  shows "n. 0 < p x y n"
using assms proof induction
  case base with T.prob_space[of x] show ?case
    by (auto intro!: exI[of _ 0])
next
  have [simp]: "x y. (if x = y then 1 else 0::real) = indicator {y} x"
    by simp
  case (step w y)
  then obtain n where "0 < p x w n" and "0 < pmf (K w) y"
    by (auto simp: set_pmf_iff less_le)
  then have "0 < p x w n * pmf (K w) y"
    by (intro mult_pos_pos)
  also have "  p x w n * p w y (Suc 0)"
    by (simp add: p_Suc' p_0 pmf.rep_eq)
  also have "  p x y (Suc n)"
    using prob_reachable_le[of n "Suc n" x w y] by simp
  finally show ?case ..
qed

lemma accI_pos: "0 < p x y n  (x, y)  acc"
proof (induct n arbitrary: x)
  case (Suc n)
  then have less: "0 < (x'. p x' y n K x)"
    by (simp add: p_Suc')
  have "x'K x. 0 < p x' y n"
  proof (rule ccontr)
    assume "¬ ?thesis"
    then have "AE x' in K x. p x' y n = 0"
      by (simp add: AE_measure_pmf_iff less_le)
    then have "(x'. p x' y n K x) = (x'. 0 K x)"
      by (intro integral_cong_AE) simp_all
    with less show False by simp
  qed
  with Suc show ?case
    by (auto intro: converse_rtrancl_into_rtrancl)
qed (simp add: p_0 split: if_split_asm)

lemma recurrent_iffI_communicating:
  assumes "(x, y)  communicating"
  shows "recurrent x  recurrent y"
proof -
  from assms obtain n m where "0 < p x y n" "0 < p y x m"
    by (force simp: communicating_def dest: accD_pos)
  moreover
  { fix x y n m assume "0 < p x y n" "0 < p y x m" "G y y = "
    then have " = ennreal (p x y n * p y x m) * G y y"
      by (auto intro: mult_pos_pos simp: ennreal_mult_top)
    also have "ennreal (p x y n * p y x m) * G y y = (i. ennreal (p x y n * p y x m) * p y y i)"
      unfolding G_eq_suminf by (rule ennreal_suminf_cmult[symmetric])
    also have "  (i. ennreal (p x x (n + i + m)))"
    proof (intro suminf_le allI)
      fix i
      have "(p x y n * p y y ((n + i) - n)) * p y x ((n + i + m) - (n + i))  p x y (n + i) * p y x ((n + i + m) - (n + i))"
        by (intro mult_right_mono prob_reachable_le) simp_all
      also have "  p x x (n + i + m)"
         by (intro prob_reachable_le) simp_all
      finally show "ennreal (p x y n * p y x m) * p y y i  ennreal (p x x (n + i + m))"
        by (simp add: ac_simps ennreal_mult'[symmetric])
    qed auto
    also have "  (i. ennreal (p x x (i + (n + m))))"
      by (simp add: ac_simps)
    also have "  (i. ennreal (p x x i))"
      by (subst suminf_offset[of "λi. ennreal (p x x i)" "n + m"]) auto
    also have "  G x x"
      unfolding G_eq_suminf by (auto intro!: suminf_le_pos)
    finally have "G x x = "
      by (simp add: top_unique) }
  ultimately show ?thesis
    using recurrent_iff_G_infinite by blast
qed

lemma recurrent_acc:
  assumes "recurrent x" "(x, y)  acc"
  shows "U y x = 1" "H y x = 1" "recurrent y" "(x, y)  communicating"
proof -
  { fix w y assume step: "(x, w)  acc" "y  K w" "U w x = 1" "H w x = 1" "recurrent w" "x  y"
    have "measure (K w) UNIV = U w x"
      using step measure_pmf.prob_space[of "K w"] by simp
    also have " = (v. indicator {x} v + U v x * indicator (- {x}) v K w)"
      unfolding U_def
      by (subst prob_T)
         (auto intro!: Bochner_Integration.integral_cong arg_cong2[where f=measure] AE_I2
               simp: ev_Stream T.prob_eq_1 split: split_indicator)
    also have " = measure (K w) {x} + (v. U v x * indicator (- {x}) v K w)"
      by (subst Bochner_Integration.integral_add)
         (auto intro!: measure_pmf.integrable_const_bound[where B=1]
               simp: abs_mult mult_le_one U_bounded(2) measure_pmf.emeasure_eq_measure)
    finally have "measure (K w) UNIV - measure (K w) {x} = (v. U v x * indicator (- {x}) v K w)"
      by simp
    also have "measure (K w) UNIV - measure (K w) {x} = measure (K w) (UNIV - {x})"
      by (subst measure_pmf.finite_measure_Diff) auto
    finally have "0 = (v. indicator (- {x}) v K w) - (v. U v x * indicator (- {x}) v K w)"
      by (simp add: measure_pmf.emeasure_eq_measure Compl_eq_Diff_UNIV)
    also have " = (v. (1 - U v x) * indicator (- {x}) v K w)"
      by (subst Bochner_Integration.integral_diff[symmetric])
         (auto intro!: measure_pmf.integrable_const_bound[where B=1] Bochner_Integration.integral_cong
               simp: abs_mult mult_le_one U_bounded(2) split: split_indicator)
    also have "  (v. (1 - U y x) * indicator {y} v K w)" (is "_  ?rhs")
      using recurrent x
      by (intro integral_mono measure_pmf.integrable_const_bound[where B=1])
         (auto simp: abs_mult mult_le_one U_bounded(2) recurrent_iff_U_eq_1 field_simps
               split: split_indicator)
    also (xtrans) have "?rhs = (1 - U y x) * pmf (K w) y"
      by (simp add: measure_pmf.emeasure_eq_measure pmf.rep_eq)
    finally have "(1 - U y x) * pmf (K w) y = 0"
      by (auto intro!: antisym simp: U_bounded(2) mult_le_0_iff)
    with y  K w have "U y x = 1"
      by (simp add: set_pmf_iff)
    then have "U y x = 1" "H y x = 1"
      using H_eq(3)[of y x] H_eq(1)[of x] by (simp_all add: recurrent x)
    then have "(y, x)  acc"
      by (intro accI_U) auto
    with step have "(x, y)  communicating"
      by (auto simp add: communicating_def intro: rtrancl_trans)
    with recurrent x have "recurrent y"
      by (simp add: recurrent_iffI_communicating)
    note this U y x = 1 H y x = 1 (x, y)  communicating }
  note enabled = this

  from (x, y)  acc
  show "U y x = 1" "H y x = 1" "recurrent y" "(x, y)  communicating"
  proof induction
    case base then show "U x x = 1" "H x x = 1" "recurrent x" "(x, x)  communicating"
      using recurrent x H_eq(1)[of x] by (auto simp: recurrent_iff_U_eq_1 communicating_def)
  next
    case (step w y)
    with enabled[of w y] recurrent x H_eq(1)[of x]
    have "U y x = 1  H y x = 1  recurrent y  (x, y)  communicating"
      by (cases "x = y") (auto simp: recurrent_iff_U_eq_1 communicating_def)
    then show "U y x = 1" "H y x = 1" "recurrent y" "(x, y)  communicating"
      by auto
  qed
qed

lemma equiv_communicating: "equiv UNIV communicating"
  by (auto simp: equiv_def sym_def communicating_def refl_on_def trans_def)

lemma recurrent_class:
  assumes "recurrent x"
  shows "acc `` {x} = communicating `` {x}"
  using recurrent_acc(4)[OF recurrent x] by (auto simp: communicating_def)

lemma irreduccible_recurrent_class:
  assumes "recurrent x" shows "acc `` {x}  UNIV // communicating"
  unfolding recurrent_class[OF recurrent x] by (rule quotientI) simp

lemma essential_classI:
  assumes C: "C  UNIV // communicating"
  assumes eq: "x y. x  C  (x, y)  acc  y  C"
  shows "essential_class C"
  by (auto simp: essential_class_def intro: C) (metis eq)

lemma essential_recurrent_class:
  assumes "recurrent x" shows "essential_class (communicating `` {x})"
  unfolding recurrent_class[OF recurrent x, symmetric]
  apply (rule essential_classI)
  apply (rule irreduccible_recurrent_class[OF assms])
  apply (auto simp: communicating_def)
  done

lemma essential_classD2:
  "essential_class C  x  C  (x, y)  acc  y  C"
  unfolding essential_class_def by auto

lemma essential_classD3:
  "essential_class C  x  C  y  C  (x, y)  communicating"
  unfolding essential_class_def
  by (auto elim!: quotientE simp: communicating_def)

lemma AE_acc:
  shows "AE ω in T x. m. (x, (x ## ω) !! m)  acc"
  using AE_T_reachable
  by eventually_elim (auto simp: alw_HLD_iff_streams streams_iff_snth Stream_snth split: nat.splits)

lemma finite_essential_class_imp_recurrent:
  assumes C: "essential_class C" "finite C" and x: "x  C"
  shows "recurrent x"
proof -
  have "AE ω in T x. yC. alw (ev (HLD {y})) ω"
    using AE_T_reachable
  proof eventually_elim
    fix ω assume "alw (HLD (acc `` {x})) ω"
    then have "alw (HLD C) ω"
      by (rule alw_mono) (auto simp: HLD_iff intro: assms essential_classD2)
    then show "yC. alw (ev (HLD {y})) ω"
      by (rule pigeonhole_stream) fact
  qed
  then have "1 = 𝒫(ω in T x. yC. alw (ev (HLD {y})) ω)"
    by (subst (asm) T.prob_Collect_eq_1[symmetric]) (auto simp: finite C)
  also have " = measure (T x) (yC. {ωspace (T x). alw (ev (HLD {y})) ω})"
    by (intro arg_cong2[where f=measure]) auto
  also have "  (yC. H x y)"
    unfolding H_def using finite C by (rule T.finite_measure_subadditive_finite) auto
  also have " = (yC. U x y * H y y)"
    by (auto intro!: sum.cong H_eq)
  finally have "yC. recurrent y"
    by (rule_tac ccontr) (simp add: H_eq(2))
  then obtain y where "y  C" "recurrent y" ..
  from essential_classD3[OF C(1) x this(1)] recurrent_acc(3)[OF this(2)]
  show "recurrent x"
    by (simp add: communicating_def)
qed

lemma irreducibleD:
  "C  UNIV // communicating  a  C  b  C  (a, b)  communicating"
  by (auto elim!: quotientE simp: communicating_def)

lemma irreducibleD2:
  "C  UNIV // communicating  a  C  (a, b)  communicating  b  C"
  by (auto elim!: quotientE simp: communicating_def)

lemma essential_class_iff_recurrent:
  "finite C  C  UNIV // communicating  essential_class C  (xC. recurrent x)"
  by (metis finite_essential_class_imp_recurrent irreducibleD2 recurrent_acc(4) essential_classI)

definition "U' x y = (+ω. eSuc (sfirst (HLD {y}) ω) T x)"

lemma U'_neq_zero[simp]: "U' x y  0"
  unfolding U'_def by (simp add: nn_integral_add)

definition "gf_U' x y z = (n. u x y n * Suc n * z ^ n)"

definition "pos_recurrent x  recurrent x  U' x x  "

lemma summable_gf_U':
  assumes z: "norm z < 1"
  shows "summable (λn. u x y n * Suc n * z ^ n)"
proof -
  have "summable (λn. n * ¦z¦ ^ n)"
  proof (rule root_test_convergence)
    have "(λn. root n n * ¦z¦)  1 * ¦z¦"
      by (intro tendsto_intros LIMSEQ_root)
    then show "(λn. root n (norm (n * ¦z¦ ^ n)))  ¦z¦"
      by (rule filterlim_cong[THEN iffD1, rotated 3])
         (auto intro!: exI[of _ 1]
               simp add: abs_mult u_nonneg real_root_mult power_abs eventually_sequentially real_root_power)
  qed (insert z, simp add: abs_less_iff)
  note summable_mult[OF this, of "1 / ¦z¦"]
  from summable_ignore_initial_segment[OF this, of 1]
  show "summable (λn. u x y n * Suc n * z ^ n)"
    apply (rule summable_comparison_test[rotated])
    using z
    apply (auto intro!: exI[of _ 1]
                simp: abs_mult u_nonneg power_abs Suc_le_eq gr0_conv_Suc field_simps le_divide_eq u_le_1
                simp del: of_nat_Suc)
    done
qed

lemma gf_U'_nonneg[simp]: "0 < z  z < 1  0  gf_U' x y z"
  unfolding gf_U'_def
  by (intro suminf_nonneg summable_gf_U') (auto simp: u_nonneg)

lemma DERIV_gf_U:
  fixes z :: real assumes z: "0 < z" "z < 1"
  shows "DERIV (gf_U x y) z :> gf_U' x y z"
  unfolding gf_U_def[abs_def]  gf_U'_def real_scaleR_def u_def[symmetric]
  using z by (intro DERIV_power_series'[where R=1] summable_gf_U') auto

lemma sfirst_finiteI_recurrent:
  "recurrent x  (x, y)  acc  AE ω in T x. sfirst (HLD {y}) ω < "
  using recurrent_acc(1)[of y x] recurrent_acc[of x y]
    T.AE_prob_1[of x "{ωspace (T x). ev (HLD {y}) ω}"]
  unfolding sfirst_finite U_def by (simp add: space_stream_space communicating_def)

lemma U'_eq_suminf:
  assumes x: "recurrent x" "(x, y)  acc"
  shows "U' x y = (i. ennreal (u x y i * Suc i))"
proof -
  have "(+ω. eSuc (sfirst (HLD {y}) ω) T x) =
      (+ω. (i. ennreal (Suc i) * indicator {ωspace (T y). ev_at (HLD {y}) i ω} ω) T x)"
    using sfirst_finiteI_recurrent[OF x]
  proof (intro nn_integral_cong_AE, eventually_elim)
    fix ω assume "sfirst (HLD {y}) ω < "
    then obtain n :: nat where [simp]: "sfirst (HLD {y}) ω = n"
      by auto
    show "eSuc (sfirst (HLD {y}) ω) = (i. ennreal (Suc i) * indicator {ωspace (T y). ev_at (HLD {y}) i ω} ω)"
      by (subst suminf_cmult_indicator[where i=n])
         (auto simp: disjoint_family_on_def ev_at_unique space_stream_space
                     sfirst_eq_enat_iff[symmetric] ennreal_of_nat_eq_real_of_nat
               split: split_indicator)
  qed
  also have " = (i. ennreal (Suc i) * emeasure (T x) {ωspace (T x). ev_at (HLD {y}) i ω})"
    by (subst nn_integral_suminf)
       (auto intro!: arg_cong[where f=suminf] nn_integral_cmult_indicator simp: fun_eq_iff)
  finally show ?thesis
    by (simp add: U'_def u_def T.emeasure_eq_measure mult_ac ennreal_mult)
qed

lemma gf_U'_tendsto_U':
  assumes x: "recurrent x" "(x, y)  acc"
  shows "((λz. ennreal (gf_U' x y z))  U' x y) (at_left 1)"
  unfolding U'_eq_suminf[OF x] gf_U'_def
  by (auto intro!: power_series_tendsto_at_left summable_gf_U' mult_nonneg_nonneg u_nonneg simp del: of_nat_Suc)

lemma one_le_integral_t:
  assumes x: "recurrent x" shows "1  U' x x"
  by (simp add: nn_integral_add T.emeasure_space_1 U'_def del: space_T)

lemma gf_U'_pos:
  fixes z :: real
  assumes z: "0 < z" "z < 1" and "U x y  0"
  shows "0 < gf_U' x y z"
  unfolding gf_U'_def
proof (subst suminf_pos_iff)
  show "summable (λn. u x y n * real (Suc n) * z ^ n)"
    using z by (intro summable_gf_U') simp
  show pos: "n. 0  u x y n * real (Suc n) * z ^ n"
    using u_nonneg z by auto
  show "n. 0 < u x y n * real (Suc n) * z ^ n"
  proof (rule ccontr)
    assume "¬ (n. 0 < u x y n * real (Suc n) * z ^ n)"
    with pos have "n. u x y n * real (Suc n) * z ^ n = 0"
      by (intro antisym allI) (simp_all add: not_less)
    with z have "u x y = (λn. 0)"
      by (intro ext) simp
    with u_sums_U[of x y, THEN sums_unique] U x y  0 show False
      by simp
  qed
qed

lemma inverse_gf_U'_tendsto:
  assumes "recurrent y"
  shows "((λx. - 1 / - gf_U' y y x)  enn2real (1 / U' y y)) (at_left (1::real))"
proof cases
  assume inf: "U' y y = "
  with gf_U'_tendsto_U'[of y y] recurrent y
  have "LIM z (at_left 1). gf_U' y y z :> at_top"
    by (auto simp: ennreal_tendsto_top_eq_at_top U'_def)
  then have "LIM z (at_left 1). gf_U' y y z :> at_infinity"
    by (rule filterlim_mono) (auto simp: at_top_le_at_infinity)
  with inf show ?thesis
    by (auto intro!: tendsto_divide_0)
next
  assume fin: "U' y y  "
  then obtain r where r: "U' y y = ennreal r" and [simp]: "0  r"
    by (cases "U' y y") (auto simp: U'_def)
  then have eq: "enn2real (1 / U' y y) = - 1 / - r" and "1  r"
    using one_le_integral_t[OF recurrent y]
    by (auto simp add: ennreal_1[symmetric] divide_ennreal simp del: ennreal_1)
  have "((λz. ennreal (gf_U' y y z))  ennreal r) (at_left 1)"
    using gf_U'_tendsto_U'[OF recurrent y, of y] r by simp
  then have gf_U': "(gf_U' y y  r) (at_left (1::real))"
    by (rule tendsto_ennrealD)
       (insert summable_gf_U', auto intro!: eventually_at_left_1 suminf_nonneg simp: gf_U'_def u_nonneg)
  show ?thesis
    using 1  r unfolding eq by (intro tendsto_intros gf_U') simp
qed

lemma gf_G_pos:
  fixes z :: real
  assumes z: "0 < z" "z < 1" and *: "(x, y)  acc"
  shows "0 < gf_G x y z"
  unfolding gf_G_def
proof (subst suminf_pos_iff)
  show "summable (λn. p x y n *R z ^ n)"
    using z by (intro convergence_G convergence_G_less_1) simp
  show pos: "n. 0  p x y n *R z ^ n"
    using z by (auto intro!: mult_nonneg_nonneg p_nonneg)
  show "n. 0 < p x y n *R z ^ n"
  proof (rule ccontr)
    assume "¬ (n. 0 < p x y n *R z ^ n)"
    with pos have "n. p x y n * z ^ n = 0"
      by (intro antisym allI) (simp_all add: not_less)
    with z have "n. p x y n = 0"
      by simp
    with *[THEN accD_pos] show False
      by simp
  qed
qed

lemma pos_recurrentI_communicating:
  assumes y: "pos_recurrent y" and x: "(y, x)  communicating"
  shows "pos_recurrent x"
proof -
  from y x have recurrent: "recurrent y" "recurrent x" and fin: "U' y y  "
    by (auto simp: pos_recurrent_def recurrent_iffI_communicating nn_integral_add)
  have pos: "0 < enn2real (1 / U' y y)"
    using one_le_integral_t[OF recurrent y] fin
    by (auto simp: U'_def enn2real_positive_iff less_top[symmetric] ennreal_zero_less_divide ennreal_divide_eq_top_iff)

  from fin obtain r where r: "U' y y = ennreal r" and [simp]: "0  r"
    by (cases "U' y y") (auto simp: U'_def)

  from x obtain n m where "0 < p x y n" "0 < p y x m"
    by (auto dest!: accD_pos simp: communicating_def)

  let ?L = "at_left (1::real)"
  have le: "eventually (λz. p x y n * p y x m * z^(n + m)  (1 - gf_U y y z) / (1 - gf_U x x z)) ?L"
  proof (rule eventually_at_left_1)
    fix z :: real assume z: "0 < z" "z < 1"
    then have conv: "x. convergence_G x x z"
      by (intro convergence_G_less_1) simp
    have sums: "(λi. (p x y n * p y x m * z^(n + m)) * (p y y i * z^i)) sums ((p x y n * p y x m * z^(n + m)) * gf_G y y z)"
      unfolding gf_G_def
      by (intro sums_mult summable_sums) (auto intro!: conv convergence_G[where 'a=real, simplified])
    have "(i. (p x y n * p y x m * z^(n + m)) * (p y y i * z^i))  (i. p x x (i + (n + m)) * z^(i + (n + m)))"
    proof (intro allI suminf_le sums_summable[OF sums] summable_ignore_initial_segment convergence_G[where 'a=real, simplified] convergence_G_less_1)
      show "norm z < 1" using z by simp
      fix i
      have "(p x y n * p y y ((n + i) - n)) * p y x ((n + i + m) - (n + i))  p x y (n + i) * p y x ((n + i + m) - (n + i))"
        by (intro mult_right_mono prob_reachable_le) simp_all
      also have "  p x x (n + i + m)"
         by (intro prob_reachable_le) simp_all
      finally show "p x y n * p y x m * z ^ (n + m) * (p y y i * z ^ i)  p x x (i + (n + m)) * z ^ (i + (n + m))"
        using z by (auto simp add: ac_simps power_add intro!: mult_left_mono)
    qed
    also have "  gf_G x x z"
      unfolding gf_G_def
      using z
      apply (subst (2) suminf_split_initial_segment[where k="n + m"])
      apply (intro convergence_G conv)
      apply (simp add: sum_nonneg)
      done
    finally have "(p x y n * p y x m * z^(n + m)) * gf_G y y z  gf_G x x z"
      using sums_unique[OF sums] by simp
    then have "(p x y n * p y x m * z^(n + m))  gf_G x x z / gf_G y y z"
      using z gf_G_pos[of z y y] by (simp add: field_simps)
    also have " = (1 - gf_U y y z) / (1 - gf_U x x z)"
      unfolding gf_G_eq_gf_U[OF conv] using gf_G_eq_gf_U(2)[OF conv] by (simp add: field_simps )
    finally show "p x y n * p y x m * z^(n + m)  (1 - gf_U y y z) / (1 - gf_U x x z)" .
  qed

  have "U' x x  "
  proof
    assume "U' x x = "
    have "((λz. (1 - gf_U y y z) / (1 - gf_U x x z))  0) ?L"
    proof (rule lhopital_left)
      show "((λz. 1 - gf_U y y z)  0) ?L"
        using gf_U[of y] recurrent_iff_U_eq_1[of y] recurrent y by (auto intro!: tendsto_eq_intros)
      show "((λz. 1 - gf_U x x z)  0) ?L"
        using gf_U[of x] recurrent_iff_U_eq_1[of x] recurrent x by (auto intro!: tendsto_eq_intros)
      show "eventually (λz. 1 - gf_U x x z  0) ?L"
        by (auto intro!: eventually_at_left_1 simp: gf_G_eq_gf_U(2) convergence_G_less_1)
      show "eventually (λz. - gf_U' x x z  0) ?L"
        using gf_U'_pos[of _ x x] recurrent_iff_U_eq_1[of x] recurrent x
        by (auto intro!: eventually_at_left_1) (metis less_le)
      show "eventually (λz. DERIV (λxa. 1 - gf_U x x xa) z :> - gf_U' x x z) ?L"
        by (auto intro!: eventually_at_left_1 derivative_eq_intros DERIV_gf_U)
      show "eventually (λz. DERIV (λxa. 1 - gf_U y y xa) z :> - gf_U' y y z) ?L"
        by (auto intro!: eventually_at_left_1 derivative_eq_intros DERIV_gf_U)

      have "(gf_U' y y  U' y y) ?L"
        using recurrent y by (rule gf_U'_tendsto_U') simp
      then have *: "(gf_U' y y  r) ?L"
        by (auto simp add: r eventually_at_left_1 dest!: tendsto_ennrealD)
      moreover
      have "(gf_U' x x  U' x x) ?L"
        using recurrent x by (rule gf_U'_tendsto_U') simp
      then have "LIM z ?L. - gf_U' x x z :> at_bot"
        by (simp add: ennreal_tendsto_top_eq_at_top U' x x =  filterlim_uminus_at_top
                 del: ennreal_of_enat_eSuc)
      then have "LIM z ?L. - gf_U' x x z :> at_infinity"
        by (rule filterlim_mono) (auto simp: at_bot_le_at_infinity)
      ultimately show "((λz. - gf_U' y y z / - gf_U' x x z)  0) ?L"
        by (intro tendsto_divide_0[where c="- r"] tendsto_intros)
    qed
    moreover
    have "((λz. p x y n * p y x m * z^(n + m))  p x y n * p y x m) ?L"
      by (auto intro!: tendsto_eq_intros)
    ultimately have "p x y n * p y x m  0"
      using le by (rule tendsto_le[OF trivial_limit_at_left_real])
    with 0 < p x y n 0 < p y x m show False
      by (auto simp add: mult_le_0_iff)
  qed
  with recurrent x show ?thesis
    by (simp add: pos_recurrent_def nn_integral_add)
qed

lemma pos_recurrent_iffI_communicating:
  "(y, x)  communicating  pos_recurrent y  pos_recurrent x"
  using pos_recurrentI_communicating[of x y] pos_recurrentI_communicating[of y x]
  by (auto simp add: communicating_def)

lemma U_le_F: "U x y  F x y"
  by (auto simp: U_def F_def intro!: T.finite_measure_mono)

lemma not_empty_irreducible: "C  UNIV // communicating  C  {}"
  by (auto simp: quotient_def Image_def communicating_def)

subsection ‹Stationary distribution›

definition stat :: "'s set  's measure" where
  "stat C = point_measure UNIV (λx. indicator C x / U' x x)"

lemma sets_stat[simp]: "sets (stat C) = sets (count_space UNIV)"
  by (simp add: stat_def sets_point_measure)

lemma space_stat[simp]: "space (stat C) = UNIV"
  by (simp add: stat_def space_point_measure)

lemma stat_subprob:
  assumes C: "essential_class C" and "countable C" and pos: "cC. pos_recurrent c"
  shows "emeasure (stat C) C  1"
proof -
  let ?L = "at_left (1::real)"
  from finite_sequence_to_countable_set[OF countable C]
  obtain A where A: "i. A i  C" "i. A i  A (Suc i)" "i. finite (A i)" " (range A) = C"
    by blast
  then have "(λn. emeasure (stat C) (A n))  emeasure (stat C) (i. A i)"
    by (intro Lim_emeasure_incseq) (auto simp: incseq_Suc_iff)
  then have "emeasure (stat C) (i. A i)  1"
  proof (rule LIMSEQ_le[OF _ tendsto_const], intro exI allI impI)
    fix n
    from A(1,3) have A_n: "finite (A n)"
      by auto

    from C have "C  {}"
      by (simp add: essential_class_def not_empty_irreducible)
    then obtain x where "x  C" by auto

    have "((λz. (yA n. gf_F x y z * ((1 - z) / (1 - gf_U y y z))))  (yA n. F x y * enn2real (1 / U' y y))) ?L"
    proof (intro tendsto_intros gf_F, rule lhopital_left)
      fix y assume "y  A n"
      with A n  C have "y  C"
        by auto
      show "((-) 1  0) ?L"
        by (intro tendsto_eq_intros) simp_all
      have "recurrent y"
        using pos[THEN bspec, OF yC] by (simp add: pos_recurrent_def)
      then have "U y y = 1"
        by (simp add: recurrent_iff_U_eq_1)

      show "((λx. 1 - gf_U y y x)  0) ?L"
        using gf_U[of y y] U y y = 1 by (intro tendsto_eq_intros) auto
      show "eventually (λx. 1 - gf_U y y x  0) ?L"
        using gf_G_eq_gf_U(2)[OF convergence_G_less_1, where 'z=real] by (auto intro!: eventually_at_left_1)
      have "eventually (λx. 0 < gf_U' y y x) ?L"
        by (intro eventually_at_left_1 gf_U'_pos) (simp_all add: U y y = 1)
      then show "eventually (λx. - gf_U' y y x  0) ?L"
        by eventually_elim simp
      show "eventually (λx. DERIV (λx. 1 - gf_U y y x) x :> - gf_U' y y x) ?L"
        by (auto intro!: eventually_at_left_1 derivative_eq_intros DERIV_gf_U)
      show "eventually (λx. DERIV ((-) 1) x :> - 1) ?L"
        by (auto intro!: eventually_at_left_1 derivative_eq_intros)
      show "((λx. - 1 / - gf_U' y y x)  enn2real (1 / U' y y)) ?L"
        using recurrent y by (rule inverse_gf_U'_tendsto)
    qed
    also have "(yA n. F x y * enn2real (1 / U' y y)) = (yA n. enn2real (1 / U' y y))"
    proof (intro sum.cong refl)
      fix y assume "y  A n"
      with A n  C have "y  C" by auto
      with x  C have "(x, y)  communicating"
        by (rule essential_classD3[OF C])
      with yC have "recurrent y" "(y, x)  acc"
        using pos[THEN bspec, of y] by (auto simp add: pos_recurrent_def communicating_def)
      then have "U x y = 1"
        by (rule recurrent_acc)
      with F_le_1[of x y] U_le_F[of x y] have "F x y = 1" by simp
      then show "F x y * enn2real (1 / U' y y) = enn2real (1 / U' y y)"
        by simp
    qed
    finally have le: "(yA n. enn2real (1 / U' y y))  1"
    proof (rule tendsto_le[OF trivial_limit_at_left_real tendsto_const], intro eventually_at_left_1)
      fix z :: real assume z: "0 < z" "z < 1"
      with x  C have "norm z < 1"
        by auto
      then have conv: "x y. convergence_G x y z"
        by (simp add: convergence_G_less_1)
      have "(yA n. gf_F x y z / (1 - gf_U y y z)) = (yA n. gf_G x y z)"
        using norm z < 1
        apply (intro sum.cong refl)
        apply (subst gf_G_eq_gf_F)
        apply assumption
        apply (subst gf_G_eq_gf_U(1)[OF conv])
        apply auto
        done
      also have " = (yA n. n. p x y n * z^n)"
        by (simp add: gf_G_def)
      also have "  = (i. yA n. p x y i *R z^i)"
        by (subst suminf_sum[OF convergence_G[OF conv]]) simp
      also have "   (i. z^i)"
      proof (intro suminf_le summable_sum convergence_G conv summable_geometric allI)
        fix l
        have "(yA n. p x y l *R z ^ l) = (yA n. p x y l) * z ^ l"
          by (simp add: sum_distrib_right)
        also have "  z ^ l"
        proof (intro mult_left_le_one_le)
          have "(yA n. p x y l) = 𝒫(ω in T x. (x ## ω) !! l  A n)"
            unfolding p_def using finite (A n)
            by (subst T.finite_measure_finite_Union[symmetric])
               (auto simp: disjoint_family_on_def intro!: arg_cong2[where f=measure])
          then show "(yA n. p x y l)  1"
            by simp
        qed (insert z, auto simp: sum_nonneg)
        finally show "(yA n. p x y l *R z ^ l)  z ^ l" .
      qed fact
      also have " = 1 / (1 - z)"
        using sums_unique[OF geometric_sums, OF norm z < 1] ..
      finally have "(yA n. gf_F x y z / (1 - gf_U y y z))  1 / (1 - z)" .
      then have "(yA n. gf_F x y z / (1 - gf_U y y z)) * (1 - z)  1"
        using z by (simp add: field_simps)
      then have "(yA n. gf_F x y z / (1 - gf_U y y z) * (1 - z))  1"
        by (simp add: sum_distrib_right)
      then show "(yA n. gf_F x y z * ((1 - z) / (1 - gf_U y y z)))  1"
        by simp
    qed

    from A_n have "emeasure (stat C) (A n) = (yA n. emeasure (stat C) {y})"
      by (intro emeasure_eq_sum_singleton) simp_all
    also have " = (yA n. inverse (U' y y))"
      unfolding stat_def U'_def using A(1)[of n]
      apply (intro sum.cong refl)
      apply (subst emeasure_point_measure_finite2)
        apply (auto simp: divide_ennreal_def Collect_conv_if)
      done
    also have " = ennreal (yA n. enn2real (1 / U' y y))"
      apply (subst sum_ennreal[symmetric], simp)
    proof (intro sum.cong refl)
      fix y assume "y  A n"
      with A n  C pos have "pos_recurrent y"
        by auto
      with one_le_integral_t[of y] obtain r where "U' y y = ennreal r" "1  U' y y" and [simp]: "0  r"
        by (cases "U' y y") (auto simp: pos_recurrent_def nn_integral_add)
      then show "inverse (U' y y) = ennreal (enn2real (1 / U' y y))"
        by (simp add: ennreal_1[symmetric] divide_ennreal inverse_ennreal inverse_eq_divide del: ennreal_1)
    qed
    also have "  1"
      using le by simp
    finally show "emeasure (stat C) (A n)  1" .
  qed
  with A show ?thesis
    by simp
qed

lemma emeasure_stat_not_C:
  assumes "y  C"
  shows "emeasure (stat C) {y} = 0"
  unfolding stat_def using y  C
  by (subst emeasure_point_measure_finite2) auto

definition stationary_distribution :: "'s pmf  bool" where
  "stationary_distribution N  N = bind_pmf N K"

lemma stationary_distributionI:
  assumes le: "y. (x. pmf (K x) y measure_pmf N)  pmf N y"
  shows "stationary_distribution N"
  unfolding stationary_distribution_def
proof (rule pmf_eqI antisym)+
  fix i
  show "pmf (bind_pmf N K) i  pmf N i"
    by (simp add: pmf_bind le)

  define Ω where "Ω = N  (iN. set_pmf (K i))"
  then have Ω: "countable Ω"
    by (auto intro: countable_set_pmf)
  then interpret N: sigma_finite_measure "count_space Ω"
    by (rule sigma_finite_measure_count_space_countable)
  interpret pN: pair_sigma_finite N "count_space Ω"
    by unfold_locales

  have measurable_pmf[measurable]: "(λ(x, y). pmf (K x) y)  borel_measurable (N M count_space Ω)"
    unfolding measurable_split_conv
    apply (rule measurable_compose_countable'[OF _ measurable_snd])
    apply (rule measurable_compose[OF measurable_fst])
    apply (simp_all add: Ω)
    done

  { assume *: "(y. pmf (K y) i N) < pmf N i"
    have "0  (y. pmf (K y) i N)"
      by (intro integral_nonneg_AE) simp
    with * have i: "i  set_pmf N" "i  Ω"
      by (auto simp: set_pmf_iff Ω_def not_le[symmetric])
    from * have "0 < pmf N i - (y. pmf (K y) i N)"
      by (simp add: field_simps)
    also have " = (t. (pmf N i - (y. pmf (K y) i N)) * indicator {i} t count_space Ω)"
      by (simp add: i)
    also have "  (t. pmf N t - y. pmf (K y) t N count_space Ω)"
      using le
      by (intro integral_mono integrable_diff)
         (auto simp: i pmf_bind[symmetric] integrable_pmf field_simps split: split_indicator)
    also have " = (t. pmf N t count_space Ω) - (t. y. pmf (K y) t N count_space Ω)"
      by (subst Bochner_Integration.integral_diff) (auto intro!: integrable_pmf simp: pmf_bind[symmetric])
    also have "(t. y. pmf (K y) t N count_space Ω) = (y. t. pmf (K y) t count_space Ω N)"
      apply (intro pN.Fubini_integral integrable_iff_bounded[THEN iffD2] conjI)
      apply (auto simp add: N.nn_integral_fst[symmetric] nn_integral_eq_integral integrable_pmf)
      unfolding less_top[symmetric] unfolding infinity_ennreal_def[symmetric]
      apply (intro integrableD)
      apply (auto intro!: measure_pmf.integrable_const_bound[where B=1]
                  simp: AE_measure_pmf_iff integral_nonneg_AE integral_pmf)
      done
    also have "(y. t. pmf (K y) t count_space Ω N) = (y. 1 N)"
      by (intro integral_cong_AE)
         (auto simp: AE_measure_pmf_iff integral_pmf Ω_def intro!: measure_pmf.prob_eq_1[THEN iffD2])
    finally have False
      using measure_pmf.prob_space[of N] by (simp add: integral_pmf field_simps not_le[symmetric]) }
  then show "pmf N i  pmf (bind_pmf N K) i"
    by (auto simp: pmf_bind not_le[symmetric])
qed

lemma stationary_distribution_iterate:
  assumes N: "stationary_distribution N"
  shows "ennreal (pmf N y) = (+x. p x y n N)"
proof (induct n arbitrary: y)
  have [simp]: "x y. ennreal (if x = y then 1 else 0) = indicator {y} x"
    by simp
  case 0 then show ?case
    by (simp add: p_0 pmf.rep_eq measure_pmf.emeasure_eq_measure)
next
  case (Suc n) with N show ?case
    apply (simp add: nn_integral_eq_integral[symmetric] p_le_1 p_Suc'
                     measure_pmf.integrable_const_bound[where B=1])
    apply (subst nn_integral_bind[symmetric, where B="count_space UNIV"])
    apply (auto simp: stationary_distribution_def measure_pmf_bind[symmetric]
                simp del: measurable_pmf_measure1)
    done
qed

lemma stationary_distribution_iterate':
  assumes "stationary_distribution N"
  shows "measure N {y} = (x. p x y n N)"
  using stationary_distribution_iterate[OF assms]
  by (subst (asm) nn_integral_eq_integral)
     (auto intro!: measure_pmf.integrable_const_bound[where B=1] simp: p_le_1 pmf.rep_eq)

lemma stationary_distributionD:
  assumes C: "essential_class C" "countable C"
  assumes N: "stationary_distribution N" "N  C"
  shows "xC. pos_recurrent x" "measure_pmf N = stat C"
proof -
  have integrable_K: "f x. integrable N (λs. pmf (K s) (f x))"
    by (rule measure_pmf.integrable_const_bound[where B=1]) (simp_all add: pmf_le_1)

  have measure_C: "measure N C = 1" and ae_C: "AE x in N. x  C"
    using N C measure_pmf.prob_eq_1[of C] by (auto simp: AE_measure_pmf_iff)

  have integrable_p: "n y. integrable N (λx. p x y n)"
    by (rule measure_pmf.integrable_const_bound[where B=1]) (simp_all add: p_le_1)

  { fix e :: real assume "0 < e"
    then have [simp]: "0  e" by simp
    have "AC. finite A  1 - e < measure N A"
    proof (rule ccontr)
      assume contr: "¬ (A  C. finite A  1 - e < measure N A)"
      from finite_sequence_to_countable_set[OF countable C]
      obtain F where F: "i. F i  C" "i. F i  F (Suc i)" "i. finite (F i)" " (range F) = C"
        by blast
      then have *: "(λn. measure N (F n))  measure N (i. F i)"
        by (intro measure_pmf.finite_Lim_measure_incseq) (auto simp: incseq_Suc_iff)
      with F contr have "measure N (i. F i)  1 - e"
        by (intro LIMSEQ_le[OF * tendsto_const]) (auto simp: not_less)
      with F 0 < e show False
        by (simp add: measure_C)
    qed
    then obtain A where "A  C" "finite A" and e: "1 - e < measure N A" by auto

    { fix y n assume "y  C"
      from N(1) have "measure N {y} = (x. p x y n N)"
        by (rule stationary_distribution_iterate')
      also have "  (x. p x y n * indicator A x + indicator (C - A) x N)"
        using ae_C A  C
        by (intro integral_mono_AE)
           (auto elim!: eventually_mono
                 intro!: integral_add integral_indicator p_le_1 integrable_real_mult_indicator
                   integrable_add
                 split: split_indicator simp: integrable_p less_top[symmetric] top_unique)
      also have " = (x. p x y n * indicator A x N) + measure N (C - A)"
        using ae_C A  C
        apply (subst Bochner_Integration.integral_add)
        apply (auto elim!: eventually_mono
                    intro!: integral_add integral_indicator p_le_1 integrable_real_mult_indicator
                    split: split_indicator simp: integrable_p less_top[symmetric] top_unique)
        done
      also have "  (x. p x y n * indicator A x N) + e"
        using e A  C  by (simp add: measure_pmf.finite_measure_Diff measure_C)
      finally have "measure N {y}  (x. p x y n * indicator A x N) + e" .
      then have "emeasure N {y}  ennreal (x. p x y n * indicator A x N) + e"
        by (simp add: measure_pmf.emeasure_eq_measure ennreal_plus[symmetric] del: ennreal_plus)
      also have " = (+x. ennreal (p x y n) * indicator A x N) + e"
        by (subst nn_integral_eq_integral[symmetric])
           (auto intro!: measure_pmf.integrable_const_bound[where B=1]
                 simp: abs_mult p_le_1 mult_le_one ennreal_indicator ennreal_mult)
      finally have "emeasure N {y}  (+x. ennreal (p x y n) * indicator A x N) + e" . }
    note v_le = this

    { fix y and z :: real assume y: "y  C" and z: "0 < z" "z < 1"
      have summable_int_p: "summable (λn. ( x. p x y n * indicator A x N) * (1 - z) * z ^ n)"
        using yC z A  C
        by (auto intro!: summable_comparison_test[OF _ summable_mult[OF summable_geometric[of z], of 1]] exI[of _ 0] mult_le_one
                            measure_pmf.integral_le_const integrable_real_mult_indicator integrable_p AE_I2 p_le_1
                    simp: abs_mult integral_nonneg_AE)

      from y z have sums_y: "(λn. measure N {y} * (1 - z) * z ^ n) sums measure N {y}"
        using sums_mult[OF geometric_sums[of z], of "measure N {y} * (1 - z)"] by simp
      then have "emeasure N {y} = ennreal (n. (measure N {y} * (1 - z)) * z ^ n)"
        by (auto simp add: sums_unique[symmetric] measure_pmf.emeasure_eq_measure)
      also have " = (n. emeasure N {y} * (1 - z) * z ^ n)"
        using z  summable_mult[OF summable_geometric[of z], of "measure_pmf.prob N {y} * (1 - z)"]
        by (subst suminf_ennreal[symmetric])
           (auto simp: measure_pmf.emeasure_eq_measure ennreal_mult[symmetric] ennreal_suminf_neq_top)
      also have "  (n. ((+x. ennreal (p x y n) * indicator A x N) + e) * (1 - z) * z ^ n)"
        using yC z A  C
        by (intro suminf_le mult_right_mono v_le allI)
           (auto simp: measure_pmf.emeasure_eq_measure)
      also have " = (n. (+x. ennreal (p x y n) * indicator A x N) * (1 - z) * z ^ n) + e"
        using 0 < e z sums_mult[OF geometric_sums[of z], of "e * (1 - z)"] 0<z z<1
        by (simp add: distrib_right suminf_add[symmetric] ennreal_suminf_cmult[symmetric]
                      ennreal_mult[symmetric] suminf_ennreal_eq sums_unique[symmetric]
                 del: ennreal_suminf_cmult)
      also have " = (n. ennreal (1 - z) * ((+x. ennreal (p x y n) * indicator A x N) * z ^ n)) + e"
        by (simp add: ac_simps)
      also have " = ennreal (1 - z) * (n. ((+x. ennreal (p x y n) * indicator A x N) * z ^ n)) + e"
        using z by (subst ennreal_suminf_cmult) simp_all
      also have "(n. ((+x. ennreal (p x y n) * indicator A x N) * z ^ n)) =
          (n. (+x. ennreal (p x y n * z ^ n) * indicator A x N))"
        using z by (simp add: ac_simps nn_integral_cmult[symmetric] ennreal_mult)
      also have " = (+x. ennreal (gf_G x y z) * indicator A x N)"
        using z
        apply (subst nn_integral_suminf[symmetric])
        apply (auto simp add: gf_G_def simp del: suminf_ennreal
                    intro!: ennreal_mult_right_cong suminf_ennreal2 nn_integral_cong)
        apply (intro summable_comparison_test[OF _ summable_mult[OF summable_geometric[of z], of 1]] impI)
        apply (simp_all add: abs_mult p_le_1 mult_le_one power_le_one split: split_indicator)
        done
      also have " = (+x. ennreal (gf_F x y z * gf_G y y z) * indicator A x N)"
        using z by (intro nn_integral_cong) (simp add: gf_G_eq_gf_F[symmetric])
      also have " = ennreal (gf_G y y z) * (+x. ennreal (gf_F x y z) * indicator A x N)"
        using z by (subst nn_integral_cmult[symmetric]) (simp_all add: gf_G_nonneg gf_F_nonneg ac_simps ennreal_mult)
      also have " = ennreal (1 / (1 - gf_U y y z)) * (+x. ennreal (gf_F x y z) * indicator A x N)"
        using z y  C by (subst gf_G_eq_gf_U) (auto intro!: convergence_G_less_1)
      finally have "emeasure N {y}  ennreal ((1 - z) / (1 - gf_U y y z)) * (+x. gf_F x y z * indicator A x N) + e"
        using z
        by (subst (asm) mult.assoc[symmetric])
           (simp add: ennreal_indicator[symmetric] ennreal_mult'[symmetric] gf_F_nonneg)
      then have "measure N {y}  (1 - z) / (1 - gf_U y y z) * (x. gf_F x y z * indicator A x N) + e"
        using z
        by (subst (asm) nn_integral_eq_integral[OF measure_pmf.integrable_const_bound[where B=1]])
           (auto simp: gf_F_nonneg gf_U_le_1 gf_F_le_1 measure_pmf.emeasure_eq_measure mult_le_one
                       ennreal_mult''[symmetric] ennreal_plus[symmetric]
                 simp del: ennreal_plus) }
    then have "A  C. finite A  (yC. z. 0 < z  z < 1  measure N {y}  (1 - z) / (1 - gf_U y y z) * (x. gf_F x y z * indicator A x N) + e)"
      using A  C finite A by auto }
  note eps = this

  { fix y A assume "y  C" "finite A" "A  C"
    then have "((λz. x. gf_F x y z * indicator A x N)  x. F x y * indicator A x N) (at_left 1)"
      by (subst (1 2) integral_measure_pmf[of A]) (auto intro!: tendsto_intros gf_F simp: indicator_eq_0_iff) }
  note int_gf_F = this

  have all_recurrent: "yC. recurrent y"
  proof (rule ccontr)
    assume "¬ (yC. recurrent y)"
    then obtain x where "x  C" "¬ recurrent x" by auto
    then have transient: "x. x  C  ¬ recurrent x"
      using C by (auto simp: essential_class_def recurrent_iffI_communicating[symmetric] elim!: quotientE)

    { fix y assume "y  C"
      with transient have "U y y < 1"
        by (metis recurrent_iff_U_eq_1 U_cases)
      have "measure N {y}  0"
      proof (rule dense_ge)
        fix e :: real assume "0 < e"
        from eps[OF this] y  C obtain A where
          A: "finite A" "A  C" and
          le: "z. 0 < z  z < 1  measure N {y}  (1 - z) / (1 - gf_U y y z) * (x. gf_F x y z * indicator A x N) + e"
          by auto
        have "((λz. (1 - z) / (1 - gf_U y y z) * (x. gf_F x y z * indicator A x N) + e) 
          (1 - 1) / (1 - U y y) * (x. F x y * indicator A x N) + e) (at_left (1::real))"
          using A U y y < 1 y  C by (intro tendsto_intros gf_U int_gf_F) auto
        then have 1: "((λz. (1 - z) / (1 - gf_U y y z) * (x. gf_F x y z * indicator A x N) + e)  e) (at_left (1::real))"
          by simp
        with le show "measure N {y}  e"
          by (intro tendsto_le[OF trivial_limit_at_left_real _ tendsto_const])
             (auto simp: eventually_at_left_1)
      qed
      then have "measure N {y} = 0"
        by (intro antisym measure_nonneg) }
    then have "emeasure N C = 0"
      by (subst emeasure_countable_singleton) (auto simp: measure_pmf.emeasure_eq_measure nn_integral_0_iff_AE ae_C C)
    then show False
      using measure N C = 1 by (simp add: measure_pmf.emeasure_eq_measure)
  qed
  then have "x. x  C  U x x = 1"
    by (metis recurrent_iff_U_eq_1)

  { fix y assume "y  C"
    then have "U y y = 1" "recurrent y"
      using y  C  U y y = 1 all_recurrent by auto
    have "measure N {y}  enn2real (1 / U' y y)"
    proof (rule field_le_epsilon)
      fix e :: real assume "0 < e"
      from eps[OF 0 < e] y  C obtain A where
        A: "finite A" "A  C" and
        le: "z. 0 < z  z < 1  measure N {y}  (1 - z) / (1 - gf_U y y z) * (x. gf_F x y z * indicator A x N) + e"
        by auto
      let ?L = "at_left (1::real)"
      have "((λz. (1 - z) / (1 - gf_U y y z) * (x. gf_F x y z * indicator A x N) + e) 
          enn2real (1 / U' y y) * (x. F x y * indicator A x N) + e) ?L"
      proof (intro tendsto_add tendsto_const tendsto_mult int_gf_F,
             rule lhopital_left[where f'="λx. - 1" and g'="λz. - gf_U' y y z"])
        show "((-) 1  0) ?L" "((λx. 1 - gf_U y y x)  0) ?L"
          using gf_U[of y y] by (auto intro!: tendsto_eq_intros simp: U y y = 1)
        show "y  C" "finite A" "A  C" by fact+
        show "eventually (λx. 1 - gf_U y y x  0) ?L"
          using gf_G_eq_gf_U(2)[OF convergence_G_less_1, where 'z=real] by (auto intro!: eventually_at_left_1)
        show "((λx. - 1 / - gf_U' y y x)  enn2real (1 / U' y y)) ?L"
          using recurrent y by (rule inverse_gf_U'_tendsto)
        have "eventually (λx. 0 < gf_U' y y x) ?L"
          by (intro eventually_at_left_1 gf_U'_pos) (simp_all add: U y y = 1)
        then show "eventually (λx. - gf_U' y y x  0) ?L"
          by eventually_elim simp
        show "eventually (λx. DERIV (λx. 1 - gf_U y y x) x :> - gf_U' y y x) ?L"
          by (auto intro!: eventually_at_left_1 derivative_eq_intros DERIV_gf_U)
        show "eventually (λx. DERIV ((-) 1) x :> - 1) ?L"
          by (auto intro!: eventually_at_left_1 derivative_eq_intros)
      qed
      then have "measure N {y}  enn2real (1 / U' y y) * (x. F x y * indicator A x N) + e"
        by (rule tendsto_le[OF trivial_limit_at_left_real _ tendsto_const]) (intro eventually_at_left_1 le)
      then have "measure N {y} - e  enn2real (1 / U' y y) * (x. F x y * indicator A x N)"
        by simp
      also have "  enn2real (1 / U' y y)"
        using A
        by (intro mult_left_le measure_pmf.integral_le_const measure_pmf.integrable_const_bound[where B=1])
           (auto simp: mult_le_one F_le_1 U'_def)
      finally show "measure N {y}  enn2real (1 / U' y y) + e"
        by simp
    qed }
  note measure_y_le = this

  show pos: "yC. pos_recurrent y"
  proof (rule ccontr)
    assume "¬ (yC. pos_recurrent y)"
    then obtain x where x: "x  C" "¬ pos_recurrent x" by auto
    { fix y assume "y  C"
      with x have "¬ pos_recurrent y"
        using C by (auto simp: essential_class_def pos_recurrent_iffI_communicating[symmetric] elim!: quotientE)
      with all_recurrent y  C have "enn2real (1 / U' y y) = 0"
        by (simp add: pos_recurrent_def nn_integral_add)
      with measure_y_le[OF y  C] have "measure N {y} = 0"
        by (auto intro!: antisym simp: pos_recurrent_def) }
    then have "emeasure N C = 0"
      by (subst emeasure_countable_singleton) (auto simp: C ae_C measure_pmf.emeasure_eq_measure nn_integral_0_iff_AE)
    then show False
      using measure N C = 1 by (simp add: measure_pmf.emeasure_eq_measure)
  qed

  { fix A :: "'s set" assume [simp]: "countable A"
    have "emeasure N A = (+x. emeasure N {x} count_space A)"
      by (intro emeasure_countable_singleton) auto
    also have "  (+x. emeasure (stat C) {x} count_space A)"
    proof (intro nn_integral_mono)
      fix y assume "y  space (count_space A)"
      show "emeasure N {y}  emeasure (stat C) {y}"
      proof cases
        assume "y  C"
        with pos have "pos_recurrent y"
          by auto
        with one_le_integral_t[of y] obtain r where r: "U' y y = ennreal r" "1  U' y y" and [simp]: "0  r"
          by (cases "U' y y") (auto simp: pos_recurrent_def nn_integral_add)

        from measure_y_le[OF y  C]
        have "emeasure N {y}  ennreal (enn2real (1 / U' y y))"
          by (simp add: measure_pmf.emeasure_eq_measure)
        also have " = emeasure (stat C) {y}"
          unfolding stat_def using y  C r
          by (subst emeasure_point_measure_finite2)
             (auto simp add: ennreal_1[symmetric] divide_ennreal inverse_ennreal inverse_eq_divide ennreal_mult[symmetric]
                   simp del: ennreal_1)
        finally show "emeasure N {y}  emeasure (stat C) {y}"
          by simp
      next
        assume "y  C"
        with ae_C have "emeasure N {y} = 0"
          by (subst AE_iff_measurable[symmetric, where P="λx. x  y"]) (auto elim!: eventually_mono)
        moreover have "emeasure (stat C) {y} = 0"
          using emeasure_stat_not_C[OF y  C] .
        ultimately show ?thesis by simp
      qed
    qed
    also have " = emeasure (stat C) A"
      by (intro emeasure_countable_singleton[symmetric]) auto
    finally have "emeasure N A  emeasure (stat C) A" . }
  note N_le_C = this

  from stat_subprob[OF C(1) countable C pos] N_le_C[OF countable C] measure N C = 1
  have stat_C_eq_1: "emeasure (stat C) C = 1"
    by (auto simp add: measure_pmf.emeasure_eq_measure one_ennreal_def)
  moreover have "emeasure (stat C) (UNIV - C) = 0"
    by (subst AE_iff_measurable[symmetric, where P="λx. x  C"])
       (auto simp: stat_def AE_point_measure sets_point_measure space_point_measure
                split: split_indicator cong del: AE_cong)
  ultimately have "emeasure (stat C) (space (stat C)) = 1"
    using plus_emeasure[of C "stat C" "UNIV - C"] by (simp add: Un_absorb1)
  interpret stat: prob_space "stat C"
    by standard fact

  show "measure_pmf N = stat C"
  proof (rule measure_eqI_countable_AE)
    show "sets N = UNIV" "sets (stat C) = UNIV"
      by auto
    show "countable C" "AE x in N. x  C" and ae_stat: "AE x in stat C. x  C"
      using C ae_C stat_C_eq_1 by (auto intro!: stat.AE_prob_1 simp: stat.emeasure_eq_measure)

    { assume "x. emeasure N {x}  emeasure (stat C) {x}"
      then obtain x where [simp]: "emeasure N {x}  emeasure (stat C) {x}" by auto
      with N_le_C[of "{x}"] have x: "emeasure N {x} < emeasure (stat C) {x}"
        by (auto simp: less_le)
      have "1 = emeasure N {x} + emeasure N (C - {x})"
        using ae_C
        by (subst plus_emeasure) (auto intro!: measure_pmf.emeasure_eq_1_AE)
      also have " < emeasure (stat C) {x} + emeasure (stat C) (C - {x})"
        using x N_le_C[of "C - {x}"] C ae_C
        by (simp add: stat.emeasure_eq_measure measure_pmf.emeasure_eq_measure
                      ennreal_plus[symmetric] ennreal_less_iff
                 del: ennreal_plus)
      also have " = 1"
        using ae_stat by (subst plus_emeasure) (auto intro!: stat.emeasure_eq_1_AE)
      finally have False by simp }
    then show "x. emeasure N {x} = emeasure (stat C) {x}" by auto
  qed
qed

lemma measure_point_measure_singleton:
  "x  A  measure (point_measure A X) {x} = enn2real (X x)"
  unfolding measure_def by (subst emeasure_point_measure_finite2) auto

lemma stationary_distribution_imp_int_t:
  assumes C: "essential_class C" "countable C" "stationary_distribution N" "N  C"
  assumes x: "x  C" shows "U' x x = 1 / ennreal (pmf N x)"
proof -
  from stationary_distributionD[OF C]
  have "measure_pmf N = stat C" and *: "xC. pos_recurrent x" by auto
  show ?thesis
    unfolding measure_pmf N = stat C pmf.rep_eq stat_def
    using *[THEN bspec, OF x] x
    apply (simp add: measure_point_measure_singleton)
    apply (cases "U' x x")
    subgoal for r
      by (cases "r = 0")
         (simp_all add: divide_ennreal_def inverse_ennreal)
    apply simp
    done
qed

definition "period_set x = {i. 0 < i  0 < p x x i }"
definition "period C = (SOME d. xC. d = Gcd (period_set x))"

lemma Gcd_period_set_invariant:
  assumes c: "(x, y)  communicating"
  shows "Gcd (period_set x) = Gcd (period_set y)"
proof -
  { fix x y n assume c: "(x, y)  communicating" "x  y" and n: "n  period_set x"
    from c obtain l k where "0 < p x y l" "0 < p y x k"
      by (auto simp: communicating_def dest!: accD_pos)
    moreover with x  y have "l  0  k  0"
      by (intro notI conjI) (auto simp: p_0)
    ultimately have pos: "0 < l" "0 < k" and l: "0 < p x y l" and k: "0 < p y x k"
      by auto

    from mult_pos_pos[OF k l] prob_reachable_le[of k "k + l" y x y] c
    have k_l: "0 < p y y (k + l)"
      by simp
    then have "Gcd (period_set y) dvd k + l"
      using pos by (auto intro!: Gcd_dvd_nat simp: period_set_def)
    moreover
    from n have "0 < p x x n" "0 < n" by (auto simp: period_set_def)
    from mult_pos_pos[OF k this(1)] prob_reachable_le[of k "k + n" y x x] c
    have "0 < p y x (k + n)"
      by simp
    from mult_pos_pos[OF this(1) l] prob_reachable_le[of "k + n" "(k + n) + l" y x y] c
    have "0 < p y y (k + n + l)"
      by simp
    then have "Gcd (period_set y) dvd (k + l) + n"
      using pos by (auto intro!: Gcd_dvd_nat simp: period_set_def ac_simps)
    ultimately have "Gcd (period_set y) dvd n"
      by (metis dvd_add_left_iff add.commute) }
  note this[of x y] this[of y x] c
  moreover have "(y, x)  communicating"
    using c by (simp add: communicating_def)
  ultimately show ?thesis
    by (auto intro: dvd_antisym Gcd_greatest Gcd_dvd)
qed

lemma period_eq:
  assumes "C  UNIV // communicating" "x  C"
  shows "period C = Gcd (period_set x)"
  unfolding period_def
  using assms
  by (rule_tac someI2[where a="Gcd (period_set x)"])
     (auto intro!: Gcd_period_set_invariant irreducibleD)

definition "aperiodic C  C  UNIV // communicating  period C = 1"

definition "not_ephemeral C  C  UNIV // communicating  ¬ (x. C = {x}  p x x 1 = 0)"

lemma not_ephemeralD:
  assumes C: "not_ephemeral C" "x  C"
  shows "n>0. 0 < p x x n"
proof cases
  assume "x. C = {x}"
  with x  C have "C = {x}" by auto
  with C p_nonneg[of x x 1] have "0 < p x x 1"
    by (auto simp: not_ephemeral_def less_le)
  with C = {x} show ?thesis by auto
next
  from C have irr: "C  UNIV // communicating"
    by (auto simp: not_ephemeral_def)
  assume "¬(x. C = {x})"
  then have "x. C  {x}" by auto
  with x  C obtain y where "y  C" "x  y"
    by blast
  with irreducibleD[OF irr, of x y] C x  C have c: "(x, y)  communicating" by auto
  with accD_pos[of x y] accD_pos[of y x]
  obtain k l where pos: "0 < p x y k" "0 < p y x l"
    by (auto simp: communicating_def)
  with x  y have "l  0"
    by (intro notI) (auto simp: p_0)
  have "0 < p x y k * p y x (k + l - k)"
    using pos by auto
  also have "p x y k * p y x (k + l - k)  p x x (k + l)"
    using prob_reachable_le[of "k" "k + l" x y x] c by auto
  finally show ?thesis
    using l  0 x  C by (auto intro!: exI[of _ "k + l"])
qed

lemma not_ephemeralD_pos_period:
  assumes C: "not_ephemeral C"
  shows "0 < period C"
proof -
  from C not_empty_irreducible[of C] obtain x where "x  C"
    by (auto simp: not_ephemeral_def)
  from not_ephemeralD[OF C this]
  obtain n where n: "0 < p x x n" "0 < n" by auto
  have C': "C  UNIV // communicating"
    using C by (auto simp: not_ephemeral_def)

  have "period C  0"
    unfolding period_eq [OF C' x  C]
    using n by (auto simp: period_set_def)
  then show ?thesis by auto
qed


lemma period_posD:
  assumes C: "C  UNIV // communicating" and "0 < period C" "x  C"
  shows "n>0. 0 < p x x n"
proof -
  from 0 < period C have "period C  0"
    by auto
  then show ?thesis
    unfolding period_eq [OF C x  C]
    unfolding period_set_def by auto
qed

lemma not_ephemeralD_pos_period':
  assumes C: "C  UNIV // communicating"
  shows "not_ephemeral C  0 < period C"
proof (auto dest!: not_ephemeralD_pos_period intro: C)
  from C not_empty_irreducible[of C] obtain x where "x  C"
    by (auto simp: not_ephemeral_def)

  assume "0 < period C"
  then show "not_ephemeral C"
    apply (auto simp: not_ephemeral_def C)
oops ― ‹should be easy to finish›


lemma eventually_periodic:
  assumes C: "C  UNIV // communicating" "0 < period C" "x  C"
  shows "eventually (λm. 0 < p x x (m * period C)) sequentially"
proof -
  from period_posD[OF assms] obtain n where n: "0 < p x x n" "0 < n" by auto
  have C': "C  UNIV // communicating"
    using C by auto

  have "period C  0"
    unfolding period_eq [OF C' x  C]
    using n by (auto simp: period_set_def)
  have "eventually (λm. m * Gcd (period_set x)  (period_set x)) sequentially"
  proof (rule eventually_mult_Gcd)
    show "n > 0" "n  period_set x"
      using n by (auto simp add: period_set_def)
    fix k l  assume "k  period_set x" "l  period_set x"
    then have "0 < p x x k * p x x l" "0 < l" "0 < k"
      by (auto simp: period_set_def)
    moreover have "p x x k * p x x l  p x x (k + l)"
      using prob_reachable_le[of k "k + l" x x x] x  C
      by auto
    ultimately show "k + l  period_set x"
      using 0 < l by (auto simp: period_set_def)
  qed
  with eventually_ge_at_top[of 1] show "eventually (λm. 0 < p x x (m * period C)) sequentially"
    by eventually_elim 
       (insert period C  0 period_eq[OF C' x  C, symmetric], auto simp: period_set_def)
qed


lemma aperiodic_eventually_recurrent:
  "aperiodic C  C  UNIV // communicating  (xC. eventually (λm. 0 < p x x m) sequentially)"
proof safe
  fix x assume "x  C" "aperiodic C"
  with eventually_periodic[of C x]
  show "eventually (λm. 0 < p x x m) sequentially"
    by (auto simp add: aperiodic_def)
next
  assume "xC. eventually (λm. 0 < p x x m) sequentially" and C: "C  UNIV // communicating"
  moreover from not_empty_irreducible[OF C] obtain x where "x  C" by auto
  ultimately obtain N where "M.  MN  0 < p x x M"
    by (auto simp: eventually_sequentially)
  then have "{N <..}  period_set x"
    by (auto simp: period_set_def)
  from C show "aperiodic C"
    unfolding period_eq [OF C x  C] aperiodic_def
  proof
    show "Gcd (period_set x) = 1"
    proof (rule Gcd_eqI)
      from one_dvd show "1 dvd q" for q :: nat .
      fix m
      assume "q. q  period_set x  m dvd q"
      moreover from {N <..}  period_set x
      have "{Suc N, Suc (Suc N)}  period_set x"
        by auto
      ultimately have "m dvd Suc (Suc N)" and "m dvd Suc N"
        by auto
      then have "m dvd Suc (Suc N) - Suc N"
        by (rule dvd_diff_nat)
      then show "is_unit m"
        by simp
    qed simp
  qed
qed (simp add: aperiodic_def)

lemma stationary_distributionD_emeasure:
  assumes N: "stationary_distribution N"
  shows "emeasure N A = (+s. emeasure (K s) A N)"
proof -
  have "prob_space (measure_pmf N)"
    by intro_locales
  then interpret subprob_space "measure_pmf N"
    by (rule prob_space_imp_subprob_space)
  show ?thesis
    unfolding measure_pmf.emeasure_eq_measure
    apply (subst N[unfolded stationary_distribution_def])
    apply (simp add: measure_pmf_bind)
    apply (subst measure_pmf.measure_bind[where N="count_space UNIV"])
    apply (rule measurable_compose[OF _ measurable_measure_pmf])
    apply (auto intro!: nn_integral_eq_integral[symmetric] measure_pmf.integrable_const_bound[where B=1])
    done
qed

lemma communicatingD1:
  "C  UNIV // communicating  (a, b)  communicating  a  C  b  C"
  by (auto elim!: quotientE) (auto simp add: communicating_def)

lemma communicatingD2:
  "C  UNIV // communicating  (a, b)  communicating  b  C  a  C"
  by (auto elim!: quotientE) (auto simp add: communicating_def)

lemma acc_iff: "(x, y)  acc  (n. 0 < p x y n)"
  by (blast intro: accD_pos accI_pos)

lemma communicating_iff: "(x, y)  communicating  (n. 0 < p x y n)  (n. 0 < p y x n)"
  by (auto simp add: acc_iff communicating_def)

end

context MC_pair
begin

lemma p_eq_p1_p2:
  "p (x1, x2) (y1, y2) n = K1.p x1 y1 n * K2.p x2 y2 n"
  unfolding p_def K1.p_def K2.p_def
  by (subst prod_eq_prob_T)
     (auto intro!: arg_cong2[where f=measure] split: nat.splits simp: Stream_snth)

lemma P_accD:
  assumes "((x1, x2), (y1, y2))  acc"shows "(x1, y1)  K1.acc" "(x2, y2)  K2.acc"
  using assms by (auto simp: acc_iff K1.acc_iff K2.acc_iff p_eq_p1_p2 zero_less_mult_iff not_le[of 0, symmetric]
                       cong: conj_cong)

lemma aperiodicI_pair:
  assumes C1: "K1.aperiodic C1" and C2: "K2.aperiodic C2"
  shows "aperiodic (C1 × C2)"
  unfolding aperiodic_eventually_recurrent
proof safe
  from C1[unfolded K1.aperiodic_eventually_recurrent] C2[unfolded K2.aperiodic_eventually_recurrent]
  have C1: "C1  UNIV // K1.communicating" and C2: "C2  UNIV // K2.communicating" and
    ev: "x. x  C1  eventually (λm. 0 < K1.p x x m) sequentially" "x. x  C2  eventually (λm. 0 < K2.p x x m) sequentially"
    by auto
  { fix x1 x2 assume x: "x1  C1" "x2  C2"
    from ev(1)[OF x(1)] ev(2)[OF x(2)]
    show "eventually (λm. 0 < p (x1, x2) (x1, x2) m) sequentially"
       by eventually_elim  (simp add: p_eq_p1_p2 x) }

  { fix x1 x2 y1 y2
    assume acc: "(x1, y1)  K1.acc" "(x2, y2)  K2.acc" "x1  C1" "y1  C1" "x2  C2" "y2  C2"
    then obtain k l where "0 < K1.p x1 y1 l" "0 < K2.p x2 y2 k"
      by (auto dest!: K1.accD_pos K2.accD_pos)
    with acc ev(1)[of y1] ev(2)[of y2]
    have "eventually (λm. 0 < K1.p x1 y1 l * K1.p y1 y1 m  0 < K2.p x2 y2 k * K2.p y2 y2 m) sequentially"
      by (auto elim: eventually_elim2)
    then have "eventually (λm. 0 < K1.p x1 y1 (m + l)  0 < K2.p x2 y2 (m + k)) sequentially"
    proof eventually_elim
      fix m assume "0 < K1.p x1 y1 l * K1.p y1 y1 m  0 < K2.p x2 y2 k * K2.p y2 y2 m"
      with acc
        K1.prob_reachable_le[of l "l + m" x1 y1 y1]
        K2.prob_reachable_le[of k "k + m" x2 y2 y2]
      show "0 < K1.p x1 y1 (m + l)  0 < K2.p x2 y2 (m + k)"
        by (auto simp add: ac_simps)
    qed
    then have "eventually (λm. 0 < K1.p x1 y1 m  0 < K2.p x2 y2 m) sequentially"
      unfolding eventually_conj_iff by (subst (asm) (1 2) eventually_sequentially_seg) (auto elim: eventually_elim2)
    then obtain N where "0 < K1.p x1 y1 N" "0 < K2.p x2 y2 N"
      by (auto simp: eventually_sequentially)
    with acc have "0 < p (x1, x2) (y1, y2) N"
      by (auto simp add: p_eq_p1_p2)
    with acc have "((x1, x2), (y1, y2))  acc"
      by (auto intro!: accI_pos) }
  note 1 = this

  { fix x1 x2 y1 y2 assume acc:"((x1, x2), (y1, y2))  acc"
    moreover from acc obtain k where "0 < p (x1, x2) (y1, y2) k" by (auto dest!: accD_pos)
    ultimately have "(x1, y1)  K1.acc  (x2, y2)  K2.acc"
      by (subst (asm) p_eq_p1_p2)
         (auto intro!: K1.accI_pos K2.accI_pos simp: zero_less_mult_iff not_le[of 0, symmetric]) }
  note 2 = this

  from K1.not_empty_irreducible[OF C1] K2.not_empty_irreducible[OF C2]
  obtain x1 x2 where xC: "x1  C1" "x2  C2" by auto
  show "C1 × C2  UNIV // communicating"
    apply (simp add: quotient_def Image_def)
    apply (safe intro!: exI[of _ x1] exI[of _ x2])
  proof -
    fix y1 y2 assume yC: "y1  C1" "y2  C2"
    from K1.irreducibleD[OF C1 x1  C1 y1  C1] K2.irreducibleD[OF C2 x2  C2 y2  C2]
    show "((x1, x2), (y1, y2))  communicating"
      using 1[of x1 y1 x2 y2] 1[of y1 x1 y2 x2] xC yC
      by (auto simp: communicating_def K1.communicating_def K2.communicating_def)
  next
    fix y1 y2 assume "((x1, x2), (y1, y2))  communicating"
    with 2[of x1 x2 y1 y2] 2[of y1 y2 x1 x2]
    have "(x1, y1)  K1.communicating" "(x2, y2)  K2.communicating"
      by (auto simp: communicating_def K1.communicating_def K2.communicating_def)
    with xC show "y1  C1" "y2  C2"
      using K1.communicatingD1[OF C1] K2.communicatingD1[OF C2] by auto
  qed
qed

lemma stationary_distributionI_pair:
  assumes N1: "K1.stationary_distribution N1"
  assumes N2: "K2.stationary_distribution N2"
  shows "stationary_distribution (pair_pmf N1 N2)"
  unfolding stationary_distribution_def
  unfolding Kp_def pair_pmf_def
  apply (subst N1[unfolded K1.stationary_distribution_def])
  apply (subst N2[unfolded K2.stationary_distribution_def])
  apply (simp add: bind_assoc_pmf bind_return_pmf)
  apply (subst bind_commute_pmf[of N2])
  apply simp
  done

end

context MC_syntax
begin

lemma stationary_distribution_imp_limit:
  assumes C: "aperiodic C" "essential_class C" "countable C" and N: "stationary_distribution N" "N  C"
  assumes [simp]: "y  C"
  shows "(λn. x. ¦p y x n - pmf N x¦ count_space C)  0"
    (is "?L  0")
proof -
  from essential_class C have C_comm: "C  UNIV // communicating"
    by (simp add: essential_class_def)

  define K' where "K' = (λSome x  map_pmf Some (K x) | None  map_pmf Some N)"

  interpret K2: MC_syntax K' .
  interpret KN: MC_pair K K' .

  from stationary_distributionD[OF C(2,3) N]
  have pos: "x. x  C  pos_recurrent x" and "measure_pmf N = stat C" by auto

  have pos: "x. x  C  0 < emeasure N {x}"
    using pos unfolding stat_def measure_pmf N = stat C
    by (subst emeasure_point_measure_finite2)
       (auto simp: U'_def pos_recurrent_def nn_integral_add ennreal_zero_less_divide less_top)
  then have rpos: "x. x  C  0 < pmf N x"
    by (simp add: measure_pmf.emeasure_eq_measure pmf.rep_eq)

  have eq: "x y. (if x = y then 1 else 0) = indicator {y} x" by auto

  have intK: "f x. (x. (f x :: real) K' (Some x)) = (x. f (Some x) K x)"
    by (simp add: K'_def integral_distr map_pmf_rep_eq)

  { fix m and x y :: 's
    have "K2.p (Some x) (Some y) m = p x y m"
      by (induct m arbitrary: x)
         (auto intro!: integral_cong simp add: K2.p_Suc' p_Suc' intK K2.p_0 p_0) }
  note K_p_eq = this

  { fix n and x :: 's have "K2.p (Some x) None n = 0"
      by (induct n arbitrary: x) (auto simp: K2.p_Suc' K2.p_0 intK cong: integral_cong) }
  note K_S_None = this

  from not_empty_irreducible[OF C_comm] obtain c0 where c0: "c0  C" by auto

  have K2_acc: "x y. (Some x, y)  K2.acc  (z. y = Some z  (x, z)  acc)"
    apply (auto simp: K2.acc_iff acc_iff K_p_eq)
    apply (case_tac y)
    apply (auto simp: K_p_eq K_S_None)
    done

  have K2_communicating: "c x. c  C  (Some c, x)  K2.communicating  (c'C. x = Some c')"
  proof safe
    fix x c assume "c  C" "(Some c, x)  K2.communicating"
    then show "c'C. x = Some c'"
      by (cases x)
         (auto simp: communicating_iff K2.communicating_iff K_p_eq K_S_None intro!: irreducibleD2[OF C_comm cC])
  next
    fix c c' x assume "c  C" "c'  C"
    with irreducibleD[OF C_comm this] show "(Some c, Some c')  K2.communicating"
      by (auto simp: K2.communicating_iff communicating_iff K_p_eq)
  qed

  have "Some ` C  UNIV // K2.communicating"
    by (auto simp add: quotient_def Image_def c0 K2_communicating
             intro!: exI[of _ "Some c0"])
  then have "K2.essential_class (Some ` C)"
    by (rule K2.essential_classI)
       (auto simp: K2_acc essential_classD2[OF essential_class C])

  have "K2.aperiodic (Some ` C)"
    unfolding K2.aperiodic_eventually_recurrent
  proof safe
    fix x assume "x  C" then show "eventually (λm. 0 < K2.p (Some x) (Some x) m) sequentially"
      using aperiodic C unfolding aperiodic_eventually_recurrent
      by (auto elim!: eventually_mono simp: K_p_eq)
  qed fact
  then have aperiodic: "KN.aperiodic (C × Some ` C)"
    by (rule KN.aperiodicI_pair[OF aperiodic C])

  have KN_essential: "KN.essential_class (C × Some ` C)"
  proof (rule KN.essential_classI)
    show "C × Some ` C  UNIV // KN.communicating"
      using aperiodic by (simp add: KN.aperiodic_def)
  next
    fix x y assume "x  C × Some ` C" "(x, y)  KN.acc"
    with KN.P_accD[of "fst x" "snd x" "fst y" "snd y"]
    show "y  C × Some ` C"
      by (cases x y rule: prod.exhaust[case_product prod.exhaust])
         (auto simp: K2_acc essential_classD2[OF essential_class C])
  qed

  { fix n and x y :: 's
    have "measure N {y} = 𝒫(ω in K2.T None. (None ## ω) !! (Suc n) = Some y)"
      unfolding stationary_distribution_iterate'[OF N(1), of y n]
      apply (subst K2.p_def[symmetric])
      apply (subst K2.p_Suc')
      apply (subst K'_def)
      apply (simp add: map_pmf_rep_eq integral_distr K_p_eq)
      done
    then have "measure N {y} = 𝒫(ω in K2.T None. ω !! n = Some y)"
      by simp }
  note measure_y_eq = this

  define D where "D = {x::'s × 's option. Some (fst x) = snd x}"

  have [measurable]:
    "P::('s × 's option  bool). P  measurable (count_space UNIV) (count_space UNIV)"
    by simp

  { fix n and x :: 's
    have "𝒫(ω in KN.T (y, None). i<n. snd (ω !! n) = Some x  ev_at (HLD D) i ω) =
      (i<n. 𝒫(ω in KN.T (y, None). snd (ω !! n) = Some x  ev_at (HLD D) i ω))"
      by (subst KN.T.finite_measure_finite_Union[symmetric])
         (auto simp: disjoint_family_on_def intro!: arg_cong2[where f=measure] dest: ev_at_unique)
    also have " = (i<n. 𝒫(ω in KN.T (y, None). fst (ω !! n) = x  ev_at (HLD D) i ω))"
    proof (intro sum.cong refl)
      fix i assume i: "i  {..< n}"
      show "𝒫(ω in KN.T (y, None). snd (ω !! n) = Some x  ev_at (HLD D) i ω) =
        𝒫(ω in KN.T (y, None). fst (ω !! n) = x  ev_at (HLD D) i ω)"
        apply (subst (1 2) KN.prob_T_split[where n="Suc i"])
        apply (simp_all add: ev_at_shift snth_Stream del: stake.simps KN.space_T)
        unfolding ev_at_shift snth_Stream
      proof (intro Bochner_Integration.integral_cong refl)
        fix ω :: "('s × 's option) stream" let ?s = "λω'. stake (Suc i) ω @- ω'"
        show "𝒫(ω' in KN.T (ω !! i). snd (?s ω' !! n) = Some x  ev_at (HLD D) i ω) =
          𝒫(ω' in KN.T (ω !! i). fst (?s ω' !! n) = x  ev_at (HLD D) i ω)"
        proof cases
          assume "ev_at (HLD D) i ω"
          from ev_at_imp_snth[OF this]
          have eq: "snd (ω !! i) = Some (fst (ω !! i))"
            by (simp add: D_def HLD_iff)

          have "𝒫(ω' in KN.T (ω !! i). fst (ω' !! (n - Suc i)) = x) =
            𝒫(ω' in T (fst (ω !! i)). ω' !! (n - Suc i) = x) * 𝒫(ω' in K2.T (snd (ω !! i)). True)"
            by (subst KN.prod_eq_prob_T) simp_all
          also have " = p (fst (ω !! i)) x (Suc (n - Suc i))"
            using K2.T.prob_space by (simp add: p_def)
          also have " = K2.p (snd (ω !! i)) (Some x) (Suc (n - Suc i))"
            by (simp add: K_p_eq eq)
          also have " = 𝒫(ω' in T (fst (ω !! i)). True) * 𝒫(ω' in K2.T (snd (ω !! i)). ω' !! (n - Suc i) = Some x)"
            using T.prob_space by (simp add: K2.p_def)
          also have " = 𝒫(ω' in KN.T (ω !! i). snd (ω' !! (n - Suc i)) = Some x)"
            by (subst KN.prod_eq_prob_T) simp_all
          finally show ?thesis using ev_at (HLD D) i ω i
            by (simp del: stake.simps)
        qed simp
      qed
    qed
    also have " = 𝒫(ω in KN.T (y, None). (i<n. fst (ω !! n) = x  ev_at (HLD D) i ω))"
      by (subst KN.T.finite_measure_finite_Union[symmetric])
         (auto simp add: disjoint_family_on_def dest: ev_at_unique
               intro!: arg_cong2[where f=measure])
    finally have eq: "𝒫(ω in KN.T (y, None). (i<n. snd (ω !! n) = Some x  ev_at (HLD D) i ω)) =
      𝒫(ω in KN.T (y, None). (i<n. fst (ω !! n) = x  ev_at (HLD D) i ω))" .

    have "p y x (Suc n) - measure N {x} = 𝒫(ω in T y. ω !! n = x) - 𝒫(ω in K2.T None. ω !! n = Some x)"
      unfolding p_def by (subst measure_y_eq) simp_all
    also have "𝒫(ω in T y. ω !! n = x) = 𝒫(ω in T y. ω !! n = x) * 𝒫(ω in K2.T None. True)"
      using K2.T.prob_space by simp
    also have " = 𝒫(ω in KN.T (y, None). fst (ω !! n) = x)"
      by (subst KN.prod_eq_prob_T) auto
    also have " = 𝒫(ω in KN.T (y, None). (i<n. fst (ω !! n) = x  ev_at (HLD D) i ω)) +
      𝒫(ω in KN.T (y, None). fst (ω !! n) = x  ¬ (i<n. ev_at (HLD D) i ω))"
      by (subst KN.T.finite_measure_Union[symmetric])
         (auto intro!: arg_cong2[where f=measure])
    also have "𝒫(ω in K2.T None. ω !! n = Some x) = 𝒫(ω in T y. True) * 𝒫(ω in K2.T None. ω !! n = Some x)"
      using T.prob_space by simp
    also have " = 𝒫(ω in KN.T (y, None). snd (ω !! n) = Some x)"
      by (subst KN.prod_eq_prob_T) auto
    also have " = 𝒫(ω in KN.T (y, None). (i<n. snd (ω !! n) = Some x  ev_at (HLD D) i ω)) +
      𝒫(ω in KN.T (y, None). snd (ω !! n) = Some x  ¬ (i<n. ev_at (HLD D) i ω))"
      by (subst KN.T.finite_measure_Union[symmetric])
         (auto intro!: arg_cong2[where f=measure])
    finally have "¦ p y x (Suc n) - measure N {x} ¦ =
      ¦ 𝒫(ω in KN.T (y, None). fst (ω !! n) = x  ¬ (i<n. ev_at (HLD D) i ω)) -
      𝒫(ω in KN.T (y, None). snd (ω !! n) = Some x  ¬ (i<n. ev_at (HLD D) i ω)) ¦"
      unfolding eq by (simp add: field_simps)
    also have "  ¦ 𝒫(ω in KN.T (y, None). fst (ω !! n) = x  ¬ (i<n. ev_at (HLD D) i ω)) ¦ +
      ¦ 𝒫(ω in KN.T (y, None). snd (ω !! n) = Some x  ¬ (i<n. ev_at (HLD D) i ω)) ¦"
      by (rule abs_triangle_ineq4)
    also have "  𝒫(ω in KN.T (y, None). fst (ω !! n) = x  ¬ (i<n. ev_at (HLD D) i ω)) +
      𝒫(ω in KN.T (y, None). snd (ω !! n) = Some x  ¬ (i<n. ev_at (HLD D) i ω))"
      by simp
    finally have "¦ p y x (Suc n) - measure N {x} ¦  " . }
  note mono = this

  { fix n :: nat
    have "(+x. ¦ p y x (Suc n) - measure N {x} ¦ count_space C) 
      (+x. ennreal (𝒫(ω in KN.T (y, None). fst (ω !! n) = x  ¬ (i<n. ev_at (HLD D) i ω))) +
      ennreal (𝒫(ω in KN.T (y, None). snd (ω !! n) = Some x  ¬ (i<n. ev_at (HLD D) i ω))) count_space C)"
      using mono by (intro nn_integral_mono) (simp add: ennreal_plus[symmetric] del: ennreal_plus)
    also have " = (+x. 𝒫(ω in KN.T (y, None). fst (ω !! n) = x  ¬ (i<n. ev_at (HLD D) i ω)) count_space C) +
      (+x. 𝒫(ω in KN.T (y, None). snd (ω !! n) = Some x  ¬ (i<n. ev_at (HLD D) i ω)) count_space C)"
      by (subst nn_integral_add) auto
    also have " = emeasure (KN.T (y, None)) (xC. {ωspace (KN.T (y, None)). fst (ω !! n) = x  ¬ (i<n. ev_at (HLD D) i ω)}) +
      emeasure (KN.T (y, None)) (xC. {ωspace (KN.T (y, None)). snd (ω !! n) = Some x  ¬ (i<n. ev_at (HLD D) i ω)})"
      by (subst (1 2) emeasure_UN_countable)
         (auto simp add: disjoint_family_on_def KN.T.emeasure_eq_measure C)
    also have "  ennreal (𝒫(ω in KN.T (y, None). ¬ (i<n. ev_at (HLD D) i ω))) + ennreal (𝒫(ω in KN.T (y, None). ¬ (i<n. ev_at (HLD D) i ω)))"
      unfolding KN.T.emeasure_eq_measure
      by (intro add_mono) (auto intro!: KN.T.finite_measure_mono)
    also have "  2 * 𝒫(ω in KN.T (y, None). ¬ (i<n. ev_at (HLD D) i ω))"
      by (simp add: ennreal_plus[symmetric] del: ennreal_plus)
    finally have "?L (Suc n)  2 * 𝒫(ω in KN.T (y, None). ¬ (i<n. ev_at (HLD D) i ω))"
      by (auto intro!: integral_real_bounded simp add: pmf.rep_eq) }
  note le_2 = this

  have c0_D: "(c0, Some c0)  D"
    by (simp add: D_def c0)

  let ?N' = "map_pmf Some N"
  interpret NP: pair_prob_space N ?N' ..

  have pos_recurrent: "xC × Some ` C. KN.pos_recurrent x"
  proof (rule KN.stationary_distributionD(1)[OF KN_essential _ KN.stationary_distributionI_pair[OF N(1)]])
    show "K2.stationary_distribution ?N'"
      unfolding K2.stationary_distribution_def
      by (subst N(1)[unfolded stationary_distribution_def])
         (auto intro!: bind_pmf_cong simp: K'_def map_pmf_def bind_assoc_pmf bind_return_pmf)
    show "countable (C × Some`C)"
      using C by auto
    show "set_pmf (pair_pmf N (map_pmf Some N))  C × Some ` C"
      using N  C by auto
  qed

  from c0_D have "𝒫(ω in KN.T (y, None). alw (not (HLD D)) ω)  𝒫(ω in KN.T (y, None). alw (not (HLD {(c0, Some c0)})) ω)"
    apply (auto intro!: KN.T.finite_measure_mono)
    apply (rule alw_mono, assumption)
    apply (auto simp: HLD_iff)
    done
  also have " = 0"
    apply (rule KN.T.prob_eq_0_AE)
    apply (simp add: not_ev_iff[symmetric])
    apply (subst KN.AE_T_iff)
    apply simp
  proof
    fix t assume t: "t  KN.Kp (y, None)"
    then obtain a b where t_eq: "t = (a, Some b)" "a  K y" "b  N"
      unfolding KN.Kp_def by (auto simp: K'_def)
    with y  C have "a  C"
      using essential_classD2[OF essential_class C y  C] by auto
    have "b  C"
      using N  C b  N by auto

    from pos_recurrent[THEN bspec, of "(c0, Some c0)"]
    have recurrent_c0: "KN.recurrent (c0, Some c0)"
      by (simp add: KN.pos_recurrent_def c0)
    have "C × Some ` C  UNIV // KN.communicating"
      using aperiodic by (simp add: KN.aperiodic_def)
    then have "((c0, Some c0), t)  KN.communicating"
      by (rule KN.irreducibleD) (simp_all add: t_eq c0 b  C a  C)
    then have "((c0, Some c0), t)  KN.acc"
      by (simp add: KN.communicating_def)
    then have "KN.U t (c0, Some c0) = 1"
      by (rule KN.recurrent_acc(1)[OF recurrent_c0])
    then show "AE ω in KN.T t. ev (HLD {(c0, Some c0)}) (t ## ω)"
      unfolding KN.U_def by (subst (asm) KN.T.prob_Collect_eq_1) (auto simp add: ev_Stream)
  qed
  finally have "𝒫(ω in KN.T (y, None). alw (not (HLD D)) ω) = 0"
    by (intro antisym measure_nonneg)

  have "(λn. 𝒫(ω in KN.T (y, None). ¬ (i<n. ev_at (HLD D) i ω))) 
    measure (KN.T (y, None)) (n. {ωspace (KN.T (y, None)). ¬ (i<n. ev_at (HLD D) i ω)})"
    by (rule KN.T.finite_Lim_measure_decseq) (auto simp: decseq_def)
  also have "(n. {ωspace (KN.T (y, None)). ¬ (i<n. ev_at (HLD D) i ω)}) =
    {ωspace (KN.T (y, None)). alw (not (HLD D)) ω}"
    by (auto simp: not_ev_iff[symmetric] ev_iff_ev_at)
  also have "𝒫(ω in KN.T (y, None). alw (not (HLD D)) ω) = 0" by fact
  finally have *: "(λn. 2 * 𝒫(ω in KN.T (y, None). ¬ (i<n. ev_at (HLD D) i ω)))  0"
    by (intro tendsto_eq_intros) auto

  show ?thesis
    apply (rule LIMSEQ_imp_Suc)
    apply (rule tendsto_sandwich[OF _ _ tendsto_const *])
    using le_2
    apply (simp_all add: integral_nonneg_AE)
    done
qed

lemma stationary_distribution_imp_p_limit:
  assumes "aperiodic C" "essential_class C" and [simp]: "countable C"
  assumes N: "stationary_distribution N" "N  C"
  assumes [simp]: "x  C" "y  C"
  shows "p x y  pmf N y"
proof -
  define D where "D y n = ¦p x y n - pmf N y¦" for y n

  from stationary_distribution_imp_limit[OF assms(1,2,3,4,5,6)]
  have INT: "(λn. y. D y n count_space C)  0"
    unfolding D_def .

  { fix n
    have "D y n  (z. D y n * indicator {y} z count_space C)"
      by simp
    also have "  (y. D y n count_space C)"
      by (intro integral_mono)
         (auto split: split_indicator simp: D_def p_def disjoint_family_on_def
               intro!: Bochner_Integration.integrable_diff integrable_pmf T.integrable_measure)
    finally have "D y n  (y. D y n count_space C)" . }
  note * = this

  have D_nonneg: "n. 0  D y n" by (simp add: D_def)

  have "D y  0"
    by (rule tendsto_sandwich[OF _ _ tendsto_const INT])
       (auto simp: eventually_sequentially * D_nonneg)
  then show ?thesis
    using Lim_null[where l="pmf N y" and net=sequentially and f="p x y"]
    by (simp add: D_def [abs_def] tendsto_rabs_zero_iff)
qed

end

lemma (in MC_syntax) essential_classI2:
  assumes "X  {}"
  assumes accI: "x y. x  X  y  X  (x, y)  acc"
  assumes ED: "x y. x  X  y  set_pmf (K x)  y  X"
  shows "essential_class X"
proof (rule essential_classI)
  { fix x y assume "(x, y)  acc" "x  X"
    then show "y  X"
      by induct (auto dest: ED)}
  note accD = this
  from X  {} obtain x where "x  X" by auto
  from x  X show "X  UNIV // communicating"
    by (auto simp add: quotient_def Image_def communicating_def accI dest: accD intro!: exI[of _ x])
qed

end