diff --git a/src/iwp.jl b/src/iwp.jl index 817caa9..59d9d39 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) = @. 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 index cf7e475..b19d3f3 100644 --- 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