/-
Copyright (c) 2026 Michael R. Douglas. All rights reserved.
Released under Apache 2.0 license as described in the file LICENSE.

# Gaussian Measure Construction

Given a CLM T : E → H from a nuclear Fréchet space to a separable
Hilbert space (finite- or infinite-dimensional), constructs a centered
Gaussian probability measure μ on the weak dual E' with covariance
C(f,g) = ⟨T(f), T(g)⟩_H.

## Construction overview

1. **Nuclear factorization** (TargetFactorization.lean):
   Extract adapted ONB {eₙ}, intermediate space K, CLM j, vectors vₙ
   with ∑ ‖vₙ‖ < ∞.

2. **Noise measure**: μ_noise = ⊗_n N(0,1) on ℝ^ℕ

3. **Series limit map**: For a.e. noise ξ, the series
   ω(f) = ∑ₙ ξₙ ⟨eₙ, T(f)⟩ = ⟨∑ₙ ξₙ vₙ, j(f)⟩_K
   converges and defines ω ∈ E' (continuous by CLM composition).

4. **Pushforward**: μ = (seriesLimit)_* μ_noise

## Main definitions

- `Configuration E` — the configuration space E' = WeakDual ℝ E
- `covariance T f g` — the covariance bilinear form C(f,g) = ⟨T(f), T(g)⟩_H
- `measure T` — the constructed probability measure

## Main theorems

- `charFun` — characteristic functional identity
- See `Properties.lean` for moments and integrability.
-/

import GaussianField.TargetFactorization
import Mathlib.Probability.Distributions.Gaussian.Real
import Mathlib.Probability.Distributions.Gaussian.Basic
import Mathlib.Probability.ProductMeasure
import Mathlib.Probability.Independence.InfinitePi
import Mathlib.Probability.Independence.Integration
import Mathlib.MeasureTheory.Measure.ProbabilityMeasure
import Mathlib.Analysis.Fourier.BoundedContinuousFunctionChar
import Mathlib.Analysis.Fourier.FourierTransform
import Mathlib.Analysis.InnerProductSpace.Dual
import Mathlib.Analysis.InnerProductSpace.ProdL2
import Mathlib.Analysis.Normed.Lp.MeasurableSpace
import Mathlib.MeasureTheory.Group.IntegralConvolution
import Mathlib.MeasureTheory.Integral.Pi
import Mathlib.MeasureTheory.Measure.FiniteMeasureExt
import Mathlib.Topology.Algebra.Module.Spaces.WeakDual

noncomputable section

namespace GaussianField

open MeasureTheory ProbabilityTheory TopologicalSpace
open scoped BigOperators
open Classical

variable {E : Type*} [AddCommGroup E] [Module ℝ E]
  [TopologicalSpace E] [IsTopologicalAddGroup E] [ContinuousSMul ℝ E]
  [DyninMityaginSpace E]

/-- Configuration space: continuous linear functionals on E. -/
abbrev Configuration (E : Type*)
    [AddCommGroup E] [Module ℝ E]
    [TopologicalSpace E] [IsTopologicalAddGroup E] [ContinuousSMul ℝ E] :=
  WeakDual ℝ E

/-- Cylindrical σ-algebra on Configuration E: generated by evaluation maps
    ω ↦ ω(f) for each test function f. -/
instance instMeasurableSpaceConfiguration :
    MeasurableSpace (Configuration E) :=
  MeasurableSpace.comap (fun ω f => ω f) MeasurableSpace.pi

omit [DyninMityaginSpace E] in
/-- A function into Configuration is measurable iff all evaluations are measurable. -/
theorem configuration_measurable_of_eval_measurable
    {X : Type*} [MeasurableSpace X] (g : X → Configuration E)
    (h : ∀ φ : E, Measurable (fun x => g x φ)) :
    @Measurable X (Configuration E) _ instMeasurableSpaceConfiguration g := by
  rw [measurable_iff_comap_le]
  show (MeasurableSpace.comap (fun ω f => ω f) MeasurableSpace.pi).comap g ≤ _
  rw [MeasurableSpace.comap_comp]
  exact (measurable_pi_lambda _ h).comap_le

omit [DyninMityaginSpace E] in
/-- Each evaluation map ω ↦ ω(f) is measurable w.r.t. the cylindrical σ-algebra. -/
theorem configuration_eval_measurable (φ : E) :
    @Measurable (Configuration E) ℝ instMeasurableSpaceConfiguration _ (fun ω => ω φ) := by
  intro s hs
  exact ⟨(fun f => f φ) ⁻¹' s, measurable_pi_apply φ hs, rfl⟩

variable {H : Type*} [NormedAddCommGroup H] [InnerProductSpace ℝ H]
  [CompleteSpace H] [SeparableSpace H]

/-- The covariance bilinear form induced by T: C(f,g) = ⟨T(f), T(g)⟩_H. -/
def covariance (T : E →L[ℝ] H) (f g : E) : ℝ :=
  @inner ℝ H _ (T f) (T g)

/-! ## Factorization Extraction

We extract the factorization data from `nuclear_clm_target_factorization` using
noncomputable choice. This gives us the adapted ONB, intermediate space K, CLM j,
and nuclear vectors v. -/

set_option linter.defProp false in
/-- The full factorization result, privately stored. -/
private noncomputable def fullFactorization
    (T : E →L[ℝ] H) (h_inf : ¬ FiniteDimensional ℝ H) :=
  nuclear_clm_target_factorization h_inf T

