diff --git a/src/iwp.jl b/src/iwp.jl index ac3c567..4786291 100644 --- 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) = Diagonal(@. 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} @@ -110,6 +124,20 @@ function transition_cov_cholf_1d(ndiff::Integer, dt::T, ::ReverseTaylor) where { return L end + +function transition_cov_cholf_precond_1d(ndiff::Integer, T, ::ReverseTaylor) + L = zeros(T, ndiff + one(ndiff), ndiff + one(ndiff)) + @simd ivdep for n = 0:ndiff + @simd ivdep for m = 0:ndiff + if m <= n + @inbounds L[n+1, m+1] = + sqrt(2 * m + 1) * factorial(n)^2 / factorial(n - m) / factorial(n + m + 1) + end + end + end + return L +end + """ transition_cov(model::IWP{T}, dt::T, method::AbstractRealizationMethod) where {T} diff --git a/test/runtests.jl b/test/runtests.jl index e55a898..7527e31 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -31,5 +31,18 @@ using Test end + @testset "preconditioning" begin + for ndiff = 0:9 + Φ = IntegratedWienerProcesses.transition_matrix_1d(ndiff, dt, ReverseTaylor()) + L = IntegratedWienerProcesses.transition_cov_cholf_1d(ndiff, dt, ReverseTaylor()) + 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 Φ ≈ prec * Φbreve / prec + @test L ≈ prec * Lbreve + end + end + end