CubicPathwidth.T2.StarFour

Source file CubicPathwidth/T2/StarFour.lean of the Lean project (how it is built).

import CubicPathwidth.T2.Defs
import CubicPathwidth.C3.StarFrontier
import Mathlib.Analysis.InnerProductSpace.GramMatrix

/-!
# The 4-star comparison (O2-O5), conditional on three hypotheses

`starFour_of : SlepianClaim → BridgeClaim → SimplexFourClaim → StarFourClaim`.

Route:
1. By the bridge, probabilities of events about `(⟨w_j, x⟩)_j` depend only on the
   Gram matrix, so vectors may be zero-padded into `Sample (m + 4)`.
2. The frontier event is the center-above event minus the all-above event.
3. Lowering every incident correlation to `r` by `c_i = b_i - (ρ_i - r) a + √(ρ_i² - r²) ξ_i`
   (fresh orthonormal `ξ_i`) lowers the leaf-leaf covariances; Slepian compares.
4. Conditioning on `⟨a, x⟩ = z` (disintegration) reduces the star `c` and the ideal star
   to orthants of unit-vector families; the simplex hypothesis compares them.
-/

noncomputable section
open scoped RealInnerProductSpace
open MeasureTheory ProbabilityTheory

namespace CubicPathwidth
namespace T2
namespace StarFour

open GaussianField (Sample)

/-- The map `x ↦ (⟨w_j, x⟩)_j`. -/
def iv {m k : ℕ} (w : Fin k → Sample m) (x : Sample m) : EuclideanSpace ℝ (Fin k) :=
  WithLp.toLp 2 fun j => ⟪w j, x⟫

theorem measurable_iv {m k : ℕ} (w : Fin k → Sample m) : Measurable (iv w) :=
  ((PiLp.continuous_toLp 2 _).comp
    (continuous_pi fun _ => continuous_const.inner continuous_id)).measurable

/-- Probability that all coordinates `⟨w_j, x⟩` exceed `τ`. -/
def allAbove {m k : ℕ} (w : Fin k → Sample m) (τ : ℝ) : ℝ :=
  (stdGaussian (Sample m)).real {x | ∀ j, τ < ⟪w j, x⟫}

/-- Strict upper orthant. -/
def Aup (k : ℕ) (θ : ℝ) : Set (EuclideanSpace ℝ (Fin k)) := {y | ∀ i, θ < y i}

theorem measurableSet_Aup (k : ℕ) (θ : ℝ) : MeasurableSet (Aup k θ) := by
  have : Aup k θ = ⋂ i, {y : EuclideanSpace ℝ (Fin k) | θ < y i} := by
    ext y; simp [Aup]
  rw [this]
  exact MeasurableSet.iInter fun i =>
    measurableSet_lt measurable_const (by fun_prop)

theorem allAbove_eq_gaussProb (hb : BridgeClaim) {m k : ℕ} (w : Fin k → Sample m) (τ : ℝ) :
    allAbove w τ = gaussProb (gram w) (Aup k τ) := by
  rw [hb m k w _ (measurableSet_Aup k τ)]
  rfl