set_option linter.defProp false in
private noncomputable def factProps (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) :=
  (fullFactorization T h_inf).choose_spec.choose_spec.choose_spec.choose_spec.choose_spec.choose_spec.choose_spec

/-- The adapted ONB from the nuclear factorization. -/
private noncomputable def adaptedONB (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) : ℕ → H :=
  (fullFactorization T h_inf).choose

/-- The intermediate Hilbert space type K from the factorization. -/
private noncomputable def FactK (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) : Type :=
  (fullFactorization T h_inf).choose_spec.choose

private noncomputable instance instNACG_FactK (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) : NormedAddCommGroup (FactK T h_inf) :=
  (fullFactorization T h_inf).choose_spec.choose_spec.choose

private noncomputable instance instIPS_FactK (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) : InnerProductSpace ℝ (FactK T h_inf) :=
  (fullFactorization T h_inf).choose_spec.choose_spec.choose_spec.choose

private noncomputable instance instCS_FactK (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) : CompleteSpace (FactK T h_inf) :=
  (fullFactorization T h_inf).choose_spec.choose_spec.choose_spec.choose_spec.choose

/-- The CLM j : E → K from the factorization. -/
private noncomputable def factJ (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) : E →L[ℝ] FactK T h_inf :=
  (fullFactorization T h_inf).choose_spec.choose_spec.choose_spec.choose_spec.choose_spec.choose

/-- The vectors v_n ∈ K from the factorization. -/
private noncomputable def factV (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) : ℕ → FactK T h_inf :=
  (fullFactorization T h_inf).choose_spec.choose_spec.choose_spec.choose_spec.choose_spec.choose_spec.choose

private lemma adaptedONB_orthonormal (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) : Orthonormal ℝ (adaptedONB T h_inf) :=
  (factProps T h_inf).1

private lemma adaptedONB_spans (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) :
    ⊤ ≤ (Submodule.span ℝ (Set.range (adaptedONB T h_inf))).topologicalClosure :=
  (factProps T h_inf).2.1

private lemma factV_summable (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) : Summable (fun n => ‖factV T h_inf n‖) :=
  (factProps T h_inf).2.2.1

private lemma factV_inner_eq (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) (f : E) (n : ℕ) :
    @inner ℝ H _ (adaptedONB T h_inf n) (T f) =
    @inner ℝ (FactK T h_inf) _ (factV T h_inf n) (factJ T h_inf f) :=
  (factProps T h_inf).2.2.2 f n

/-- The adapted HilbertBasis for real Parseval identity. -/
private noncomputable def adaptedBasis (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) : HilbertBasis ℕ ℝ H :=
  HilbertBasis.mk (adaptedONB_orthonormal T h_inf) (adaptedONB_spans T h_inf)

private lemma adaptedBasis_eq_ONB (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) :
    ⇑(adaptedBasis T h_inf) = adaptedONB T h_inf :=
  HilbertBasis.coe_mk _ _

/-- The coefficient of T f in the n-th ONB direction. -/
private noncomputable def coeff (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) (n : ℕ) (f : E) : ℝ :=
  @inner ℝ H _ (adaptedONB T h_inf n) (T f)

/-- Summability of squared coefficients (by Parseval). -/
private lemma coeff_sq_summable (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) (f : E) :
    Summable (fun n => (coeff T h_inf n f) ^ 2) := by
  set c := (adaptedBasis T h_inf).repr (T f) with hc_def
  have hsumm := Memℓp.summable (by norm_num : (0 : ℝ) < (2 : ENNReal).toReal) c.property
  simp only [ENNReal.toReal_ofNat] at hsumm
  have hsumm' : Summable (fun n => ‖(c : ℕ → ℝ) n‖^(2 : ℕ)) := by
    simp only [← Real.rpow_natCast] at hsumm ⊢
    convert hsumm using 2; rfl
  convert hsumm' using 2 with n
  rw [coeff, show adaptedONB T h_inf n = adaptedBasis T h_inf n from
    (adaptedBasis_eq_ONB T h_inf).symm ▸ rfl]
  rw [show @inner ℝ H _ (adaptedBasis T h_inf n) (T f) =
    (adaptedBasis T h_inf).repr (T f) n from
    ((adaptedBasis T h_inf).repr_apply_apply (T f) n).symm]
  simp [sq_abs, ← hc_def]

/-- Parseval identity: sum of squared coefficients = ‖T f‖². -/
private lemma coeff_parseval (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) (f : E) :
    ∑' n, (coeff T h_inf n f) ^ 2 = ‖T f‖ ^ 2 := by
  set x := T f with hx_def
  set c := (adaptedBasis T h_inf).repr x with hc_def
  have hiso : ‖x‖ = ‖c‖ := ((adaptedBasis T h_inf).repr.norm_map x).symm
  have hlp := lp.norm_rpow_eq_tsum (p := 2) (by norm_num : (0 : ℝ) < (2 : ENNReal).toReal) c
  simp only [ENNReal.toReal_ofNat] at hlp
  have hlp' : ‖c‖^(2 : ℕ) = ∑' n, ‖(c : ℕ → ℝ) n‖^(2 : ℕ) := by
    simp only [← Real.rpow_natCast] at hlp ⊢
    convert hlp using 2 <;> rfl
  calc ∑' n, (coeff T h_inf n f)^2
      = ∑' n, ‖(c : ℕ → ℝ) n‖^2 := by
        congr 1; ext n
        rw [coeff, show adaptedONB T h_inf n = adaptedBasis T h_inf n from
          (adaptedBasis_eq_ONB T h_inf).symm ▸ rfl]
        rw [show @inner ℝ H _ (adaptedBasis T h_inf n) x =
          (adaptedBasis T h_inf).repr x n from
          ((adaptedBasis T h_inf).repr_apply_apply x n).symm]
        simp [sq_abs, ← hc_def]
    _ = ‖c‖^2 := hlp'.symm
    _ = ‖x‖^2 := by rw [hiso]

