import RNLAPhaseB.B3Base
import RNLAPhaseB.Complexity

open MeasureTheory
open scoped Matrix
open scoped Matrix.Norms.L2Operator
open scoped BigOperators
open scoped ENNReal

namespace RNLAPhaseB.HMT2011.PropBidiag

/-- The source's `m × n` lower-bidiagonal matrix `L` (0-based indices).
Row `0` holds `X_n` at column `0`. For `1 ≤ i ≤ m-1`, row `i` holds `Y_{m-i}` at
column `i-1` and `X_{n-i}` at column `i`. Every other entry is zero.
Here `x i` is the diagonal entry of row `i`, which is the source's `X_{n-i}`.
`y k` is the subdiagonal entry at `(k+1, k)`, which is the source's `Y_{m-1-k}`. -/
def bidiagonalL (m n : ℕ) (x : Fin m → ℝ) (y : Fin (m - 1) → ℝ) :
    Matrix (Fin m) (Fin n) ℝ :=
  fun i j =>
    if (j : ℕ) = (i : ℕ) then x i
    else if h : (j : ℕ) + 1 = (i : ℕ) then
      y ⟨(j : ℕ), by have := i.isLt; omega⟩
    else 0

/-- HMT 2011, Proposition (prop:bidiag).
Let `n ≥ m` and let `G` be an `m × n` standard Gaussian matrix.
Then `G` is orthogonally equivalent to the random bidiagonal matrix `L`, that is,
`G = U L Vᵀ` with `U` and `V` orthogonal.
For each `j`, `X_j^2` and `Y_j^2` are `χ²` with `j` degrees of freedom.
The `χ²_j` law is written as `Gamma(j/2, rate 1/2)`.
In the Lean indexing, `x i` is `X_{n-i}` with `n - i` degrees of freedom.
`y k` is `Y_{m-1-k}` with `m - 1 - k` degrees of freedom.
All the variates `X`, `Y` are mutually independent. -/
theorem prop_bidiag {Ω : Type*} [MeasurableSpace Ω] (P : Measure Ω)
    (m n : ℕ) (hmn : m ≤ n)
    (G : Ω → Matrix (Fin m) (Fin n) ℝ)
    (hG : RNLAPhaseB.StandardGaussianTestMatrix P G) :
    ∃ (U : Ω → Matrix (Fin m) (Fin m) ℝ) (V : Ω → Matrix (Fin n) (Fin n) ℝ)
      (X : Fin m → Ω → ℝ) (Y : Fin (m - 1) → Ω → ℝ),
      (∀ ω, (U ω)ᵀ * U ω = 1 ∧ (V ω)ᵀ * V ω = 1 ∧
        G ω = U ω * bidiagonalL m n (fun i => X i ω) (fun k => Y k ω) * (V ω)ᵀ) ∧
      (∀ i : Fin m, ProbabilityTheory.HasLaw (fun ω => (X i ω) ^ 2)
        (ProbabilityTheory.gammaMeasure (((n - (i : ℕ) : ℕ) : ℝ) / 2) (1 / 2)) P) ∧
      (∀ k : Fin (m - 1), ProbabilityTheory.HasLaw (fun ω => (Y k ω) ^ 2)
        (ProbabilityTheory.gammaMeasure (((m - 1 - (k : ℕ) : ℕ) : ℝ) / 2) (1 / 2)) P) ∧
      ProbabilityTheory.iIndepFun (Sum.elim X Y) P := by
  sorry

end RNLAPhaseB.HMT2011.PropBidiag