theorem allAbove_eq_of_gram_eq (hb : BridgeClaim) {m m' k : ℕ} (w : Fin k → Sample m)
    (w' : Fin k → Sample m') (h : gram w = gram w') (τ : ℝ) :
    allAbove w τ = allAbove w' τ := by
  rw [allAbove_eq_gaussProb hb, allAbove_eq_gaussProb hb, h]

/-! ## Frontier versus all-above -/

/-- The frontier event for a 5-family: center above, some leaf at or below. -/
def frontierProb {M : ℕ} (w : Fin 5 → Sample M) (τ : ℝ) : ℝ :=
  (stdGaussian (Sample M)).real {x | τ < ⟪w 0, x⟫ ∧ ∃ i : Fin 4, ⟪w i.succ, x⟫ ≤ τ}

theorem frontierProb_eq {M : ℕ} (w : Fin 5 → Sample M) (τ : ℝ) :
    frontierProb w τ =
      (stdGaussian (Sample M)).real {x | τ < ⟪w 0, x⟫} - allAbove w τ := by
  unfold frontierProb allAbove
  have hB : MeasurableSet {x : Sample M | ∀ j, τ < ⟪w j, x⟫} := by
    have : {x : Sample M | ∀ j, τ < ⟪w j, x⟫} = ⋂ j, {x | τ < ⟪w j, x⟫} := by
      ext x; simp
    rw [this]
    exact MeasurableSet.iInter fun j => measurableSet_lt measurable_const (by fun_prop)
  have hsub : {x : Sample M | ∀ j, τ < ⟪w j, x⟫} ⊆ {x | τ < ⟪w 0, x⟫} := fun x hx => hx 0
  rw [← measureReal_sdiff hsub hB]
  congr 1
  ext x
  simp only [Set.mem_ofPred_eq, Set.mem_sdiff, not_forall, not_lt]
  constructor
  · rintro ⟨h0, i, hi⟩
    exact ⟨h0, i.succ, hi⟩
  · rintro ⟨h0, j, hj⟩
    refine ⟨h0, ?_⟩
    induction j using Fin.cases with
    | zero => exact absurd h0 (not_lt.mpr hj)
    | succ i => exact ⟨i, hj⟩

theorem prob_unit_gt {M : ℕ} {a : Sample M} (ha : ‖a‖ = 1) (τ : ℝ) :
    (stdGaussian (Sample M)).real {x | τ < ⟪a, x⟫} = (gaussianReal 0 1).real (Set.Ioi τ) := by
  have haa : ⟪a, a⟫ = 1 := by rw [real_inner_self_eq_norm_sq, ha]; norm_num
  have e : {x : Sample M | τ < ⟪a, x⟫} = (fun x => ⟪a, x⟫) ⁻¹' Set.Ioi τ := rfl
  rw [e, ← map_measureReal_apply (by fun_prop) measurableSet_Ioi,
    GaussianAngle.unit_projection_map haa]

/-! ## The ideal star as explicit vectors -/

/-- `√((1 - r²)/3)`. -/
def cs (r : ℝ) : ℝ := √((1 - r ^ 2) / 3)

/-- Explicit vectors realizing the ideal 4-star: the center `e₀` and
`r e₀ + √(1-r²) u_i` with `u_i` the tetrahedral unit vectors in coordinates 1..3. -/
def starVec (r : ℝ) : Fin 5 → Sample 4
  | 0 => !₂[1, 0, 0, 0]
  | 1 => !₂[r, cs r, cs r, cs r]
  | 2 => !₂[r, cs r, -cs r, -cs r]
  | 3 => !₂[r, -cs r, cs r, -cs r]
  | 4 => !₂[r, -cs r, -cs r, cs r]

theorem inner_s4 (x y : Sample 4) :
    ⟪x, y⟫ = x 0 * y 0 + x 1 * y 1 + x 2 * y 2 + x 3 * y 3 := by
  simp [PiLp.inner_apply, Fin.sum_univ_four, mul_comm]

theorem gram_starVec {r : ℝ} (h0 : 0 ≤ r) (h1 : r ≤ 1) :
    gram (starVec r) = idealStarCov 4 r := by
  have hc : cs r ^ 2 = (1 - r ^ 2) / 3 := by
    unfold cs
    rw [Real.sq_sqrt]
    have : r ^ 2 ≤ 1 := by nlinarith
    linarith
  ext i j
  fin_cases i <;> fin_cases j <;>
    (simp only [gram]; rw [inner_s4]; simp [idealStarCov, starVec]) <;> norm_num <;>
      nlinarith [hc]

theorem measurableSet_frontierEvent (τ : ℝ) : MeasurableSet (frontierEvent 4 τ) := by
  have : frontierEvent 4 τ = {x : EuclideanSpace ℝ (Fin 5) | τ < x 0} ∩
      ⋃ i : Fin 5, ({x | i ≠ 0} ∩ {x | x i ≤ τ}) := by
    ext x; simp [frontierEvent]
  rw [this]
  refine (measurableSet_lt measurable_const (by fun_prop)).inter
    (MeasurableSet.iUnion fun i => ?_)
  exact (MeasurableSet.const _).inter (measurableSet_le (by fun_prop) measurable_const)

theorem idealStarFn_eq (hb : BridgeClaim) {r : ℝ} (h0 : 0 ≤ r) (h1 : r ≤ 1) (τ : ℝ) :
    idealStarFn 4 r τ = frontierProb (starVec r) τ := by
  unfold idealStarFn
  rw [← gram_starVec h0 h1, hb 4 5 (starVec r) _ (measurableSet_frontierEvent τ)]
  unfold frontierProb
  congr 1
  ext x
  simp only [frontierEvent, Set.mem_ofPred_eq]
  constructor
  · rintro ⟨h0, i, hi, hx⟩
    refine ⟨h0, ?_⟩
    induction i using Fin.cases with
    | zero => exact absurd rfl hi
    | succ i => exact ⟨i, hx⟩
  · rintro ⟨h0, i, hi⟩
    exact ⟨h0, i.succ, Fin.succ_ne_zero i, hi⟩

/-! ## Conditioning a star family on its center -/

/-- A star family: unit center `w 0`, unit leaves, center-leaf correlation `r`. -/
def IsStar {M : ℕ} (r : ℝ) (w : Fin 5 → Sample M) : Prop :=
  ‖w 0‖ = 1 ∧ (∀ i : Fin 4, ‖w i.succ‖ = 1) ∧ ∀ i : Fin 4, ⟪w 0, w i.succ⟫ = r

/-- The unit residual of each leaf of a star family. -/
def Ufam {M : ℕ} (r : ℝ) (w : Fin 5 → Sample M) (i : Fin 4) : Sample M :=
  (√(1 - r ^ 2))⁻¹ • (w i.succ - r • w 0)

/-- The conditional orthant probability at center value `z`. -/
def Fz (C : Matrix (Fin 4) (Fin 4) ℝ) (r τ z : ℝ) : ℝ :=
  if τ < z then gaussProb C (Aup 4 ((τ - r * z) / √(1 - r ^ 2))) else 0

section star

variable {M : ℕ} {r : ℝ} {w : Fin 5 → Sample M}

theorem sqrt_pos' (hr1 : r < 1) (hr : 0 < r) : 0 < √(1 - r ^ 2) :=
  Real.sqrt_pos.mpr (by nlinarith)

theorem sq_sqrt' (hr1 : r < 1) (hr : 0 < r) : √(1 - r ^ 2) ^ 2 = 1 - r ^ 2 :=
  Real.sq_sqrt (by nlinarith)

theorem inner_self_one {x : Sample M} (h : ‖x‖ = 1) : ⟪x, x⟫ = 1 := by
  rw [real_inner_self_eq_norm_sq, h]; norm_num

theorem inner_w0_Ufam (hw : IsStar r w) (i : Fin 4) : ⟪w 0, Ufam r w i⟫ = 0 := by
  obtain ⟨h0, -, hc⟩ := hw
  simp only [Ufam, real_inner_smul_right, inner_sub_right, hc, inner_self_one h0]
  ring

theorem inner_Ufam (hr : 0 < r) (hr1 : r < 1) (hw : IsStar r w) (i j : Fin 4) :
    ⟪Ufam r w i, Ufam r w j⟫ = (⟪w i.succ, w j.succ⟫ - r ^ 2) / (1 - r ^ 2) := by
  obtain ⟨h0, -, hc⟩ := hw
  have hs := sqrt_pos' hr1 hr
  have h2 := sq_sqrt' hr1 hr
  simp only [Ufam, real_inner_smul_left, real_inner_smul_right, inner_sub_left, inner_sub_right,
    inner_self_one h0, real_inner_comm (w 0) (w i.succ), hc]
  generalize √(1 - r ^ 2) = s at hs h2 ⊢
  rw [← h2]
  field_simp
  ring

theorem norm_Ufam (hr : 0 < r) (hr1 : r < 1) (hw : IsStar r w) (i : Fin 4) :
    ‖Ufam r w i‖ = 1 := by
  have h := inner_Ufam hr hr1 hw i i
  rw [inner_self_one (hw.2.1 i)] at h
  have h1 : ⟪Ufam r w i, Ufam r w i⟫ = 1 := by
    have hp : 0 < 1 - r ^ 2 := by nlinarith
    rw [h]; exact div_self hp.ne'
  rw [real_inner_self_eq_norm_sq] at h1
  nlinarith [norm_nonneg (Ufam r w i)]

theorem inner_succ_eq (hr : 0 < r) (hr1 : r < 1) (i : Fin 4) (x : Sample M) :
    ⟪w i.succ, x⟫ = r * ⟪w 0, x⟫ + √(1 - r ^ 2) * ⟪Ufam r w i, x⟫ := by
  have hs := sqrt_pos' hr1 hr
  simp only [Ufam, real_inner_smul_left, inner_sub_left]
  field_simp
  ring

end star

section cond

variable {M : ℕ} {r : ℝ} {w : Fin 5 → Sample M}

theorem allAbove_star (hb : BridgeClaim) (hr : 0 < r) (hr1 : r < 1) (hw : IsStar r w) (τ : ℝ) :
    allAbove w τ = ∫ z, Fz (gram (Ufam r w)) r τ z ∂gaussianReal 0 1 ∧
      Integrable (fun z => Fz (gram (Ufam r w)) r τ z) (gaussianReal 0 1) := by
  have hs := sqrt_pos' hr1 hr
  have hworth := inner_w0_Ufam hw
  set a := w 0 with ha_def
  have ha : ‖a‖ = 1 := hw.1
  set s := √(1 - r ^ 2) with hs_def
  set U := Ufam r w with hU
  set S : Set (ℝ × Sample M) :=
    {p | τ < p.1 ∧ ∀ i : Fin 4, τ < r * p.1 + s * ⟪U i, p.2⟫} with hS_def
  have hS : MeasurableSet S := by
    have : S = {p : ℝ × Sample M | τ < p.1} ∩
        ⋂ i : Fin 4, {p : ℝ × Sample M | τ < r * p.1 + s * ⟪U i, p.2⟫} := by
      ext p; simp [S]
    rw [this]
    exact (measurableSet_lt measurable_const measurable_fst).inter
      (MeasurableSet.iInter fun i => measurableSet_lt measurable_const (by fun_prop))
  have hev : {x : Sample M | ∀ j, τ < ⟪w j, x⟫} = {x | (⟪a, x⟫, C3.remainder a x) ∈ S} := by
    ext x
    simp only [Set.mem_ofPred_eq, S, C3.inner_remainder_of_orth (hworth _)]
    rw [Fin.forall_fin_succ]
    simp only [inner_succ_eq (w := w) hr hr1, ha_def, hU, hs_def]
  have hint : ∀ z, (stdGaussian (Sample M)).real {x | (z, C3.remainder a x) ∈ S} =
      Fz (gram U) r τ z := by
    intro z
    unfold Fz
    by_cases hz : τ < z
    · simp only [hz, ↓reduceIte]
      rw [hb M 4 U _ (measurableSet_Aup 4 _)]
      congr 1
      ext x
      simp only [S, Set.mem_ofPred_eq, Aup, hz, true_and, C3.inner_remainder_of_orth (hworth _)]
      change _ ↔ ∀ i, _ < ⟪U i, x⟫
      refine forall_congr' fun i => ?_
      rw [div_lt_iff₀ hs]
      constructor <;> intro h <;> linarith
    · simp only [hz, ↓reduceIte]
      have : {x : Sample M | (z, C3.remainder a x) ∈ S} = ∅ := by
        ext x; simp [S, hz]
      rw [this, measureReal_empty]
  have hfun : (fun z => Fz (gram U) r τ z) =
      fun z => (stdGaussian (Sample M)).real {x | (z, C3.remainder a x) ∈ S} := by
    funext z; exact (hint z).symm
  refine ⟨?_, ?_⟩
  · unfold allAbove
    rw [hev, C3.real_disintegrate ha hS]
    exact integral_congr_ae (ae_of_all _ fun z => hint z)
  · rw [hfun]
    exact C3.integrable_cond hS

end cond

/-! ## Orthant comparison from the simplex hypothesis -/

theorem allGe_eq_allAbove {M k : ℕ} (U : Fin k → Sample M) (hU : ∀ i, ‖U i‖ = 1) (θ : ℝ) :
    (stdGaussian (Sample M)).real {x | ∀ i, θ ≤ ⟪U i, x⟫} = allAbove U θ := by
  unfold allAbove
  set μ := stdGaussian (Sample M)
  have hN : μ.real (⋃ i, {x : Sample M | ⟪U i, x⟫ = θ}) = 0 := by
    refine le_antisymm ?_ measureReal_nonneg
    refine (measureReal_iUnion_fintype_le _).trans ?_
    exact (Finset.sum_eq_zero fun i _ => C3.prob_inner_eq (hU i) θ).le
  apply le_antisymm
  · calc μ.real {x | ∀ i, θ ≤ ⟪U i, x⟫}
        ≤ μ.real ({x | ∀ i, θ < ⟪U i, x⟫} ∪ ⋃ i, {x | ⟪U i, x⟫ = θ}) := by
          refine measureReal_mono ?_ (measure_ne_top _ _)
          intro x hx
          by_cases h : ∀ i, θ < ⟪U i, x⟫
          · exact Or.inl h
          · right
            push Not at h
            obtain ⟨i, hi⟩ := h
            exact Set.mem_iUnion.mpr ⟨i, le_antisymm hi (hx i)⟩
      _ ≤ μ.real {x | ∀ i, θ < ⟪U i, x⟫} + μ.real (⋃ i, {x : Sample M | ⟪U i, x⟫ = θ}) :=
          measureReal_union_le _ _
      _ = μ.real {x | ∀ i, θ < ⟪U i, x⟫} := by rw [hN, add_zero]
  · exact measureReal_mono (fun x hx i => (hx i).le) (measure_ne_top _ _)

theorem measurableSet_lowerOrthant (k : ℕ) (t : ℝ) :
    MeasurableSet (lowerOrthant (m := k) t) := by
  have : lowerOrthant (m := k) t = ⋂ i, {y : EuclideanSpace ℝ (Fin k) | y i ≤ t} := by
    ext y; simp [lowerOrthant]
  rw [this]
  exact MeasurableSet.iInter fun i => measurableSet_le (by fun_prop) measurable_const

theorem upper_eq_lower (hb : BridgeClaim) {M k : ℕ} (U : Fin k → Sample M)
    (hU : ∀ i, ‖U i‖ = 1) (θ : ℝ) :
    gaussProb (gram U) (Aup k θ) = gaussProb (gram U) (lowerOrthant (-θ)) := by
  have hg : gram (fun i => -U i) = gram U := by ext i j; simp [gram]
  rw [← allAbove_eq_gaussProb hb, ← hg,
    hb M k (fun i => -U i) _ (measurableSet_lowerOrthant k _), ← allGe_eq_allAbove U hU θ]
  congr 1
  ext x
  simp [lowerOrthant, inner_neg_left]

theorem upper_cmp (hb : BridgeClaim) (h1 : SimplexFourClaim) {M M' : ℕ}
    (U : Fin 4 → Sample M) (V : Fin 4 → Sample M') (hU : ∀ i, ‖U i‖ = 1)
    (hV : ∀ i, ‖V i‖ = 1) (hVg : gram V = regularSimplexGram 4) (θ : ℝ) :
    gaussProb (gram V) (Aup 4 θ) ≤ gaussProb (gram U) (Aup 4 θ) := by
  rw [upper_eq_lower hb U hU, upper_eq_lower hb V hV, hVg]
  refine h1 (gram U) ⟨Matrix.posSemidef_gram ℝ U, fun i => ?_⟩ (-θ)
  exact inner_self_one (hU i)

/-- Comparison of two star families through their conditional orthants. -/
theorem allAbove_star_le (hb : BridgeClaim) (h1 : SimplexFourClaim) {M M' : ℕ} {r : ℝ}
    (hr : 0 < r) (hr1 : r < 1) {w : Fin 5 → Sample M} {w' : Fin 5 → Sample M'}
    (hw : IsStar r w) (hw' : IsStar r w') (hg : gram (Ufam r w') = regularSimplexGram 4)
    (τ : ℝ) : allAbove w' τ ≤ allAbove w τ := by
  obtain ⟨e1, i1⟩ := allAbove_star hb hr hr1 hw τ
  obtain ⟨e2, i2⟩ := allAbove_star hb hr hr1 hw' τ
  rw [e1, e2]
  refine integral_mono i2 i1 fun z => ?_
  unfold Fz
  split_ifs with hz
  · exact upper_cmp hb h1 (Ufam r w) (Ufam r w') (norm_Ufam hr hr1 hw) (norm_Ufam hr hr1 hw')
      hg _
  · exact le_rfl

/-! ## The ideal star is a star family with simplex residuals -/

theorem norm_eq_one_of_inner_self {M : ℕ} {x : Sample M} (h : ⟪x, x⟫ = 1) : ‖x‖ = 1 := by
  rw [real_inner_self_eq_norm_sq] at h
  nlinarith [norm_nonneg x]

theorem inner_starVec {r : ℝ} (h0 : 0 ≤ r) (h1 : r ≤ 1) (i j : Fin 5) :
    ⟪starVec r i, starVec r j⟫ = idealStarCov 4 r i j :=
  congrFun (congrFun (gram_starVec h0 h1) i) j

theorem isStar_starVec {r : ℝ} (hr : 0 < r) (hr1 : r < 1) : IsStar r (starVec r) := by
  have h0 := hr.le
  have h1 := hr1.le
  refine ⟨norm_eq_one_of_inner_self ?_, fun i => norm_eq_one_of_inner_self ?_, fun i => ?_⟩
  · rw [inner_starVec h0 h1]; simp [idealStarCov]
  · rw [inner_starVec h0 h1]; simp [idealStarCov]
  · rw [inner_starVec h0 h1]
    simp [idealStarCov, (Fin.succ_ne_zero i).symm]

theorem gram_Ufam_starVec {r : ℝ} (hr : 0 < r) (hr1 : r < 1) :
    gram (Ufam r (starVec r)) = regularSimplexGram 4 := by
  have hp : (1 - r ^ 2) ≠ 0 := by nlinarith
  ext i j
  show ⟪Ufam r (starVec r) i, Ufam r (starVec r) j⟫ = _
  rw [inner_Ufam hr hr1 (isStar_starVec hr hr1), inner_starVec hr.le hr1.le]
  by_cases hij : i = j
  · subst hij
    simp [idealStarCov, regularSimplexGram]
    rw [div_self hp]; norm_num
  · have : i.succ ≠ j.succ := fun h => hij (Fin.succ_injective _ h)
    simp [idealStarCov, regularSimplexGram, hij, this, Fin.succ_ne_zero]
    field_simp

/-! ## Zero-padding and fresh orthonormal directions -/

/-- Zero-pad into `K` extra coordinates. -/
def pad {m : ℕ} (K : ℕ) (x : Sample m) : Sample (m + K) :=
  WithLp.toLp 2 (Fin.append x.ofLp (0 : Fin K → ℝ))

/-- The `j`-th fresh coordinate direction. -/
def xi (m K : ℕ) (j : Fin K) : Sample (m + K) :=
  EuclideanSpace.single (Fin.natAdd m j) 1

theorem inner_pad {m : ℕ} (K : ℕ) (x y : Sample m) : ⟪pad K x, pad K y⟫ = ⟪x, y⟫ := by
  simp only [PiLp.inner_apply, pad, Fin.sum_univ_add, Fin.append_left, Fin.append_right]
  simp

theorem inner_pad_xi {m : ℕ} (K : ℕ) (x : Sample m) (j : Fin K) :
    ⟪pad K x, xi m K j⟫ = 0 := by
  simp [xi, EuclideanSpace.inner_single_right, pad]

theorem inner_xi_pad {m : ℕ} (K : ℕ) (x : Sample m) (j : Fin K) :
    ⟪xi m K j, pad K x⟫ = 0 := by
  rw [real_inner_comm, inner_pad_xi]

theorem inner_xi_xi (m K : ℕ) (i j : Fin K) :
    ⟪xi m K i, xi m K j⟫ = if i = j then 1 else 0 := by
  simp [xi, EuclideanSpace.inner_single_right]
  exact if_congr eq_comm rfl rfl

theorem gram_pad {m k : ℕ} (K : ℕ) (w : Fin k → Sample m) :
    gram (fun i => pad K (w i)) = gram w := by
  ext i j; exact inner_pad K _ _

/-! ## Lowering the incident correlations to `r` -/

/-- The comparison star `c`: center `A`, leaves `B_i - (ρ_i - r) A + √(ρ_i² - r²) ξ_i`. -/
def cfam {N : ℕ} (r : ℝ) (A : Sample N) (B ξ : Fin 4 → Sample N) (ρ : Fin 4 → ℝ) :
    Fin 5 → Sample N :=
  Fin.cons A (fun i => B i - (ρ i - r) • A + √(ρ i ^ 2 - r ^ 2) • ξ i)

section cfam

variable {N : ℕ} {r : ℝ} {A : Sample N} {B ξ : Fin 4 → Sample N} {ρ : Fin 4 → ℝ}

theorem inner_cfam_succ (hA : ⟪A, A⟫ = 1) (hAB : ∀ i, ⟪A, B i⟫ = ρ i)
    (hAξ : ∀ i, ⟪A, ξ i⟫ = 0) (hBξ : ∀ i j, ⟪B i, ξ j⟫ = 0)
    (hξ : ∀ i j, ⟪ξ i, ξ j⟫ = if i = j then 1 else 0) (i j : Fin 4) :
    ⟪cfam r A B ξ ρ i.succ, cfam r A B ξ ρ j.succ⟫ =
      ⟪B i, B j⟫ - ρ i * ρ j + r ^ 2 +
        (if i = j then √(ρ i ^ 2 - r ^ 2) * √(ρ j ^ 2 - r ^ 2) else 0) := by
  have e1 : ⟪B i, A⟫ = ρ i := by rw [real_inner_comm, hAB]
  have e2 : ⟪ξ i, A⟫ = 0 := by rw [real_inner_comm, hAξ]
  have e3 : ⟪ξ i, B j⟫ = 0 := by rw [real_inner_comm, hBξ]
  simp only [cfam, Fin.cons_succ, inner_add_left, inner_add_right, inner_sub_left, inner_sub_right,
    real_inner_smul_left, real_inner_smul_right, hA, hAB, hAξ, hBξ, hξ, e1, e2, e3]
  by_cases h : i = j
  · subst h; simp; ring
  · simp [h]; ring

theorem inner_cfam_zero_succ (hA : ⟪A, A⟫ = 1) (hAB : ∀ i, ⟪A, B i⟫ = ρ i)
    (hAξ : ∀ i, ⟪A, ξ i⟫ = 0) (i : Fin 4) :
    ⟪cfam r A B ξ ρ 0, cfam r A B ξ ρ i.succ⟫ = r := by
  simp only [cfam, Fin.cons_zero, Fin.cons_succ, inner_add_right, inner_sub_right,
    real_inner_smul_right, hA, hAB, hAξ]
  ring

theorem isStar_cfam (hr : 0 < r) (hA : ⟪A, A⟫ = 1) (hAB : ∀ i, ⟪A, B i⟫ = ρ i)
    (hAξ : ∀ i, ⟪A, ξ i⟫ = 0) (hBξ : ∀ i j, ⟪B i, ξ j⟫ = 0)
    (hξ : ∀ i j, ⟪ξ i, ξ j⟫ = if i = j then 1 else 0) (hBn : ∀ i, ⟪B i, B i⟫ = 1)
    (hρ : ∀ i, r ≤ ρ i) : IsStar r (cfam r A B ξ ρ) := by
  refine ⟨?_, fun i => ?_, fun i => ?_⟩
  · apply norm_eq_one_of_inner_self
    simpa [cfam] using hA
  · apply norm_eq_one_of_inner_self
    rw [inner_cfam_succ hA hAB hAξ hBξ hξ, hBn]
    simp only [↓reduceIte]
    have : 0 ≤ ρ i ^ 2 - r ^ 2 := by nlinarith [hρ i]
    rw [Real.mul_self_sqrt this]
    ring
  · exact inner_cfam_zero_succ hA hAB hAξ i

end cfam

/-! ## The core inequality for exactly four neighbors -/

theorem core4 (hsl : SlepianClaim) (hb : BridgeClaim) (h1 : SimplexFourClaim)
    {m : ℕ} {r : ℝ} (hr : 0 < r) (hr1 : r < 1) (a : Sample m) (b : Fin 4 → Sample m)
    (ha : ‖a‖ = 1) (hbn : ∀ i, ‖b i‖ = 1) (hab : ∀ i, r ≤ ⟪a, b i⟫) (τ : ℝ) :
    allAbove (starVec r) τ ≤ allAbove (Fin.cons a b : Fin 5 → Sample m) τ := by
  set A : Sample (m + 4) := pad 4 a with hA_def
  set B : Fin 4 → Sample (m + 4) := fun i => pad 4 (b i) with hB_def
  set ρ : Fin 4 → ℝ := fun i => ⟪a, b i⟫ with hρ_def
  set ξ : Fin 4 → Sample (m + 4) := xi m 4 with hξ_def
  have haa : ⟪A, A⟫ = 1 := by rw [hA_def, inner_pad]; exact inner_self_one ha
  have hAB : ∀ i, ⟪A, B i⟫ = ρ i := fun i => inner_pad 4 a (b i)
  have hAξ : ∀ i, ⟪A, ξ i⟫ = 0 := fun i => inner_pad_xi 4 a i
  have hBξ : ∀ i j, ⟪B i, ξ j⟫ = 0 := fun i j => inner_pad_xi 4 (b i) j
  have hξξ : ∀ i j, ⟪ξ i, ξ j⟫ = if i = j then 1 else 0 := fun i j => inner_xi_xi m 4 i j
  have hBn : ∀ i, ⟪B i, B i⟫ = 1 := fun i => by
    rw [hB_def]; simp only; rw [inner_pad]; exact inner_self_one (hbn i)
  have hstar : IsStar r (cfam r A B ξ ρ) :=
    isStar_cfam hr haa hAB hAξ hBξ hξξ hBn hab
  -- the ideal star is below the comparison star
  have step1 : allAbove (starVec r) τ ≤ allAbove (cfam r A B ξ ρ) τ :=
    allAbove_star_le hb h1 hr hr1 hstar (isStar_starVec hr hr1) (gram_Ufam_starVec hr hr1) τ
  -- Slepian
  have hA1 : ‖A‖ = 1 := norm_eq_one_of_inner_self haa
  have step2 : allAbove (cfam r A B ξ ρ) τ ≤
      allAbove (Fin.cons A B : Fin 5 → Sample (m + 4)) τ := by
    have hn1 : ∀ i, ‖cfam r A B ξ ρ i‖ = 1 := by
      intro i
      induction i using Fin.cases with
      | zero => exact hstar.1
      | succ i => exact hstar.2.1 i
    have hn2 : ∀ i, ‖(Fin.cons A B : Fin 5 → Sample (m + 4)) i‖ = 1 := by
      intro i
      induction i using Fin.cases with
      | zero => simpa using hA1
      | succ i => simpa using norm_eq_one_of_inner_self (hBn i)
    have hle : ∀ i j, i ≠ j → ⟪cfam r A B ξ ρ i, cfam r A B ξ ρ j⟫ ≤
        ⟪(Fin.cons A B : Fin 5 → Sample (m + 4)) i, (Fin.cons A B : Fin 5 → Sample (m + 4)) j⟫ := by
      intro i j hij
      induction i using Fin.cases with
      | zero =>
        induction j using Fin.cases with
        | zero => exact absurd rfl hij
        | succ j =>
          rw [inner_cfam_zero_succ haa hAB hAξ]
          simpa [hAB] using hab j
      | succ i =>
        induction j using Fin.cases with
        | zero =>
          rw [real_inner_comm, inner_cfam_zero_succ haa hAB hAξ, real_inner_comm]
          simpa [hAB] using hab i
        | succ j =>
          have hij' : i ≠ j := fun h => hij (h ▸ rfl)
          rw [inner_cfam_succ haa hAB hAξ hBξ hξξ]
          simp only [hij', ↓reduceIte, add_zero, Fin.cons_succ]
          have : r ^ 2 ≤ ρ i * ρ j := by
            have := hab i
            have := hab j
            nlinarith
          simp only [hB_def, inner_pad]
          nlinarith
    exact hsl (m + 4) 5 (cfam r A B ξ ρ) (Fin.cons A B) (fun i => by rw [hn1, hn2]) hle
      (fun _ => τ)
  -- transfer back to `Sample m`
  have step3 : allAbove (Fin.cons A B : Fin 5 → Sample (m + 4)) τ =
      allAbove (Fin.cons a b : Fin 5 → Sample m) τ := by
    have hfam : (Fin.cons A B : Fin 5 → Sample (m + 4)) =
        fun i => pad 4 ((Fin.cons a b : Fin 5 → Sample m) i) := by
      funext i
      induction i using Fin.cases with
      | zero => simp [hA_def]
      | succ i => simp [hB_def]
    rw [hfam]
    exact allAbove_eq_of_gram_eq hb _ _ (gram_pad 4 _) τ
  exact (step1.trans step2).trans step3.le

/-! ## The main theorem -/

theorem prob_pad (hb : BridgeClaim) {m k : ℕ} (w : Fin k → Sample m) (K : ℕ)
    (A : Set (EuclideanSpace ℝ (Fin k))) (hA : MeasurableSet A) :
    (stdGaussian (Sample m)).real {x | iv w x ∈ A} =
      (stdGaussian (Sample (m + K))).real {x | iv (fun i => pad K (w i)) x ∈ A} := by
  have h1 := hb m k w A hA
  have h2 := hb (m + K) k (fun i => pad K (w i)) A hA
  rw [gram_pad] at h2
  exact h1.symm.trans h2

/-- The ideal-star center-above probability is `Φ(-τ)`, independently of the space. -/
theorem starVec_center (hb : BridgeClaim) {r : ℝ} (hr : 0 < r) (hr1 : r < 1) (τ : ℝ) :
    idealStarFn 4 r τ = (gaussianReal 0 1).real (Set.Ioi τ) - allAbove (starVec r) τ := by
  rw [idealStarFn_eq hb hr.le hr1.le, frontierProb_eq,
    prob_unit_gt (isStar_starVec hr hr1).1]

end StarFour

open GaussianField (Sample) in
open StarFour in

/-- **The 4-star comparison**, conditional on Slepian, the Gaussian bridge and the
regular-simplex hypothesis. -/
theorem starFour_of (hsl : SlepianClaim) (hb : BridgeClaim) (h1 : SimplexFourClaim) :
    StarFourClaim := by
  intro m r hr hr1 a d hd b ha hbn hab τ
  set A : Sample (m + 4) := pad 4 a with hA_def
  have haa : ⟪A, A⟫ = 1 := by rw [hA_def, inner_pad]; exact inner_self_one ha
  have hA1 : ‖A‖ = 1 := norm_eq_one_of_inner_self haa
  have hAξ : ∀ i, ⟪A, xi m 4 i⟫ = 0 := fun i => inner_pad_xi 4 a i
  have hξA : ∀ i, ⟪xi m 4 i, A⟫ = 0 := fun i => inner_xi_pad 4 a i
  have hs2 := sq_sqrt' hr1 hr
  set b' : Fin 4 → Sample (m + 4) := fun i =>
    if h : i.val < d then pad 4 (b ⟨i.val, h⟩) else r • A + √(1 - r ^ 2) • xi m 4 i with hb'_def
  have hb'n : ∀ i, ‖b' i‖ = 1 := by
    intro i
    by_cases h : i.val < d
    · simp only [hb'_def, h, ↓reduceDIte]
      exact norm_eq_one_of_inner_self (by rw [inner_pad]; exact inner_self_one (hbn _))
    · simp only [hb'_def, h, ↓reduceDIte]
      apply norm_eq_one_of_inner_self
      simp only [inner_add_left, inner_add_right, real_inner_smul_left, real_inner_smul_right, haa,
        hAξ, hξA, inner_xi_xi, ↓reduceIte]
      nlinarith
  have hb'ab : ∀ i, r ≤ ⟪A, b' i⟫ := by
    intro i
    by_cases h : i.val < d
    · simp only [hb'_def, h, ↓reduceDIte, hA_def, inner_pad]
      exact hab _
    · simp only [hb'_def, h, ↓reduceDIte]
      simp only [inner_add_right, real_inner_smul_right, haa, hAξ]
      simp
  have core := core4 hsl hb h1 hr hr1 A b' hA1 hb'n hb'ab τ
  -- left side as a probability in `Sample (m + 4)`
  set Aev : Set (EuclideanSpace ℝ (Fin (d + 1))) :=
    {y | τ < y 0 ∧ ∃ i : Fin d, y i.succ ≤ τ} with hAev
  have hAev_meas : MeasurableSet Aev := by
    have : Aev = {y : EuclideanSpace ℝ (Fin (d + 1)) | τ < y 0} ∩
        ⋃ i : Fin d, {y | y i.succ ≤ τ} := by
      ext y; simp [Aev]
    rw [this]
    exact (measurableSet_lt measurable_const (by fun_prop)).inter
      (MeasurableSet.iUnion fun i => measurableSet_le (by fun_prop) measurable_const)
  have hL : (stdGaussian (Sample m)).real {x | τ < ⟪a, x⟫ ∧ ∃ i, ⟪b i, x⟫ ≤ τ} =
      (stdGaussian (Sample (m + 4))).real {x | τ < ⟪A, x⟫ ∧ ∃ i : Fin d, ⟪pad 4 (b i), x⟫ ≤ τ} := by
    have := prob_pad hb (Fin.cons a b : Fin (d + 1) → Sample m) 4 Aev hAev_meas
    have e1 : {x : Sample m | τ < ⟪a, x⟫ ∧ ∃ i, ⟪b i, x⟫ ≤ τ} =
        {x | iv (Fin.cons a b : Fin (d + 1) → Sample m) x ∈ Aev} := by
      ext x; simp [iv, Aev]
    have e2 : {x : Sample (m + 4) | τ < ⟪A, x⟫ ∧ ∃ i : Fin d, ⟪pad 4 (b i), x⟫ ≤ τ} =
        {x | iv (fun i => pad 4 ((Fin.cons a b : Fin (d + 1) → Sample m) i)) x ∈ Aev} := by
      ext x; simp [iv, Aev, hA_def]
    rw [e1, e2]
    exact this
  have hsub : {x : Sample (m + 4) | τ < ⟪A, x⟫ ∧ ∃ i : Fin d, ⟪pad 4 (b i), x⟫ ≤ τ} ⊆
      {x | τ < ⟪(Fin.cons A b' : Fin 5 → Sample (m + 4)) 0, x⟫ ∧
        ∃ i : Fin 4, ⟪(Fin.cons A b' : Fin 5 → Sample (m + 4)) i.succ, x⟫ ≤ τ} := by
    rintro x ⟨h0, i, hi⟩
    refine ⟨by simpa using h0, ⟨i.val, by omega⟩, ?_⟩
    have hbi : b' ⟨i.val, by omega⟩ = pad 4 (b i) := by simp [hb'_def, i.isLt]
    rw [Fin.cons_succ, hbi]
    exact hi
  calc (stdGaussian (Sample m)).real {x | τ < ⟪a, x⟫ ∧ ∃ i, ⟪b i, x⟫ ≤ τ}
      = _ := hL
    _ ≤ frontierProb (Fin.cons A b' : Fin 5 → Sample (m + 4)) τ :=
        measureReal_mono hsub (measure_ne_top _ _)
    _ = (gaussianReal 0 1).real (Set.Ioi τ) - allAbove (Fin.cons A b' : Fin 5 → Sample (m + 4)) τ := by
        rw [frontierProb_eq]
        simp only [Fin.cons_zero]
        rw [prob_unit_gt hA1]
    _ ≤ idealStarFn 4 r τ := by
        rw [starVec_center hb hr hr1]
        linarith

end T2
end CubicPathwidth

#print axioms CubicPathwidth.T2.starFour_of