/-! ## Noise Measure -/

/-- The space of noise sequences ξ : ℕ → ℝ. -/
abbrev NoiseSpace := ℕ → ℝ

/-- The infinite product measure of standard Gaussians. -/
private noncomputable def noiseMeasure : Measure NoiseSpace :=
  Measure.infinitePi (fun _ => gaussianReal 0 1)

private instance noiseMeasure_isProbability : IsProbabilityMeasure noiseMeasure :=
  Measure.instIsProbabilityMeasureForallInfinitePi _

/-! ## Series Limit Map -/

/-- Predicate: ξ is "good" if the K-valued Gaussian series converges. -/
private def IsGoodNoise (T : E →L[ℝ] H) (h_inf : ¬ FiniteDimensional ℝ H)
    (ξ : NoiseSpace) : Prop :=
  Summable (fun n => ξ n • factV T h_inf n)

/-- K-valued summability implies scalar summability for each test function. -/
private lemma IsGoodNoise.scalar_summable
    {T : E →L[ℝ] H} {h_inf : ¬ FiniteDimensional ℝ H} {ξ : NoiseSpace}
    (h : IsGoodNoise T h_inf ξ) (f : E) :
    Summable (fun n => ξ n * coeff T h_inf n f) := by
  have hK := (innerSL ℝ (factJ T h_inf f)).summable h
  exact hK.congr fun n => by
    rw [innerSL_apply_apply, real_inner_smul_right, real_inner_comm,
        ← factV_inner_eq T h_inf f n]; rfl

/-- The series limit as a raw function. -/
private noncomputable def seriesLimitFun
    (T : E →L[ℝ] H) (h_inf : ¬ FiniteDimensional ℝ H)
    (ξ : NoiseSpace) (f : E) : ℝ :=
  if IsGoodNoise T h_inf ξ
  then ∑' n, ξ n * coeff T h_inf n f
  else 0

/-- The series limit as a linear map. -/
private noncomputable def seriesLimitLinear
    (T : E →L[ℝ] H) (h_inf : ¬ FiniteDimensional ℝ H)
    (ξ : NoiseSpace) : E →ₗ[ℝ] ℝ where
  toFun := seriesLimitFun T h_inf ξ
  map_add' := fun f g => by
    simp only [seriesLimitFun]
    split_ifs with h
    · have hf := h.scalar_summable f
      have hg := h.scalar_summable g
      have heq : (fun n => ξ n * coeff T h_inf n (f + g)) =
                 (fun n => ξ n * coeff T h_inf n f + ξ n * coeff T h_inf n g) := by
        ext n; simp [coeff, inner_add_right, mul_add]
      conv_lhs => rw [heq]
      exact hf.tsum_add hg
    · simp
  map_smul' := fun c f => by
    simp only [seriesLimitFun, RingHom.id_apply, smul_eq_mul]
    split_ifs with h
    · have hf := h.scalar_summable f
      have heq : (fun n => ξ n * coeff T h_inf n (c • f)) =
                 (fun n => c • (ξ n * coeff T h_inf n f)) := by
        ext n; simp [coeff, real_inner_smul_right, smul_eq_mul]; ring
      conv_lhs => rw [heq]
      rw [hf.tsum_const_smul c, smul_eq_mul]
    · rw [mul_zero]

/-- The series limit as an element of Configuration E (weak dual).
    Continuity is proved via CLM composition through K. -/
private noncomputable def seriesLimit
    (T : E →L[ℝ] H) (h_inf : ¬ FiniteDimensional ℝ H)
    (ξ : NoiseSpace) : Configuration E :=
  ⟨seriesLimitLinear T h_inf ξ, by
    by_cases h : IsGoodNoise T h_inf ξ
    · -- Good noise: functional = ⟨w, j(·)⟩_K, continuous via CLM composition
      let w : FactK T h_inf := ∑' n, ξ n • factV T h_inf n
      have h_eq : ∀ f, ∑' n, ξ n * coeff T h_inf n f =
          @inner ℝ (FactK T h_inf) _ w (factJ T h_inf f) := by
        intro f
        simp_rw [coeff, factV_inner_eq T h_inf f]
        have h_clm := (innerSL ℝ (factJ T h_inf f)).map_tsum h
        rw [show @inner ℝ (FactK T h_inf) _ w (factJ T h_inf f) =
              @inner ℝ (FactK T h_inf) _ (factJ T h_inf f) w from
                real_inner_comm _ _]
        rw [show @inner ℝ (FactK T h_inf) _ (factJ T h_inf f) w =
              (innerSL ℝ (factJ T h_inf f)) w from rfl]
        rw [h_clm]
        congr 1; ext n
        simp [innerSL_apply_apply, real_inner_comm (factV T h_inf n) (factJ T h_inf f)]
      have h_cont : Continuous (fun f => @inner ℝ (FactK T h_inf) _ w (factJ T h_inf f)) :=
        ((innerSL ℝ w).comp (factJ T h_inf)).continuous
      convert h_cont using 1
      ext f
      show seriesLimitFun T h_inf ξ f = _
      simp only [seriesLimitFun, if_pos h]
      exact h_eq f
    · -- Bad noise: functional is 0, trivially continuous
      have h0 : (seriesLimitLinear T h_inf ξ).toFun = fun _ => 0 := by
        ext f; show seriesLimitFun T h_inf ξ f = 0
        simp only [seriesLimitFun, if_neg h]
      simp only [h0]
      exact continuous_const⟩

