preconditioned transition matrix

filtron · Oct 13, 2023 · 5ff2c04 · 5ff2c04
1 parent 7c3db97
commit 5ff2c04
Show file tree

Hide file tree

Showing 2 changed files with 25 additions and 7 deletions.
diff --git a/src/iwp.jl b/src/iwp.jl
@@ -64,10 +64,10 @@ function transition_matrix(
     method::AbstractRealizationMethod,
 ) where {T}
     V = promote_type(real(T), typeof(dt))
-    return kron(_transition_matrix_1d(ndiff(model), V(dt), method), one(T)I(dim(model)))
+    return kron(transition_matrix_1d(ndiff(model), V(dt), method), one(T)I(dim(model)))
 end
 
-function _transition_matrix_1d(ndiff::Integer, dt::T, ::ReverseTaylor) where {T<:Real}
+function transition_matrix_1d(ndiff::Integer, dt::T, ::ReverseTaylor) where {T<:Real}
     irange = collect(0:ndiff)
     #g = @. δt^irange / factorial(irange)
     g = @. exp(irange * log(dt) - logfactorial(irange))
@@ -82,6 +82,20 @@ function _transition_matrix_1d(ndiff::Integer, dt::T, ::ReverseTaylor) where {T<
 end
 
 
+precond(ndiff::Integer, dt::Real, ::ReverseTaylor) =  @. sqrt(dt) * dt^(0:ndiff) / factorial(0:ndiff)
+
+function transition_matrix_precond_1d(ndiff::Integer, T, ::ReverseTaylor)
+    Φ = zeros(T, ndiff + one(ndiff), ndiff + one(ndiff))
+    @simd ivdep for col = 0:ndiff
+        @simd ivdep for row = col:ndiff
+            @inbounds Φ[row+1, col+1] = binomial(row, col)
+        end
+    end
+    return Φ
+end
+
+
+
 """
     transition_cov_cholf(model::IWP{T}, dt::T, method::AbstractRealizationMethod) where {T}
 
@@ -111,9 +125,8 @@ function transition_cov_cholf_1d(ndiff::Integer, dt::T, ::ReverseTaylor) where {
 end
 
 
-function transition_cov_cholf_precond_1d(ndiff::Integer, dt::T, ::ReverseTaylor) where {T}
+function transition_cov_cholf_precond_1d(ndiff::Integer, T, ::ReverseTaylor)
     L = zeros(T, ndiff + one(ndiff), ndiff + one(ndiff))
-    precond = @. sqrt(dt) * dt^(0:ndiff) / factorial(0:ndiff)
     @simd ivdep for n = 0:ndiff
         @simd ivdep for m = 0:ndiff
             if m <= n
@@ -122,7 +135,7 @@ function transition_cov_cholf_precond_1d(ndiff::Integer, dt::T, ::ReverseTaylor)
             end
         end
     end
-    return precond, L
+    return L
 end
 
 """

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -33,9 +33,14 @@ using Test
 
     @testset "preconditioning" begin 
         for ndiff = 0:9
+            Φ = IntegratedWienerProcesses.transition_matrix_1d(ndiff, dt, ReverseTaylor())
             L = IntegratedWienerProcesses.transition_cov_cholf_1d(ndiff, dt, ReverseTaylor())
-            precond, Lbreve = IntegratedWienerProcesses.transition_cov_cholf_precond_1d(ndiff, dt, ReverseTaylor())
-            @test L ≈ Diagonal(precond) * Lbreve
+            prec = IntegratedWienerProcesses.precond(ndiff, dt, ReverseTaylor())
+            Φbreve = IntegratedWienerProcesses.transition_matrix_precond_1d(ndiff, Float64, ReverseTaylor()) 
+            Lbreve = IntegratedWienerProcesses.transition_cov_cholf_precond_1d(ndiff, Float64, ReverseTaylor())
+
+            @test Φ ≈ Diagonal(prec) * Φbreve / Diagonal(prec)
+            @test L ≈ Diagonal(prec) * Lbreve
         end
     end