AutoSamplingTheory.ExampleCases.SmoothedPicardHMC.EnhancedTerminalExecution

No declaration docstring.

private theorem compatible {K d q Δ L : ℝ} {J : ℕ}
    (hK : 1 ≤ K) (hd : 1 ≤ d) (hq : 2 ≤ q)
    (hΔ : 0 < Δ) (hΔ1 : Δ ≤ 1/2)
    (hL : L = q+Real.log (K*d*q/Δ)) (hL2 : 2 ≤ L)
    (hJ : (J:ℝ) ≤ (24+4*Real.log 1024)*L) :
    let eps := Δ/Real.sqrt ((J:ℝ)+1)
    let H := q+Real.log (1/eps)
    0 < eps ∧ eps ≤ 1/2 ∧ H ≤ 3*L ∧
      64*(Real.sqrt (d*H)+H) ≤ 1024*(Real.sqrt (d*L)+L) := by
  intro eps H
  have hL0 : 0 < L := by linarith
  have hj0 : 0 < (J:ℝ)+1 := by positivity
  have hs0 : 0 < Real.sqrt ((J:ℝ)+1) := Real.sqrt_pos.2 hj0
  have hs1 : 1 ≤ Real.sqrt ((J:ℝ)+1) := by
    have := Real.sqrt_le_sqrt (show (1:ℝ) ≤ (J:ℝ)+1 by exact le_add_of_nonneg_left (Nat.cast_nonneg J))
    simpa only [Real.sqrt_one] using this
  have heps0 : 0 < eps := div_pos hΔ hs0
  have heps1 : eps ≤ 1/2 := by
    apply le_trans (div_le_self hΔ.le hs1) hΔ1
  have hlog2 : Real.log 2 ≤ 1 := by
    have := Real.log_le_sub_one_of_pos (show (0:ℝ)<2 by norm_num)
    linarith
  have hlog1024 : Real.log 1024 ≤ 10 := by
    have hp := Real.log_pow (2:ℝ) 10
    norm_num at hp
    linarith
  have hj64 : (J:ℝ) ≤ 64*L := by
    nlinarith
  have hj65 : (J:ℝ)+1 ≤ 65*L := by linarith
  have hlog65 : Real.log 65 ≤ 7 := by
    have hm := Real.log_le_log (show (0:ℝ)<65 by norm_num) (show (65:ℝ) ≤ 2^7 by norm_num)
    rw [Real.log_pow] at hm
    norm_num at hm
    linarith
  have hjlog : Real.log ((J:ℝ)+1) ≤ 4*L := by
    have hm := Real.log_le_log hj0 hj65
    rw [Real.log_mul (by norm_num : (65:ℝ)≠0) hL0.ne'] at hm
    have hl := Real.log_le_sub_one_of_pos hL0
    linarith
  have hprod : 1 ≤ K*d*q := by nlinarith [mul_le_mul hK hd (by norm_num : (0:ℝ)≤1) (by linarith : 0≤K)]
  have hbase : q+Real.log (1/Δ) ≤ L := by
    rw [hL]
    apply add_le_add le_rfl
    apply Real.log_le_log (one_div_pos.2 hΔ)
-- Source excerpt truncated; follow the exact source link.

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/EnhancedTerminalExecution.lean:46published source at 0e31a3cda412

Excerpt truncated; the exact source link is authoritative.

No declaration docstring.

private theorem termination
    (b η : S → ℝ) (hb : Measurable b) (F : S × Y → S) (hFm : Measurable F)
    (κ c q Δ γ C : ℝ) (d : ℕ) (s0 : S)
    (hκ : 1 ≤ κ) (hc : 0 < c) (hc1 : c < 1/4) (hb0 : 0 ≤ b s0)
    (hη : ∀ s, 0 < η s ∧ η s ≤ c) (hd : 0 < d) (hq : 2 ≤ q)
    (hΔ : 0 < Δ) (hΔ1 : Δ ≤ 1/2) (hγ : 0 < γ) (hγ1 : γ ≤ 1) (hC : 8 ≤ C)
    (Q : Kernel S Y) [IsMarkovKernel Q] (P : Kernel S S) [IsMarkovKernel P]
    (Lt : Kernel S Z) [IsMarkovKernel Lt] :
    let K := fun r : ℝ => (1+r)/(κ⁻¹+r)
    let τ := fun r : ℝ => if 2 ≤ K r then K r else c
    let L := q+Real.log (K (b s0)*d*q/Δ)
    let B := γ/(Real.sqrt (d*L)+L)
    let J := Nat.ceil (C*Real.log (Real.exp 1*K (b s0)/B))
    (∀ s y, b (F (s,y)) = b s+((η s+τ (b s))/(1+b s))⁻¹) →
    (∀ s, P s=if B⁻¹ ≤ b s then Measure.dirac s else (Q s).map (fun y => F (s,y))) →
    (P^J) s0 {s | B⁻¹ ≤ b s}=1 ∧
    (∀ m, (Lt ∘ₖ (P^(J+m))) s0=(Lt ∘ₖ (P^J)) s0) := by
  classical
  let statePath (G : S × Y → S) (s : S) (ys : ℕ → Y) : ℕ → S :=
    Nat.rec s (fun n x => G (x,ys n))
  have statePath_shift (G : (S) × Y → (S)) (n : ℕ) (s : (S)) (ys : ℕ → Y) :
      statePath G s ys (n+1) = statePath G (G (s,ys 0)) (fun k => ys (k+1)) n := by
    induction n with
    | zero => rfl
    | succ n ih =>
      change G (statePath G s ys (n+1),ys (n+1)) =
        G (statePath G (G (s,ys 0)) (fun k => ys (k+1)) n,ys (n+1))
      rw [ih]
  have path_endpoint_mass (G : (S) × Y → (S)) (hG : Measurable G)
      (Q : Kernel (S) Y) [IsMarkovKernel Q] (P : Kernel (S) (S)) [IsMarkovKernel P]
      (hP : ∀ s, P s = (Q s).map (fun y => G (s,y)))
      (D : Set (S)) (hD : MeasurableSet D) (n : ℕ) (s : (S))
      (hpaths : ∀ ys : ℕ → Y, statePath G s ys n ∈ D) [Nonempty Y] :
      (P^n) s D = 1 := by
    have hp (k : ℕ) : IsMarkovKernel (P^k) := by
      induction k with
      | zero => change IsMarkovKernel Kernel.id; infer_instance
      | succ k ih =>
        let := ih
        rw [pow_succ]
        change IsMarkovKernel ((P^k) ∘ₖ P)
        infer_instance
    induction n generalizing s with
    | zero =>
-- Source excerpt truncated; follow the exact source link.

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/EnhancedTerminalExecution.lean:111published source at 0e31a3cda412

Excerpt truncated; the exact source link is authoritative.

No declaration docstring.

private def firstIndex {X : Type*} (q : ℕ → X → ℝ) (b : X → ℝ) (s : X) : ℕ := by
  classical
  exact if h : ∃ n, q n s ≤ b s then Nat.find h else 0

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/EnhancedTerminalExecution.lean:299published source at 0e31a3cda412

No declaration docstring.

private theorem execution {V : E → ℝ} (κ : ℝ≥0) (hκ : 1 ≤ κ)
    (hV : ContDiff ℝ 2 V)
    (hH : ∀ x w : E, (κ:ℝ)⁻¹*‖w‖^2 ≤ fderiv ℝ (fderiv ℝ V) x w w ∧
      fderiv ℝ (fderiv ℝ V) x w w ≤ ‖w‖^2)
    (hd : 0 < Module.finrank ℝ E)
    (η : S → ℝ) (hηm : Measurable η) (c q Δ : ℝ)
    (hc : 0 < c) (hc1 : c < 1/4) (hη : ∀ s, 0 < η s ∧ η s ≤ c)
    (hq : 2 ≤ q) (hΔ : 0 < Δ) (hΔ1 : Δ ≤ 1/2)
    (s0 : S) (M Lt : Kernel S E) [IsMarkovKernel M] [IsMarkovKernel Lt] :
    let d : ℝ := Module.finrank ℝ E
    let K := fun r : ℝ => (1+r)/((κ:ℝ)⁻¹+r)
    let τ := fun s : S => if 2 ≤ K s.1 then K s.1 else c
    let L := q+Real.log (K s0.1*d*q/Δ)
    let B := (1/1024)/(Real.sqrt (d*L)+L)
    let J := Nat.ceil (8*Real.log (Real.exp 1*K s0.1/B))
    let eps := Δ/Real.sqrt ((J:ℝ)+1)
    let v := fun s : S => (η s+τ s)/(1+s.1)
    let bp := fun s : S => s.1+Real.toNNReal (v s)⁻¹
    let obs := fun p : S × (E × E) => p.2.1+Real.sqrt (τ p.1/(1+p.1.1)) • p.2.2
    let center := fun p : S × (E × E) => (bp p.1:ℝ)⁻¹ •
      ((p.1.1:ℝ) • p.1.2.1+(v p.1)⁻¹ • obs p)
    let F := fun p x => V x+(bp p.1:ℝ)/2*‖x-center p‖^2
    let T := fun p x => x-(1+(bp p.1:ℝ))⁻¹ • gradient (F p) x
    let Qn := fun n p => ‖gradient (F p) ((T p)^[n] p.2.1)‖^2
    let N := firstIndex Qn (fun p => ((κ:ℝ)⁻¹+bp p.1)*d)
    let update := fun p : S × (E × E) =>
      (bp p.1,center p,(T p)^[N p] p.2.1,p.1.2.2.2.1+1,
        fun n => Nat.casesOn n (obs p) p.1.2.2.2.2)
    0 < B ∧ 0 < eps ∧ eps ≤ 1/2 ∧
    64*(Real.sqrt (d*(q+Real.log (1/eps)))+(q+Real.log (1/eps))) ≤ B⁻¹ ∧
    ∃ P : Kernel S S, IsMarkovKernel P ∧
      (∀ s, P s=if B⁻¹ ≤ (s.1:ℝ) then Measure.dirac s else
        ((M s).prod (stdGaussian E)).map (fun z => update (s,z))) ∧
      (P^J) s0 {s | B⁻¹ ≤ (s.1:ℝ)}=1 ∧
      ∀ m, (Lt ∘ₖ (P^(J+m))) s0=(Lt ∘ₖ (P^J)) s0 := by
  classical
  intro d K τ L B J eps v bp obs center F T Qn N update
  have hκr : (1:ℝ) ≤ κ := by exact_mod_cast hκ
  have hκ0 : (0:ℝ) < κ := by linarith
  have hκn0 : (0:ℝ≥0) < κ := by exact_mod_cast hκ0
  have hα : (0:ℝ≥0) < κ⁻¹ := inv_pos.2 hκn0
  have hαβ : κ⁻¹ ≤ (1:ℝ≥0) := (inv_le_one₀ hκn0).2 hκ
  have hHm : ∀ x w : E, ((κ⁻¹:ℝ≥0):ℝ)*‖w‖^2 ≤ fderiv ℝ (fderiv ℝ V) x w w ∧
      fderiv ℝ (fderiv ℝ V) x w w ≤ ((1:ℝ≥0):ℝ)*‖w‖^2 := by simpa using hH
-- Source excerpt truncated; follow the exact source link.

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/EnhancedTerminalExecution.lean:303published source at 0e31a3cda412

Excerpt truncated; the exact source link is authoritative.

No declaration docstring.

private def firstIndex (q : ℕ → X → ℝ) (b : X → ℝ) (s : X) : ℕ := by
  classical
  exact if h : ∃ n, q n s ≤ b s then Nat.find h else 0

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/EnhancedTerminalExecution.lean:405published source at 0e31a3cda412

No declaration docstring.

private def scaledCached (V : E → ℝ) (A : S → ℝ) (h g : S → E)
    (p : (S × E) × (ℝ × E)) : ℝ :=
  let s := p.1.1
  let x := p.1.2
  let t := p.2.1
  let z := Real.sqrt (A s) • p.2.2
  inner ℝ ((Real.pi/2) • (Real.cos (Real.pi/2*t) • (x-h s)-
    Real.sin (Real.pi/2*t) • z))
    (g s-gradient V (h s+Real.sin (Real.pi/2*t) • (x-h s)+
      Real.cos (Real.pi/2*t) • z))

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/EnhancedTerminalExecution.lean:409published source at 0e31a3cda412

No declaration docstring.

private abbrev Attempt (X A : Type*) := X × (ℕ × ((ℕ → A) × ℝ))

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/EnhancedTerminalExecution.lean:420published source at 0e31a3cda412

No declaration docstring.

private def attemptLaw {X A : Type*} [MeasurableSpace X] [MeasurableSpace A]
    (q : Measure X) (nu : Measure A) (B : ℝ) (hB : 0 < B) : Measure (Attempt X A) :=
  q.prod ((poissonMeasure (⟨2*B, by positivity⟩ : ℝ≥0)).prod
    ((Measure.infinitePi (fun _ : ℕ => nu)).prod (volume.restrict (Set.Icc 0 1))))

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/EnhancedTerminalExecution.lean:422published source at 0e31a3cda412

No declaration docstring.

private def accepted {X A : Type*} (W : X × A → ℝ) (B : ℝ) : Set (Attempt X A) :=
  {p | p.2.2.2 ≤ ∏ i : Fin p.2.1, (B+W (p.1,p.2.2.1 i.val))/(2*B)}

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/EnhancedTerminalExecution.lean:427published source at 0e31a3cda412

No declaration docstring.

private def output {X A : Type*} (W : X × A → ℝ) (B : ℝ) (x0 : X)
    (omega : ℕ → Attempt X A) : X := by
  classical
  exact if h : ∃ n, omega n ∈ accepted W B then (omega (Nat.find h)).1 else x0

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/EnhancedTerminalExecution.lean:430published source at 0e31a3cda412

No declaration docstring.

private def queryCount {X A : Type*} (W : X × A → ℝ) (B : ℝ)
    (omega : ℕ → Attempt X A) : ℝ≥0∞ := by
  classical
  exact ∑' n : ℕ, if ∀ i < n, omega i ∉ accepted W B then (omega n).2.1 else 0

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/EnhancedTerminalExecution.lean:435published source at 0e31a3cda412

No declaration docstring.

private theorem terminal_extension {V : E → ℝ} {α β : ℝ≥0}
    (hV : ContDiff ℝ 2 V)
    (hH : ∀ x w : E, (α:ℝ)*‖w‖^2 ≤ fderiv ℝ (fderiv ℝ V) x w w ∧
      fderiv ℝ (fderiv ℝ V) x w w ≤ (β:ℝ)*‖w‖^2)
    (hαβ : α ≤ β) (hβ : 0 < (β:ℝ)) (hd : 0 < (Module.finrank ℝ E:ℝ))
    (b : S → ℝ) (u initial : S → E) (hb : Measurable b) (hu : Measurable u) (hi : Measurable initial)
    (Θ ell eps : ℝ) (hΘ : 0 < Θ) (hell : 2 ≤ ell) (heps : 0 < eps) (heps1 : eps ≤ 1/2)
    (hstep : 64*(β:ℝ)*(Real.sqrt ((Module.finrank ℝ E:ℝ)*(ell+Real.log (1/eps)))+
      (ell+Real.log (1/eps))) ≤ Θ) :
    let bb := fun s => max (b s) Θ
    let d : ℝ := Module.finrank ℝ E
    let F := fun s x => V x+bb s/2*‖x-u s‖^2
    let T := fun s x => x-((β:ℝ)+bb s)⁻¹ • gradient (F s) x
    let Qn := fun n s => ‖gradient (F s) ((T s)^[n] (initial s))‖^2
    let N := firstIndex Qn (fun s => d*bb s)
    let xp := fun s => (T s)^[N s] (initial s)
    let A := fun s => (bb s)⁻¹
    let g := fun s => gradient V (xp s)
    let h := fun s => u s-A s • g s
    let ν := (volume.restrict (Ioc (0 : ℝ) 1)).prod (stdGaussian E)
    let W := fun s (p : E × (ℝ × E)) => min 1 (max (-1) (scaledCached V A h g ((s,p.1),p.2)))
    let q := fun s => (stdGaussian E).map (fun z => h s+Real.sqrt (A s) • z)
    let ρ := fun s => Measure.infinitePi (fun _ : ℕ => attemptLaw (q s) ν 1 (by norm_num))
    let π := fun s => ((stdGaussian E).map (fun z => u s+Real.sqrt (A s) • z)).tilted (fun x => -V x)
    (∀ s, Θ ≤ b s → bb s=b s) ∧ Measurable N ∧ Measurable xp ∧
    (∀ s, Qn (N s) s ≤ d*bb s ∧ (∀ j < N s, d*bb s < Qn j s) ∧
      ‖h s-xp s‖ ≤ Real.sqrt (d*A s)) ∧
    ∃ R : Kernel S E, IsMarkovKernel R ∧ ∀ s,
      (Measurable (output (W s) 1 0) ∧ (ρ s).map (output (W s) 1 0)=R s ∧
        ρ s {ω | ∀ n, ω n ∉ accepted (W s) 1}=0 ∧
        (∫⁻ ω, 1+queryCount (W s) 1 ω ∂ρ s) ≤ ENNReal.ofReal (1+2*Real.exp 2)) ∧
      IsProbabilityMeasure (π s) ∧ π s ≪ R s ∧ R s ≪ π s ∧
      Integrable (fun x => ((π s).rnDeriv (R s) x).toReal^ell) (R s) ∧
      Integrable (fun x => ((R s).rnDeriv (π s) x).toReal^ell) (π s) ∧
      0 < (∫ x, ((π s).rnDeriv (R s) x).toReal^ell ∂R s) ∧
      0 < (∫ x, ((R s).rnDeriv (π s) x).toReal^ell ∂π s) ∧
      (∫ x, ((π s).rnDeriv (R s) x).toReal^ell ∂R s) ≤ 1+eps^2 ∧
      (∫ x, ((R s).rnDeriv (π s) x).toReal^ell ∂π s) ≤ 1+eps^2 ∧
      Real.log (∫ x, ((π s).rnDeriv (R s) x).toReal^ell ∂R s)/(ell-1) ≤ eps^2 ∧
      Real.log (∫ x, ((R s).rnDeriv (π s) x).toReal^ell ∂π s)/(ell-1) ≤ eps^2 := by
  intro bb d F T Qn N xp A g h ν W q ρ π
  have hbb : Measurable bb := hb.max measurable_const
  have hbb0 (s : S) : 0 < bb s := lt_of_lt_of_le hΘ (le_max_right _ _)
  have hbstep (s : S) : 64*(β:ℝ)*(Real.sqrt ((Module.finrank ℝ E:ℝ)*(ell+Real.log (1/eps)))+
-- Source excerpt truncated; follow the exact source link.

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/EnhancedTerminalExecution.lean:440published source at 0e31a3cda412

Excerpt truncated; the exact source link is authoritative.

No declaration docstring.

theorem enhanced_terminal_execution {V : E → ℝ} (κ : ℝ≥0) (hκ : 1 ≤ κ)
    (hV : ContDiff ℝ 2 V)
    (hH : ∀ x w : E, (κ:ℝ)⁻¹*‖w‖^2 ≤ fderiv ℝ (fderiv ℝ V) x w w ∧
      fderiv ℝ (fderiv ℝ V) x w w ≤ ‖w‖^2)
    (hd : 0 < Module.finrank ℝ E)
    (η : S → ℝ) (hηm : Measurable η) (c q Δ : ℝ)
    (hc : 0 < c) (hc1 : c < 1/4) (hη : ∀ s, 0 < η s ∧ η s ≤ c)
    (hq : 2 ≤ q) (hΔ : 0 < Δ) (hΔ1 : Δ ≤ 1/2)
    (s0 : S) (M : Kernel S E) [IsMarkovKernel M] :
    let d : ℝ := Module.finrank ℝ E
    let K := fun r : ℝ => (1+r)/((κ:ℝ)⁻¹+r)
    let τ := fun s : S => if 2 ≤ K s.1 then K s.1 else c
    let L := q+Real.log (K s0.1*d*q/Δ)
    let B := (1/1024)/(Real.sqrt (d*L)+L)
    let J := Nat.ceil (8*Real.log (Real.exp 1*K s0.1/B))
    let eps := Δ/Real.sqrt ((J:ℝ)+1)
    let v := fun s : S => (η s+τ s)/(1+s.1)
    let bp := fun s : S => s.1+Real.toNNReal (v s)⁻¹
    let obs := fun p : S × (E × E) => p.2.1+Real.sqrt (τ p.1/(1+p.1.1)) • p.2.2
    let center := fun p : S × (E × E) => (bp p.1:ℝ)⁻¹ •
      ((p.1.1:ℝ) • p.1.2.1+(v p.1)⁻¹ • obs p)
    let F := fun p x => V x+(bp p.1:ℝ)/2*‖x-center p‖^2
    let T := fun p x => x-(1+(bp p.1:ℝ))⁻¹ • gradient (F p) x
    let Qn := fun n p => ‖gradient (F p) ((T p)^[n] p.2.1)‖^2
    let N := firstIndex Qn (fun p => ((κ:ℝ)⁻¹+bp p.1)*d)
    let update := fun p : S × (E × E) =>
      (bp p.1,center p,(T p)^[N p] p.2.1,p.1.2.2.2.1+1,
        fun n => Nat.casesOn n (obs p) p.1.2.2.2.2)
    let bb := fun s : S => max (s.1:ℝ) B⁻¹
    let TF := fun s x => V x+bb s/2*‖x-s.2.1‖^2
    let TT := fun s x => x-(1+bb s)⁻¹ • gradient (TF s) x
    let TQn := fun n s => ‖gradient (TF s) ((TT s)^[n] (s.2.2.1))‖^2
    let TN := firstIndex TQn (fun s => d*bb s)
    let txp := fun s => (TT s)^[TN s] (s.2.2.1)
    let TA := fun s => (bb s)⁻¹
    let tg := fun s => gradient V (txp s)
    let th := fun s => s.2.1-TA s • tg s
    let tnu := (volume.restrict (Ioc (0 : ℝ) 1)).prod (stdGaussian E)
    let tW := fun s (p : E × (ℝ × E)) => min 1 (max (-1) (scaledCached V TA th tg ((s,p.1),p.2)))
    let tproposal := fun s => (stdGaussian E).map (fun z => th s+Real.sqrt (TA s) • z)
    let trho := fun s => Measure.infinitePi (fun _ : ℕ => attemptLaw (tproposal s) tnu 1 (by norm_num))
    let tpi := fun s => ((stdGaussian E).map (fun z => s.2.1+Real.sqrt (TA s) • z)).tilted (fun x => -V x)
    0 < B ∧ 0 < eps ∧ eps ≤ 1/2 ∧
    (∀ s : S, B⁻¹ ≤ (s.1:ℝ) → bb s=(s.1:ℝ)) ∧
-- Source excerpt truncated; follow the exact source link.

AutoSamplingTheory/ExampleCases/SmoothedPicardHMC/EnhancedTerminalExecution.lean:494published source at 0e31a3cda412

Excerpt truncated; the exact source link is authoritative.