set_option maxHeartbeats 1600000 in
/-- **Gaussian series convergence**: If ∑ ‖vₙ‖ < ∞ and ξₙ ~ iid N(0,1),
    then ∑ ξₙ • vₙ converges almost surely in K.

    Proof uses Tonelli's theorem:
    E[∑ ‖ξₙ vₙ‖] = ∑ E[|ξₙ|] ‖vₙ‖ = C · ∑ ‖vₙ‖ < ∞
    so ∑ ‖ξₙ vₙ‖ < ∞ a.s., hence ∑ ξₙ vₙ converges a.s. -/
private theorem hilbert_gaussian_series_converges
    {K : Type*} [NormedAddCommGroup K] [InnerProductSpace ℝ K] [CompleteSpace K]
    (v : ℕ → K) (hv : Summable (fun n => ‖v n‖)) :
    ∀ᵐ ξ ∂(Measure.infinitePi (fun _ : ℕ => gaussianReal 0 1)),
      Summable (fun n => ξ n • v n) := by
  set μ := fun _ : ℕ => gaussianReal (0 : ℝ) 1
  set P := Measure.infinitePi μ
  -- The first absolute moment C = E[|ξ|] of standard Gaussian is finite
  set C := ∫⁻ x : ℝ, ↑‖x‖₊ ∂(μ 0) with hC_def
  have hC : C ≠ ⊤ := by
    have h := (IsGaussian.integrable_dual (μ 0)
      (ContinuousLinearMap.id ℝ ℝ)).hasFiniteIntegral
    rw [hasFiniteIntegral_def] at h
    simp only [ContinuousLinearMap.id_apply, enorm_eq_nnnorm] at h
    exact h.ne
  -- Convert ∑ ‖vₙ‖ < ∞ to NNReal and ENNReal summability
  have hv_nn : Summable (fun n => ‖v n‖₊) := by
    have h1 : (fun n => (‖v n‖₊ : ℝ)) = fun n => ‖v n‖ := funext fun n => coe_nnnorm (v n)
    exact NNReal.summable_coe.mp (h1 ▸ hv)
  have hv_ennreal : tsum (fun n => ((‖v n‖₊ : NNReal) : ENNReal)) ≠ ⊤ :=
    ENNReal.tsum_coe_ne_top_iff_summable.mpr hv_nn
  -- Measurability of each summand ‖ξ n‖₊ * ‖v n‖₊
  have h_meas_f : ∀ n, Measurable (fun ξ : ℕ → ℝ => (‖ξ n‖₊ : ENNReal)) :=
    fun n => (measurable_pi_apply n).nnnorm.coe_nnreal_ennreal
  have h_meas : ∀ n, Measurable (fun ξ : ℕ → ℝ => (‖ξ n‖₊ : ENNReal) * ‖v n‖₊) :=
    fun n => (h_meas_f n).mul_const _
  -- Each marginal integral equals C (by measure-preserving projection)
  have h_eval : ∀ n, ∫⁻ ξ, (‖ξ n‖₊ : ENNReal) ∂P = C := fun n =>
    (measurePreserving_eval_infinitePi μ n).lintegral_comp
      measurable_nnnorm.coe_nnreal_ennreal
  -- Main estimate: E[∑ |ξₙ| ‖vₙ‖] = C · ∑ ‖vₙ‖ < ∞
  have h_finite : ∫⁻ ξ, ∑' n, ((‖ξ n‖₊ : ENNReal) * ‖v n‖₊) ∂P ≠ ⊤ := by
    rw [lintegral_tsum (fun n => (h_meas n).aemeasurable)]
    simp_rw [lintegral_mul_const _ (h_meas_f _), h_eval, ENNReal.tsum_mul_left]
    exact ENNReal.mul_ne_top hC hv_ennreal
  -- A.e. finiteness from finite lintegral
  have h_ae := ae_lt_top (Measurable.tsum h_meas) h_finite
  -- Convert ∑ ‖ξₙ‖₊ * ‖vₙ‖₊ < ⊤ to Summable (fun n => ξ n • v n)
  exact h_ae.mono fun ξ hξ => by
    simp_rw [← ENNReal.coe_mul] at hξ
    have hξ' : Summable (fun n => ‖ξ n‖₊ * ‖v n‖₊) :=
      ENNReal.tsum_coe_ne_top_iff_summable.mp hξ.ne
    have hξ'' : Summable (fun n => ‖ξ n • v n‖) := by
      refine (NNReal.summable_coe.mpr hξ').congr fun n => ?_
      rw [NNReal.coe_mul, coe_nnnorm, coe_nnnorm, norm_smul]
    exact hξ''.of_norm

/-- IsGoodNoise holds a.e. under noiseMeasure. -/
private theorem isGoodNoise_ae (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) :
    ∀ᵐ ξ ∂noiseMeasure, IsGoodNoise T h_inf ξ :=
  hilbert_gaussian_series_converges (factV T h_inf) (factV_summable T h_inf)

/-- The limUnder of partial sums is measurable (for a fixed test function). -/
private theorem limUnder_partialSum_measurable
    (T : E →L[ℝ] H) (h_inf : ¬ FiniteDimensional ℝ H)
    (φ : E) :
    Measurable (fun ξ : NoiseSpace =>
      Filter.limUnder Filter.atTop (fun N => ∑ k ∈ Finset.range N,
        ξ k * coeff T h_inf k φ)) := by
  apply MeasureTheory.StronglyMeasurable.measurable
  apply MeasureTheory.StronglyMeasurable.limUnder
  intro N
  have : (fun ξ : NoiseSpace => ∑ k ∈ Finset.range N, ξ k * coeff T h_inf k φ) =
    ∑ k ∈ Finset.range N, (fun ξ => ξ k * coeff T h_inf k φ) := by
    ext ξ; simp [Finset.sum_apply]
  rw [this]
  apply Finset.stronglyMeasurable_sum
  intro k _
  exact ((measurable_pi_apply k).stronglyMeasurable).mul stronglyMeasurable_const

/-- The series limit is AE-measurable. -/
private theorem aemeasurable_seriesLimit (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) :
    AEMeasurable (seriesLimit T h_inf) noiseMeasure := by
  have h_null : noiseMeasure {ξ | ¬IsGoodNoise T h_inf ξ} = 0 :=
    ae_iff.mp (isGoodNoise_ae T h_inf)
  obtain ⟨N, hN_sub, hN_meas, hN_null⟩ := exists_measurable_superset_of_null h_null
  set S := Nᶜ with hS_def
  have hS_meas : MeasurableSet S := hN_meas.compl
  have hS_sub : ∀ ξ ∈ S, IsGoodNoise T h_inf ξ := fun ξ hξ => by
    by_contra h; exact hξ (hN_sub h)
  have hS_null : noiseMeasure Sᶜ = 0 := by rwa [hS_def, compl_compl]
  let g : NoiseSpace → Configuration E := S.indicator (seriesLimit T h_inf)
  have hg_meas : @Measurable NoiseSpace (Configuration E) _ instMeasurableSpaceConfiguration g := by
    apply configuration_measurable_of_eval_measurable
    intro φ
    have h_eval : (fun ξ => g ξ φ) = S.indicator (fun ξ => seriesLimit T h_inf ξ φ) := by
      ext ξ; simp only [g, Set.indicator]; split_ifs <;> rfl
    rw [h_eval]
    have h_eq : S.indicator (fun ξ => seriesLimit T h_inf ξ φ) =
        S.indicator (fun ξ => Filter.limUnder Filter.atTop
          (fun N => ∑ k ∈ Finset.range N, ξ k * coeff T h_inf k φ)) := by
      ext ξ; simp only [Set.indicator]; split_ifs with h
      · have hgood := hS_sub ξ h
        show seriesLimitFun T h_inf ξ φ = _
        simp only [seriesLimitFun, if_pos hgood]
        exact ((hgood.scalar_summable φ).hasSum.tendsto_sum_nat).limUnder_eq.symm
      · rfl
    rw [h_eq]
    exact (limUnder_partialSum_measurable T h_inf φ).indicator hS_meas
  have hg_ae : seriesLimit T h_inf =ᵐ[noiseMeasure] g := by
    rw [Filter.EventuallyEq, Filter.eventually_iff_exists_mem]
    exact ⟨S, mem_ae_iff.mpr hS_null, fun ξ hξ => by simp [g, hξ]⟩
  exact ⟨g, hg_meas, hg_ae⟩

/-! ## Isometric Embedding for Finite-Dimensional H

When H is finite-dimensional, we embed it isometrically into ℓ² and apply
the infinite-dimensional construction there. This allows the public API
to work without requiring `¬ FiniteDimensional ℝ H`. -/

/-- Isometric embedding of a finite-dimensional H into ℓ² via an ONB.
    Maps x ↦ ∑_i ⟨b_i, x⟩ • e_i where b is the stdOrthonormalBasis of H. -/
private noncomputable def hilbertEmbedding (hfin : FiniteDimensional ℝ H) : H →L[ℝ] ell2' := by
  haveI := hfin
  exact ∑ i : Fin (Module.finrank ℝ H),
    (innerSL ℝ ((stdOrthonormalBasis ℝ H) i)).smulRight (ell2_basis i.val)

omit [CompleteSpace H] [SeparableSpace H] in
/-- The embedding maps x to ∑_i ⟨b_i, x⟩ • e_i. -/
private lemma hilbertEmbedding_apply (hfin : FiniteDimensional ℝ H) (x : H) :
    hilbertEmbedding hfin x = ∑ i : Fin (Module.finrank ℝ H),
      @inner ℝ H _ ((stdOrthonormalBasis ℝ H) i) x • ell2_basis i.val := by
  simp only [hilbertEmbedding, sum_apply,
             ContinuousLinearMap.smulRight_apply, innerSL_apply_apply]

omit [CompleteSpace H] [SeparableSpace H] in
/-- The embedding preserves inner products.
    Both sides equal Σ_i ⟨b_i, x⟩ * ⟨b_i, y⟩: the RHS by Parseval for the
    finite ONB, the LHS by orthonormality of ell2_basis. -/
private lemma hilbertEmbedding_inner (hfin : FiniteDimensional ℝ H) (x y : H) :
    @inner ℝ ell2' _ (hilbertEmbedding hfin x) (hilbertEmbedding hfin y) =
    @inner ℝ H _ x y := by
  haveI := hfin
  rw [hilbertEmbedding_apply, hilbertEmbedding_apply]
  -- Expand using bilinearity of inner product
  simp only [sum_inner, inner_sum, inner_smul_left, inner_smul_right,
             starRingEnd_apply, star_trivial]
  -- Use orthonormality of ell2_basis to collapse the double sum
  simp_rw [show ∀ (i j : Fin (Module.finrank ℝ H)),
      @inner ℝ ell2' _ (ell2_basis i.val) (ell2_basis j.val) =
      if i = j then 1 else 0 from fun i j => by
    rcases eq_or_ne i j with h | h
    · rw [if_pos h, h, real_inner_self_eq_norm_sq,
          ell2_basis_orthonormal.1 j.val, one_pow]
    · rw [if_neg h]; exact ell2_basis_orthonormal.2 (show i.val ≠ j.val from
        fun h' => h (Fin.ext h'))]
  simp only [mul_ite, mul_one, mul_zero, Finset.sum_ite_eq', Finset.mem_univ, if_true]
  -- RHS: Parseval for the finite ONB
  rw [← (stdOrthonormalBasis ℝ H).repr.inner_map_map x y, PiLp.inner_apply]
  simp only [RCLike.inner_apply, conj_trivial,
             OrthonormalBasis.repr_apply_apply, mul_comm]

/-! ## Measure Definition -/

/-- The Gaussian measure on E' constructed from a CLM T : E → H.

    This is the pushforward of the infinite product of iid N(0,1) random
    variables under the series limit map. Works for both finite- and
    infinite-dimensional H: in the finite-dimensional case, H is embedded
    isometrically into ℓ² and the construction is applied there. -/
def measure (T : E →L[ℝ] H) :
    @Measure (Configuration E) instMeasurableSpaceConfiguration :=
  if hfin : FiniteDimensional ℝ H then
    @Measure.map NoiseSpace (Configuration E) _ instMeasurableSpaceConfiguration
      (seriesLimit ((hilbertEmbedding hfin).comp T) ell2_not_finiteDimensional) noiseMeasure
  else
    @Measure.map NoiseSpace (Configuration E) _ instMeasurableSpaceConfiguration
      (seriesLimit T hfin) noiseMeasure

instance measure_isProbability (T : E →L[ℝ] H) :
    @IsProbabilityMeasure (Configuration E) instMeasurableSpaceConfiguration
      (measure T) := by
  unfold measure; split_ifs with hfin
  · exact Measure.isProbabilityMeasure_map
      (aemeasurable_seriesLimit ((hilbertEmbedding hfin).comp T) ell2_not_finiteDimensional)
  · exact Measure.isProbabilityMeasure_map (aemeasurable_seriesLimit T hfin)

/-! ## Characteristic Functional -/

/-- For a.e. ξ, partial sums converge to the series limit. -/
private theorem partialSum_tendsto (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) (f : E) :
    ∀ᵐ ξ ∂noiseMeasure, Filter.Tendsto
      (fun N => ∑ k ∈ Finset.range N, ξ k * coeff T h_inf k f)
      Filter.atTop
      (nhds (seriesLimit T h_inf ξ f)) := by
  filter_upwards [isGoodNoise_ae T h_inf] with ξ hgood
  have hsummable := hgood.scalar_summable f
  have hlimit : (seriesLimit T h_inf ξ f) =
      ∑' n, ξ n * coeff T h_inf n f := by
    show seriesLimitFun T h_inf ξ f = _
    simp only [seriesLimitFun, if_pos hgood]
  rw [hlimit]
  exact hsummable.hasSum.tendsto_sum_nat

/-- Coordinate projections are iIndepFun under the noise measure. -/
private theorem coord_projections_indep :
    iIndepFun (m := fun _ => inferInstance)
    (fun i (ξ : NoiseSpace) ↦ ξ i) noiseMeasure := by
  have h := @iIndepFun_infinitePi ℕ (fun _ => ℝ) (fun _ => inferInstance)
    (fun _ => ℝ) (fun _ => inferInstance) (fun _ => gaussianReal 0 1) _
    (fun _ => id) (fun _ => measurable_id)
  simp only [id_eq] at h
  exact h

/-- Independence factorization for integrals over the noise measure. -/
private theorem noise_integral_factor (n : ℕ)
    (g : NoiseSpace → ℂ) (h : ℝ → ℂ)
    (hg : ∀ ξ ξ', (∀ k < n, ξ k = ξ' k) → g ξ = g ξ')
    (hg_meas : Measurable g)
    (hh_meas : Measurable h) :
    ∫ ξ, g ξ * h (ξ n) ∂noiseMeasure =
      (∫ ξ, g ξ ∂noiseMeasure) * (∫ x, h x ∂gaussianReal 0 1) := by
  have hdisj : Disjoint (Finset.range n) ({n} : Finset ℕ) := by
    rw [Finset.disjoint_left]; intro a ha hn
    simp only [Finset.mem_singleton] at hn
    simp only [Finset.mem_range] at ha; omega
  let X : NoiseSpace → (Finset.range n → ℝ) := fun ξ i => ξ i
  let Y : NoiseSpace → ℝ := fun ξ => ξ n
  have hX_meas : Measurable X := measurable_pi_lambda _ (fun i => measurable_pi_apply _)
  have hY_meas : Measurable Y := measurable_pi_apply n
  let Y' : NoiseSpace → ({n} : Finset ℕ) → ℝ := fun ξ i => ξ i
  have hindep' : IndepFun X Y' noiseMeasure :=
    coord_projections_indep.indepFun_finset (Finset.range n) {n} hdisj
      (fun _ => measurable_pi_apply _)
  have hindep : IndepFun X Y noiseMeasure := by
    have hY_eq : Y = (fun f => f ⟨n, Finset.mem_singleton_self n⟩) ∘ Y' := rfl
    rw [hY_eq]
    exact hindep'.comp measurable_id (measurable_pi_apply _)
  let extend_to_full : (Finset.range n → ℝ) → NoiseSpace :=
    fun f k => if hk : k < n then f ⟨k, Finset.mem_range.mpr hk⟩ else 0
  have hext_meas : Measurable extend_to_full := by
    apply measurable_pi_lambda; intro k; by_cases hk : k < n
    · simp only [extend_to_full, dif_pos hk]; exact measurable_pi_apply _
    · simp only [extend_to_full, dif_neg hk]; exact measurable_const
  let g' : (Finset.range n → ℝ) → ℂ := g ∘ extend_to_full
  have hg'_meas : Measurable g' := hg_meas.comp hext_meas
  have hg_factor : ∀ ξ, g ξ = g' (X ξ) := by
    intro ξ; simp only [g', Function.comp_apply]
    apply hg; intro k hk; simp only [extend_to_full, X, dif_pos hk]
  have heval : MeasurePreserving Y noiseMeasure (gaussianReal 0 1) :=
    measurePreserving_eval_infinitePi (fun _ => gaussianReal 0 1) n
  have hY_int : ∫ x, h x ∂gaussianReal 0 1 = ∫ ξ, h (Y ξ) ∂noiseMeasure := by
    simp only [← heval.map_eq]
    exact integral_map hY_meas.aemeasurable hh_meas.aestronglyMeasurable
  rw [hY_int]
  conv_lhs => arg 2; ext ξ; rw [hg_factor ξ]
  conv_rhs => arg 1; arg 2; ext ξ; rw [hg_factor ξ]
  conv_lhs => arg 2; ext ξ; rw [show h (ξ n) = h (Y ξ) from rfl]
  exact hindep.integral_fun_comp_mul_comp hX_meas.aemeasurable hY_meas.aemeasurable
    hg'_meas.aestronglyMeasurable hh_meas.aestronglyMeasurable

/-- Characteristic functional of the partial sum. -/
private theorem partial_sum_charFun (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) (f : E) (N : ℕ) :
    ∫ ξ, Complex.exp (Complex.I * ↑(∑ k ∈ Finset.range N,
      ξ k * coeff T h_inf k f)) ∂noiseMeasure =
    Complex.exp (-(1/2 : ℂ) * ∑ k ∈ Finset.range N,
      (coeff T h_inf k f : ℂ) ^ 2) := by
  simp_rw [Complex.ofReal_sum, Complex.ofReal_mul]
  conv_lhs => arg 2; ext ξ; rw [Finset.mul_sum]
  simp_rw [Complex.exp_sum]
  induction N with
  | zero =>
    simp only [Finset.range_zero, Finset.prod_empty, Finset.sum_empty, mul_zero,
      Complex.exp_zero, integral_const]
    have h : noiseMeasure.real Set.univ = 1 := by
      simp only [Measure.real, measure_univ, ENNReal.toReal_one]
    simp only [h, one_smul]
  | succ n ih =>
    simp only [Finset.range_add_one, Finset.prod_insert Finset.notMem_range_self,
               Finset.sum_insert Finset.notMem_range_self]
    conv_lhs => arg 2; ext xi; rw [mul_comm]
    let g := fun xi : NoiseSpace =>
      ∏ x ∈ Finset.range n, Complex.exp (Complex.I * (↑(xi x) * ↑(coeff T h_inf x f)))
    let h := fun x : ℝ => Complex.exp (Complex.I * (↑x * ↑(coeff T h_inf n f)))
    have hfactor : ∫ xi, g xi * h (xi n) ∂noiseMeasure =
        (∫ xi, g xi ∂noiseMeasure) * (∫ x, h x ∂gaussianReal 0 1) := by
      apply noise_integral_factor
      · intro xi xi' hcoords; simp only [g]
        apply Finset.prod_congr rfl
        intro k hk; rw [hcoords k (Finset.mem_range.mp hk)]
      · exact (by continuity : Continuous g).measurable
      · exact (by continuity : Continuous h).measurable
    rw [hfactor, ih]
    have hcf : ∫ x, h x ∂gaussianReal 0 1 =
        Complex.exp (-(1/2 : ℂ) * (coeff T h_inf n f : ℂ) ^ 2) := by
      simp only [h]
      have hcf' := charFun_gaussianReal (μ := 0) (v := 1) (coeff T h_inf n f)
      simp only [charFun_apply_real] at hcf'
      convert hcf' using 1
      · apply integral_congr_ae; filter_upwards with x; congr 1; ring
      · congr 1; simp only [Complex.ofReal_zero, mul_zero, zero_mul, zero_sub]
        norm_cast; push_cast; ring
    rw [hcf, ← Complex.exp_add]; congr 1; ring

/-- Core computation: the noise integral equals the Gaussian characteristic functional.
    Proved for infinite-dimensional H where the nuclear SVD applies. -/
private theorem charFun_noise_integral (T : E →L[ℝ] H)
    (h_inf : ¬ FiniteDimensional ℝ H) (f : E) :
    ∫ ξ : NoiseSpace, Complex.exp (Complex.I * ↑(seriesLimit T h_inf ξ f)) ∂noiseMeasure =
    Complex.exp (-(1/2 : ℂ) * ↑(@inner ℝ H _ (T f) (T f))) := by
  -- Parseval identity for the target
  have hparseval : (@inner ℝ H _ (T f) (T f) : ℂ) =
      ∑' n, (coeff T h_inf n f : ℂ) ^ 2 := by
    have h := coeff_parseval T h_inf f
    rw [← real_inner_self_eq_norm_sq] at h
    simp only [← Complex.ofReal_pow] at h ⊢
    rw [← h]; simp only [Complex.ofReal_tsum]
  rw [hparseval]
  -- Summability of the complex coefficient series
  have hsummable : Summable (fun n => (coeff T h_inf n f : ℂ) ^ 2) := by
    have hr := coeff_sq_summable T h_inf f
    simp only [← Complex.ofReal_pow]
    exact Complex.summable_ofReal.mpr hr
  -- The target limit
  have hcf_limit : Filter.Tendsto
      (fun N => Complex.exp (-(1/2 : ℂ) * ∑ k ∈ Finset.range N,
        (coeff T h_inf k f : ℂ) ^ 2))
      Filter.atTop
      (nhds (Complex.exp (-(1/2 : ℂ) * ∑' n, (coeff T h_inf n f : ℂ) ^ 2))) := by
    apply Complex.continuous_exp.continuousAt.tendsto.comp
    apply Filter.Tendsto.const_mul
    exact hsummable.hasSum.tendsto_sum_nat
  -- By partial_sum_charFun, for each N:
  have hpartial : ∀ N, ∫ ξ, Complex.exp (Complex.I *
      ↑(∑ k ∈ Finset.range N, ξ k * coeff T h_inf k f)) ∂noiseMeasure =
      Complex.exp (-(1/2 : ℂ) * ∑ k ∈ Finset.range N,
        (coeff T h_inf k f : ℂ) ^ 2) :=
    fun N => partial_sum_charFun T h_inf f N
  -- DCT argument
  have hdct : Filter.Tendsto
      (fun N => ∫ ξ, Complex.exp (Complex.I *
        ↑(∑ k ∈ Finset.range N, ξ k * coeff T h_inf k f)) ∂noiseMeasure)
      Filter.atTop
      (nhds (∫ ξ, Complex.exp (Complex.I * ↑(seriesLimit T h_inf ξ f)) ∂noiseMeasure)) := by
    apply MeasureTheory.tendsto_integral_of_dominated_convergence (fun _ => (1 : ℝ))
    · intro N; apply Continuous.aestronglyMeasurable; continuity
    · exact integrable_const 1
    · intro N; filter_upwards with ξ
      rw [mul_comm, Complex.norm_exp_ofReal_mul_I]
    · filter_upwards [partialSum_tendsto T h_inf f] with ξ htends
      apply Complex.continuous_exp.continuousAt.tendsto.comp
      apply Filter.Tendsto.const_mul
      exact Complex.continuous_ofReal.continuousAt.tendsto.comp htends
  -- Combine via uniqueness of limits
  exact tendsto_nhds_unique hdct
    (Filter.Tendsto.congr (fun N => (hpartial N).symm) hcf_limit)

/-- **Characteristic functional identity**: The constructed measure has
    characteristic functional exp(-½ ‖T(f)‖²).

    E[exp(i⟨ω, f⟩)] = exp(-½ ‖T(f)‖²_H) = exp(-½ C(f,f))

    Proof: independence of the ξₙ gives a product of 1D Gaussian CFs,
    which telescopes via Parseval's identity. For finite-dimensional H,
    we embed H isometrically into ℓ² and apply the infinite-dimensional
    construction there. -/
theorem charFun (T : E →L[ℝ] H) (f : E) :
    ∫ ω : Configuration E,
      Complex.exp (Complex.I * ↑(ω f)) ∂(measure T) =
    Complex.exp (-(1/2 : ℂ) * ↑(@inner ℝ H _ (T f) (T f))) := by
  have hmeas_integrand : ∀ (μ : @Measure (Configuration E) instMeasurableSpaceConfiguration),
      AEStronglyMeasurable (fun ω : Configuration E =>
        Complex.exp (Complex.I * ↑(ω f))) μ :=
    fun _ => (Complex.continuous_exp.measurable.comp
      (measurable_const.mul (Complex.continuous_ofReal.measurable.comp
        (configuration_eval_measurable f)))).aestronglyMeasurable
  unfold measure
  by_cases hfin : FiniteDimensional ℝ H
  · -- Finite-dimensional case: embed H into ℓ² isometrically
    rw [dif_pos hfin]
    rw [integral_map
      (aemeasurable_seriesLimit ((hilbertEmbedding hfin).comp T) ell2_not_finiteDimensional)
      (hmeas_integrand _)]
    -- Apply core computation for (embed ∘ T) on ℓ²
    have h := charFun_noise_integral ((hilbertEmbedding hfin).comp T)
      ell2_not_finiteDimensional f
    -- Rewrite ⟨embed(Tf), embed(Tf)⟩_ℓ² = ⟨Tf, Tf⟩_H
    simp only [ContinuousLinearMap.comp_apply] at h
    rwa [hilbertEmbedding_inner] at h
  · -- Infinite-dimensional case: direct argument
    rw [dif_neg hfin]
    rw [integral_map (aemeasurable_seriesLimit T hfin) (hmeas_integrand _)]
    exact charFun_noise_integral T hfin f

end GaussianField