Skip to content

Commit

Permalink
Run JuliaFormatter.jl
Browse files Browse the repository at this point in the history
  • Loading branch information
nathanaelbosch committed Oct 20, 2023
1 parent 9d2b481 commit e30fc44
Show file tree
Hide file tree
Showing 16 changed files with 161 additions and 63 deletions.
14 changes: 9 additions & 5 deletions src/caches.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,8 @@
########################################################################################
mutable struct EKCache{
RType,ProjType,SolProjType,PType,PIType,EType,uType,duType,xType,PriorType,AType,QType,
HType,matType,bkType,diffusionType,diffModelType,measModType,measType,puType,llType,dtType,
rateType,UF,JC,uNoUnitsType,
HType,matType,bkType,diffusionType,diffModelType,measModType,measType,puType,llType,
dtType,rateType,UF,JC,uNoUnitsType,
} <: AbstractODEFilterCache
# Constants
d::Int # Dimension of the problem
Expand Down Expand Up @@ -93,7 +93,9 @@ function OrdinaryDiffEq.alg_cache(

FAC = get_covariance_factorization(alg)
if FAC isa KroneckerCovariance && !(f.mass_matrix isa UniformScaling)
error("The selected algorithm uses an efficient Kronecker-factorized implementation which is incompatible with the provided mass matrix. Try using the `EK1` instead.")
error(
"The selected algorithm uses an efficient Kronecker-factorized implementation which is incompatible with the provided mass matrix. Try using the `EK1` instead.",
)
end

uType = typeof(u)
Expand Down Expand Up @@ -136,7 +138,9 @@ function OrdinaryDiffEq.alg_cache(
# Initial State
initial_variance = ones(uElType, q + 1)
μ0 = zeros(uElType, D)
Σ0 = PSDMatrix(to_factorized_matrix(FAC, IsoKroneckerProduct(d, diagm(sqrt.(initial_variance)))))
Σ0 = PSDMatrix(
to_factorized_matrix(FAC, IsoKroneckerProduct(d, diagm(sqrt.(initial_variance)))),
)
x0 = Gaussian(μ0, Σ0)

# Diffusion Model
Expand Down Expand Up @@ -173,7 +177,7 @@ function OrdinaryDiffEq.alg_cache(
backward_kernel = AffineNormalKernel(
factorized_zeros(FAC, uElType, D, D; d, q),
zeros(uElType, D),
PSDMatrix(factorized_zeros(FAC, uElType, 2D, D; d, q))
PSDMatrix(factorized_zeros(FAC, uElType, 2D, D; d, q)),
)

u_pred = copy(u)
Expand Down
1 change: 0 additions & 1 deletion src/diffusions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,6 @@ function estimate_global_diffusion(::FixedDiffusion, integ)
end
diffusion_t = dot(v, e) / d


if integ.success_iter == 0
# @assert length(sol_diffusions) == 0
global_diffusion = diffusion_t
Expand Down
1 change: 0 additions & 1 deletion src/fast_linalg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,6 @@ _matmul!(
B::AbstractVecOrMat{T},
) where {T<:LinearAlgebra.BlasFloat} = matmul!(C, A, B)


"""
getupperright!(A)
Expand Down
18 changes: 15 additions & 3 deletions src/filtering/markov_kernel.jl
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,11 @@ function marginalize!(xout, x, K; C_DxD, C_3DxD)
marginalize_mean!(xout.μ, x.μ, K)
marginalize_cov!(xout.Σ, x.Σ, K; C_DxD, C_3DxD)
end
function marginalize_mean!(μout::AbstractVecOrMat, μ::AbstractVecOrMat, K::AffineNormalKernel)
function marginalize_mean!(
μout::AbstractVecOrMat,
μ::AbstractVecOrMat,
K::AffineNormalKernel,
)
_matmul!(μout, K.A, μ)
if !ismissing(K.b)
μout .+= K.b
Expand Down Expand Up @@ -224,11 +228,19 @@ function compute_backward_kernel!(
D = full_state_dim = length(x.μ)
d = ode_dimension_dim = K.A.ldim
Q = n_derivatives_dim = D ÷ d
_Kout = AffineNormalKernel(Kout.A.B, reshape_no_alloc(Kout.b, Q, d), PSDMatrix(Kout.C.R.B))
_Kout =
AffineNormalKernel(Kout.A.B, reshape_no_alloc(Kout.b, Q, d), PSDMatrix(Kout.C.R.B))
_x_pred = Gaussian(reshape_no_alloc(xpred.μ, Q, d), PSDMatrix(xpred.Σ.R.B))
_x = Gaussian(reshape_no_alloc(x.μ, Q, d), PSDMatrix(x.Σ.R.B))
_K = AffineNormalKernel(K.A.B, reshape_no_alloc(K.b, Q, d), PSDMatrix(K.C.R.B))
_D = size(_Kout.A, 1)
_C_DxD = view(C_DxD, 1:_D, 1:_D)
return compute_backward_kernel!(_Kout, _x_pred, _x, _K; C_DxD=_C_DxD, diffusion=diffusion)
return compute_backward_kernel!(
_Kout,
_x_pred,
_x,
_K;
C_DxD=_C_DxD,
diffusion=diffusion,
)
end
9 changes: 7 additions & 2 deletions src/filtering/predict.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,8 @@ predict(x::Gaussian, A::AbstractMatrix, Q::AbstractMatrix) =
Gaussian(predict_mean(x.μ, A), predict_cov(x.Σ, A, Q))
predict_mean::AbstractVector, A::AbstractMatrix) = A * μ
predict_cov::AbstractMatrix, A::AbstractMatrix, Q::AbstractMatrix) = A * Σ * A' + Q
predict_cov::PSDMatrix, A::AbstractMatrix, Q::PSDMatrix) = PSDMatrix(qr([Σ.R * A'; Q.R]).R)
predict_cov::PSDMatrix, A::AbstractMatrix, Q::PSDMatrix) =
PSDMatrix(qr([Σ.R * A'; Q.R]).R)
predict_cov::PSDMatrix{T,<:IKP}, A::IKP, Q::PSDMatrix{T,<:IKP}) where {T} = begin
P_pred_breve = predict_cov(PSDMatrix.R.B), A.B, PSDMatrix(Q.R.B))
return PSDMatrix(IsoKroneckerProduct.R.ldim, P_pred_breve.R))
Expand Down Expand Up @@ -45,7 +46,11 @@ function predict!(
return x_out
end

function predict_mean!(m_out::AbstractVecOrMat, m_curr::AbstractVecOrMat, Ah::AbstractMatrix)
function predict_mean!(
m_out::AbstractVecOrMat,
m_curr::AbstractVecOrMat,
Ah::AbstractMatrix,
)
_matmul!(m_out, Ah, m_curr)
return m_out
end
Expand Down
5 changes: 4 additions & 1 deletion src/filtering/smooth.jl
Original file line number Diff line number Diff line change
Expand Up @@ -144,7 +144,10 @@ function smooth!(
C_DxD=view(cache.C_DxD, 1:_D, 1:_D),
C_2DxD=view(cache.C_2DxD, 1:2*_D, 1:_D),
C_3DxD=view(cache.C_3DxD, 1:3*_D, 1:_D),
x_pred=Gaussian(reshape_no_alloc(cache.x_pred.μ, Q, d), PSDMatrix(cache.x_pred.Σ.R.B)),
x_pred=Gaussian(
reshape_no_alloc(cache.x_pred.μ, Q, d),
PSDMatrix(cache.x_pred.Σ.R.B),
),
)
smooth!(
_x_curr,
Expand Down
2 changes: 1 addition & 1 deletion src/gaussians.jl
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ std(g::Gaussian) = sqrt.(diag(g.Σ))
############################################################################################
# `SRGaussian`: Gaussians with PDFMatrix covariances
############################################################################################
const SRGaussian{T,S} = Gaussian{VM,PSDMatrix{T,S}} where {VM <: AbstractVecOrMat{T}}
const SRGaussian{T,S} = Gaussian{VM,PSDMatrix{T,S}} where {VM<:AbstractVecOrMat{T}}
Base.:*(M::AbstractMatrix, g::SRGaussian) = Gaussian(M * g.μ, X_A_Xt(g.Σ, M))
# GaussianDistributions.whiten(Σ::PSDMatrix, z) = Σ.L\z

Expand Down
2 changes: 1 addition & 1 deletion src/initialization/taylormode.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ function initial_update!(integ, cache, init::TaylorModeInit)
# This is hacky and should definitely be removed. But it also works so 🤷
MM = if f.mass_matrix isa UniformScaling
f.mass_matrix
else
else
_MM = copy(f.mass_matrix)
if any(iszero.(diag(_MM)))
_MM = typeof(promote(_MM[1], 1e-20)[1]).(_MM)
Expand Down
75 changes: 59 additions & 16 deletions src/kronecker.jl
Original file line number Diff line number Diff line change
@@ -1,11 +1,13 @@
mutable struct IsoKroneckerProduct{T<:Number,TB<:AbstractMatrix} <: Kronecker.AbstractKroneckerProduct{T}
mutable struct IsoKroneckerProduct{T<:Number,TB<:AbstractMatrix} <:
Kronecker.AbstractKroneckerProduct{T}
ldim::Int64
B::TB
function IsoKroneckerProduct(ldim::Int64, B::AbstractMatrix{T}) where {T}
return new{T,typeof(B)}(ldim, B)
end
end
IsoKroneckerProduct(ldim::Integer, B::AbstractVector) = IsoKroneckerProduct(ldim, reshape(B, :, 1))
IsoKroneckerProduct(ldim::Integer, B::AbstractVector) =
IsoKroneckerProduct(ldim, reshape(B, :, 1))
const IKP = IsoKroneckerProduct

Kronecker.getmatrices(K::IKP) = (I(K.ldim), K.B)
Expand Down Expand Up @@ -58,8 +60,13 @@ _matmul!(A::IKP, B::IKP, C::IKP, alpha::Number, beta::Number) = begin
_matmul!(A.B, B.B, C.B)
return A
end
_matmul!(A::IKP{T}, B::IKP{T}, C::IKP{T}, alpha::Number, beta::Number
) where {T<:LinearAlgebra.BlasFloat} = begin
_matmul!(
A::IKP{T},
B::IKP{T},
C::IKP{T},
alpha::Number,
beta::Number,
) where {T<:LinearAlgebra.BlasFloat} = begin
@assert A.ldim == B.ldim == C.ldim
_matmul!(A.B, B.B, C.B, alpha, beta)
return A
Expand Down Expand Up @@ -107,7 +114,7 @@ function mul_vectrick!(
v::AbstractVecOrMat,
alpha::Number,
beta::Number,
)
)
N = A.B
c, d = size(N)

Expand All @@ -117,18 +124,54 @@ function mul_vectrick!(
return x
end

_matmul!(C::AbstractVecOrMat, A::IsoKroneckerProduct, B::AbstractVecOrMat) = mul_vectrick!(C, A, B)
_matmul!(C::AbstractVecOrMat, A::IsoKroneckerProduct, B::AbstractVecOrMat) =
mul_vectrick!(C, A, B)
mul!(C::AbstractMatrix, A::IsoKroneckerProduct, B::AbstractMatrix) = mul_vectrick!(C, A, B)
mul!(C::AbstractMatrix, A::IsoKroneckerProduct, B::Adjoint{T, <:AbstractMatrix{T}}) where {T} = mul_vectrick!(C, A, B)
mul!(
C::AbstractMatrix,
A::IsoKroneckerProduct,
B::Adjoint{T,<:AbstractMatrix{T}},
) where {T} = mul_vectrick!(C, A, B)
mul!(C::AbstractVector, A::IsoKroneckerProduct, B::AbstractVector) = mul_vectrick!(C, A, B)

_matmul!(C::AbstractVecOrMat{T}, A::IsoKroneckerProduct{T}, B::AbstractVecOrMat{T}) where {T<:LinearAlgebra.BlasFloat} = mul_vectrick!(C, A, B)
_matmul!(C::AbstractVecOrMat, A::AbstractVecOrMat, B::IsoKroneckerProduct) = _matmul!(C', B', A')
_matmul!(C::AbstractVecOrMat{T}, A::AbstractVecOrMat{T}, B::IsoKroneckerProduct{T}) where {T<:LinearAlgebra.BlasFloat} = _matmul!(C', B', A')
_matmul!(
C::AbstractVecOrMat{T},
A::IsoKroneckerProduct{T},
B::AbstractVecOrMat{T},
) where {T<:LinearAlgebra.BlasFloat} = mul_vectrick!(C, A, B)
_matmul!(C::AbstractVecOrMat, A::AbstractVecOrMat, B::IsoKroneckerProduct) =
_matmul!(C', B', A')
_matmul!(
C::AbstractVecOrMat{T},
A::AbstractVecOrMat{T},
B::IsoKroneckerProduct{T},
) where {T<:LinearAlgebra.BlasFloat} = _matmul!(C', B', A')

_matmul!(C::AbstractVecOrMat, A::IsoKroneckerProduct, B::AbstractVecOrMat, alpha::Number, beta::Number) = mul_vectrick!(C, A, B, alpha, beta)
_matmul!(C::AbstractVecOrMat{T}, A::IsoKroneckerProduct{T}, B::AbstractVecOrMat{T}, alpha::Number, beta::Number
) where {T<:LinearAlgebra.BlasFloat} = mul_vectrick!(C, A, B, alpha, beta)
_matmul!(C::AbstractVecOrMat, A::AbstractVecOrMat, B::IsoKroneckerProduct, alpha::Number, beta::Number) = mul_vectrick!(C', B', A', alpha, beta)
_matmul!(C::AbstractVecOrMat{T}, A::AbstractVecOrMat{T}, B::IsoKroneckerProduct{T}, alpha::Number, beta::Number
) where {T<:LinearAlgebra.BlasFloat} = mul_vectrick!(C', B', A', alpha, beta)
_matmul!(
C::AbstractVecOrMat,
A::IsoKroneckerProduct,
B::AbstractVecOrMat,
alpha::Number,
beta::Number,
) = mul_vectrick!(C, A, B, alpha, beta)
_matmul!(
C::AbstractVecOrMat{T},
A::IsoKroneckerProduct{T},
B::AbstractVecOrMat{T},
alpha::Number,
beta::Number,
) where {T<:LinearAlgebra.BlasFloat} = mul_vectrick!(C, A, B, alpha, beta)
_matmul!(
C::AbstractVecOrMat,
A::AbstractVecOrMat,
B::IsoKroneckerProduct,
alpha::Number,
beta::Number,
) = mul_vectrick!(C', B', A', alpha, beta)
_matmul!(
C::AbstractVecOrMat{T},
A::AbstractVecOrMat{T},
B::IsoKroneckerProduct{T},
alpha::Number,
beta::Number,
) where {T<:LinearAlgebra.BlasFloat} = mul_vectrick!(C', B', A', alpha, beta)
7 changes: 4 additions & 3 deletions src/perform_step.jl
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,8 @@ function OrdinaryDiffEq.perform_step!(integ, cache::EKCache, repeat_step=false)

# Predict the mean
predict_mean!(x_pred.μ, xprev.μ, Ah)
write_into_solution!(integ.u, x_pred.μ; cache, is_secondorder_ode=integ.f isa DynamicalODEFunction)
write_into_solution!(
integ.u, x_pred.μ; cache, is_secondorder_ode=integ.f isa DynamicalODEFunction)

# Measure
evaluate_ode!(integ, x_pred, tnew)
Expand Down Expand Up @@ -116,7 +117,8 @@ function OrdinaryDiffEq.perform_step!(integ, cache::EKCache, repeat_step=false)

# Update state and save the ODE solution value
x_filt = update!(cache, x_pred)
write_into_solution!(integ.u, x_filt.μ; cache, is_secondorder_ode=integ.f isa DynamicalODEFunction)
write_into_solution!(
integ.u, x_filt.μ; cache, is_secondorder_ode=integ.f isa DynamicalODEFunction)

# Update the global diffusion MLE (if applicable)
if !isdynamic(cache.diffusionmodel)
Expand Down Expand Up @@ -234,7 +236,6 @@ function estimate_errors!(cache::AbstractODEFilterCache)
error_estimate .= sqrt.(error_estimate)
return error_estimate
elseif local_diffusion isa Number

_matmul!(R, Qh.R, H')

# error_estimate = diag(PSDMatrix(R))
Expand Down
12 changes: 10 additions & 2 deletions src/priors/common.jl
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,11 @@ See also: [`make_transition_matrices`](@ref).
[1] N. Krämer, P. Hennig: **Stable Implementation of Probabilistic ODE Solvers** (2020)
"""
function initialize_transition_matrices(::DenseCovariance, p::AbstractODEFilterPrior{T}, dt) where {T}
function initialize_transition_matrices(
::DenseCovariance,
p::AbstractODEFilterPrior{T},
dt,
) where {T}
d, q = p.wiener_process_dimension, p.num_derivatives
D = d * (q + 1)
Ah, Qh = zeros(T, D, D), PSDMatrix(zeros(T, D, D))
Expand All @@ -40,7 +44,11 @@ function initialize_transition_matrices(::DenseCovariance, p::AbstractODEFilterP
Q = copy(Qh)
return A, Q, Ah, Qh, P, PI
end
function initialize_transition_matrices(fac::CovarianceFactorization, p::AbstractODEFilterPrior, dt)
function initialize_transition_matrices(
fac::CovarianceFactorization,
p::AbstractODEFilterPrior,
dt,
)
error("The chosen prior can not be implemented with a $fac factorization")
end

Expand Down
9 changes: 8 additions & 1 deletion src/projection.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
function projection(
::DenseCovariance,
d::Integer,
q::Integer,
::Type{elType}=typeof(1.0),
Expand All @@ -14,6 +13,14 @@ function projection(
end
return Proj
end
function projection(
::DenseCovariance,
d::Integer,
q::Integer,
::Type{elType}=typeof(1.0),
) where {elType}
projection(d, q, elType)
end

function projection(
::KroneckerCovariance,
Expand Down
25 changes: 19 additions & 6 deletions test/core/filtering.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ import ProbNumDiffEq: IsoKroneckerProduct
x_out = copy(x_curr)

_fstr(F) = F ? "Kronecker" : "None"
@testset "Factorization: $(_fstr(KRONECKER))"for KRONECKER in (false, true)
@testset "Factorization: $(_fstr(KRONECKER))" for KRONECKER in (false, true)
if KRONECKER
K = 2
m = kron(ones(K), m)
Expand Down Expand Up @@ -78,7 +78,13 @@ import ProbNumDiffEq: IsoKroneckerProduct
PSDMatrix([Q_R; zero(Q_R)])
end
K = ProbNumDiffEq.AffineNormalKernel(A, Q_SR)
ProbNumDiffEq.marginalize!(x_out, x_curr, K; C_DxD=zeros(d, d), C_3DxD=zeros(3d, d))
ProbNumDiffEq.marginalize!(
x_out,
x_curr,
K;
C_DxD=zeros(d, d),
C_3DxD=zeros(3d, d),
)
@test m_p x_out.μ
@test P_p Matrix(x_out.Σ)
end
Expand Down Expand Up @@ -298,7 +304,7 @@ end
@test P_smoothed Matrix(x_out.Σ)
end
@testset "smooth!" begin
_d = KRONECKER ? K*d : d
_d = KRONECKER ? K * d : d
x_curr_psd = Gaussian(m, PSDMatrix(P_R)) |> copy
x_next_psd = Gaussian(m_s, PSDMatrix(P_s_R)) |> copy
cache = (
Expand All @@ -322,7 +328,8 @@ end
@testset "smooth via backward kernels" begin
K_forward = ProbNumDiffEq.AffineNormalKernel(copy(A), copy(Q_SR))
K_backward = ProbNumDiffEq.AffineNormalKernel(
copy(A), copy(m_p), if KRONECKER
copy(A), copy(m_p),
if KRONECKER
PSDMatrix(IsoKroneckerProduct(K, zeros(2d, d)))
else
PSDMatrix(zeros(2d, d))
Expand All @@ -333,7 +340,7 @@ end
x_next_smoothed = Gaussian(m_s, PSDMatrix(P_s_R)) |> copy

C_DxD = if KRONECKER
zeros(K*d, K*d)
zeros(K * d, K * d)
else
zeros(d, d)
end
Expand All @@ -349,7 +356,13 @@ end

C_3DxD = zeros(3d, d)
ProbNumDiffEq.marginalize_mean!(x_curr.μ, x_next_smoothed.μ, K_backward)
ProbNumDiffEq.marginalize_cov!(x_curr.Σ, x_next_smoothed.Σ, K_backward; C_DxD, C_3DxD)
ProbNumDiffEq.marginalize_cov!(
x_curr.Σ,
x_next_smoothed.Σ,
K_backward;
C_DxD,
C_3DxD,
)

@test m_smoothed x_curr.μ
@test P_smoothed Matrix(x_curr.Σ)
Expand Down
Loading

0 comments on commit e30fc44

Please sign in to comment.