diff --git a/src/kalman_base.jl b/src/kalman_base.jl index fd1f0b6aae2a7f404e8f031db0a882d9ed44090b..7b5b68baac250ecc79a16a00d0d91e6a30757e9e 100644 --- a/src/kalman_base.jl +++ b/src/kalman_base.jl @@ -1,4 +1,4 @@ -function det_from_cholesky(achol::AbstractMatrix{T}) where T <: AbstractFloat +function det_from_cholesky(achol::AbstractArray{T}) where T <: AbstractFloat x = 1.0 @inbounds @simd for i = 1:size(achol,1) x *= achol[i,i] @@ -14,7 +14,7 @@ end # alphah_t = a_t + Pstar_t*r0_{t-1} + Pinf_t*r1_{t-1} (DK 5.24) function get_alphah!(valphah::AbstractVector{T}, va::AbstractVector{T}, - vPstar::AbstractMatrix{T}, vPinf::AbstractMatrix{T}, + vPstar::AbstractArray{T}, vPinf::AbstractArray{T}, r0::AbstractVector{T}, r1::AbstractVector{T}) where T <: AbstractFloat copy!(valphah, va) mul!(valphah, vPstar, r0, 1.0, 1.0) @@ -45,21 +45,21 @@ end # D = inv(F_t) + K(DK)_t'*N_t*K(DK)_t (DK 4.69) -function get_D!(D::AbstractMatrix{U}, iF::AbstractMatrix{U}, K::AbstractMatrix{U}, N::AbstractMatrix{U}, tmp::AbstractMatrix{U}) where U <: AbstractFloat +function get_D!(D::AbstractArray{U}, iF::AbstractArray{U}, K::AbstractArray{U}, N::AbstractArray{U}, tmp::AbstractArray{U}) where U <: AbstractFloat copy!(D,iF) mul!(tmp, transpose(K), N) mul!(D, tmp, K, 1.0, 1.0) end # D_t = KDKinf_t'*N0_t*KDKinf_t (DK p. 135) -function get_D!(D::AbstractMatrix{T}, KDKinf::AbstractMatrix{T}, N0::AbstractMatrix{T}, Tmp::AbstractMatrix{T}) where T <: AbstractFloat +function get_D!(D::AbstractArray{T}, KDKinf::AbstractArray{T}, N0::AbstractArray{T}, Tmp::AbstractArray{T}) where T <: AbstractFloat mul!(Tmp, Transpose(KDKinf), N0) mul!(D, Tmp, KDKinf) end # epsilon_h = H*(iF_t*v_t - K(DK)_t'*r_t) (DK 4.69) -function get_epsilonh!(epsilon::AbstractVector{U}, H::AbstractMatrix{U}, - iFv::AbstractVector{U}, K::AbstractMatrix{U}, +function get_epsilonh!(epsilon::AbstractVector{U}, H::AbstractArray{U}, + iFv::AbstractVector{U}, K::AbstractArray{U}, r::AbstractVector{U}, tmp1::AbstractVector{U}, tmp2::AbstractVector{U}) where U <: AbstractFloat @@ -69,9 +69,9 @@ function get_epsilonh!(epsilon::AbstractVector{U}, H::AbstractMatrix{U}, end # epsilonh_t = H*(iF_t*v_t - K_t*T*r_t) Different K !!! -function get_epsilonh!(epsilon::AbstractVector{U}, H::AbstractMatrix{U}, - iFv::AbstractVector{U}, K::AbstractMatrix{U}, - T::AbstractMatrix{U}, r::AbstractVector{U}, +function get_epsilonh!(epsilon::AbstractVector{U}, H::AbstractArray{U}, + iFv::AbstractVector{U}, K::AbstractArray{U}, + T::AbstractArray{U}, r::AbstractVector{U}, tmp1::AbstractVector{U}, tmp2::AbstractVector{U}) where U <: AbstractFloat copy!(tmp1, iFv) @@ -81,16 +81,16 @@ function get_epsilonh!(epsilon::AbstractVector{U}, H::AbstractMatrix{U}, end # epsilon_t = -H_t*KDKinf'*r0_t (DK p. 135) -function get_epsilonh!(epsilon::AbstractVector{T}, H::AbstractMatrix{T}, - KDKinf::AbstractMatrix{T}, r0::AbstractVector{T}, +function get_epsilonh!(epsilon::AbstractVector{T}, H::AbstractArray{T}, + KDKinf::AbstractArray{T}, r0::AbstractVector{T}, tmp::AbstractVector{T}) where T <: AbstractFloat mul!(tmp, transpose(KDKinf), r0) mul!(epsilon, H, tmp, -1.0, 0.0) end # etah = Q*R'*r_t -function get_etah!(etah::AbstractVector{T}, Q::AbstractMatrix{T}, - R::AbstractMatrix{T}, r::AbstractVector{T}, +function get_etah!(etah::AbstractVector{T}, Q::AbstractArray{T}, + R::AbstractArray{T}, r::AbstractVector{T}, tmp::AbstractVector{T}) where T <: AbstractFloat mul!(tmp, transpose(R), r) mul!(etah, Q, tmp) @@ -102,25 +102,25 @@ function get_F(Zi, P, h, tmp) return F end -function get_F!(f::AbstractArray{T}, zp::AbstractArray{T}, z::AbstractArray{T}, p::AbstractArray{T}) where T <: AbstractFloat +function get_F!(f::AbstractMatrix{T}, zp::AbstractArray{T}, z::AbstractArray{T}, p::AbstractMatrix{T}) where T <: AbstractFloat mul!(zp, z, p) mul!(f, zp, transpose(z)) end -function get_F!(f::AbstractArray{T}, zp::AbstractArray{T}, z::AbstractVector{U}, p::AbstractArray{T}) where {T <: AbstractFloat, U <: Integer} +function get_F!(f::AbstractMatrix{T}, zp::AbstractArray{T}, z::AbstractVector{U}, p::AbstractMatrix{T}) where {T <: AbstractFloat, U <: Integer} zp .= view(p, z, :) f .= view(zp, :, z) end # F = Z*P*Z' + H -function get_F!(f::AbstractArray{T}, zp::AbstractArray{T}, z::AbstractArray{T}, p::AbstractArray{T}, h::AbstractArray{T}) where T <: AbstractFloat +function get_F!(f::AbstractMatrix{T}, zp::AbstractArray{T}, z::AbstractArray{T}, p::AbstractMatrix{T}, h::AbstractMatrix{T}) where T <: AbstractFloat copy!(f, h) mul!(zp, z, p) gemm!('N', 'T', 1.0, zp, z, 1.0, f) end # F = P(z,z) + H -function get_F!(f::AbstractArray{T}, zp::AbstractArray{T}, z::AbstractVector{I}, p::AbstractArray{T}, h::AbstractArray{T}) where {I <: Integer, T <: AbstractFloat} +function get_F!(f::AbstractMatrix{T}, zp::AbstractArray{T}, z::AbstractVector{I}, p::AbstractMatrix{T}, h::AbstractMatrix{T}) where {I <: Integer, T <: AbstractFloat} copy!(f, h) zp .= view(p, z, :) f .+= view(zp, :, z) @@ -240,11 +240,11 @@ function get_L!(L::AbstractArray{U}, T::AbstractArray{U}, K::AbstractArray{U}, z end # L1_t = - KDK*Z (DK 5.12) -function get_L1!(L1::AbstractMatrix{T}, KDK::AbstractMatrix{T}, Z::AbstractMatrix{T}) where T <: AbstractFloat +function get_L1!(L1::AbstractArray{T}, KDK::AbstractArray{T}, Z::AbstractArray{T}) where T <: AbstractFloat mul!(L1, KDK, Z, -1.0, 0.0) end -function get_L1!(L1::AbstractMatrix{T}, KDK::AbstractMatrix{T}, z::AbstractVector{U}) where {T <: AbstractFloat, U <: Integer} +function get_L1!(L1::AbstractArray{T}, KDK::AbstractArray{T}, z::AbstractVector{U}) where {T <: AbstractFloat, U <: Integer} vL1 = view(L1, :, z) vL1 .= -KDK end @@ -273,40 +273,40 @@ function get_prediction_error(Y::AbstractArray{T}, z::AbstractVector{U}, a::Abst prediction_error = Y[i, t] - a[z[i]] # nu_{t,i} in 6.13 in DK (2012) end -function get_QQ!(c::AbstractMatrix{T}, a::AbstractMatrix{T}, b::AbstractMatrix{T}, work::Matrix{T}) where T <: AbstractFloat +function get_QQ!(c::AbstractArray{T}, a::AbstractArray{T}, b::AbstractArray{T}, work::Array{T}) where T <: AbstractFloat mul!(work, a, b) mul!(c, work, transpose(a)) end # att = a + K'*v -function get_updated_a!(att::AbstractVector{T}, a::AbstractVector{T}, K::AbstractMatrix{T}, v::AbstractVector{T}) where T <: AbstractFloat +function get_updated_a!(att::AbstractVector{T}, a::AbstractVector{T}, K::AbstractArray{T}, v::AbstractVector{T}) where T <: AbstractFloat copy!(att, a) mul!(att, transpose(K), v, 1.0, 1.0) end # Ptt = P - K'*Z*P -function get_updated_Ptt!(Ptt::AbstractMatrix{T}, P::AbstractMatrix{T}, K::AbstractMatrix{T}, ZP::AbstractMatrix{T}) where T <: AbstractFloat +function get_updated_Ptt!(Ptt::AbstractArray{T}, P::AbstractArray{T}, K::AbstractArray{T}, ZP::AbstractArray{T}) where T <: AbstractFloat copy!(Ptt, P) mul!(Ptt, transpose(K), ZP, -1.0, 1.0) end # Pstartt = Pstar-Pstar*Z'*Kinf-Pinf*Z'*Kstar %(5.14) DK(2012) -function get_updated_Pstartt!(Pstartt::AbstractMatrix{T}, Pstar::AbstractMatrix{T}, ZPstar::AbstractMatrix{T}, - Kinf::AbstractMatrix{T}, ZPinf::AbstractMatrix{T}, - Kstar::AbstractMatrix{T}, vPinftt::AbstractMatrix{T}, PTmp::AbstractMatrix{T}) where T <: AbstractFloat +function get_updated_Pstartt!(Pstartt::AbstractArray{T}, Pstar::AbstractArray{T}, ZPstar::AbstractArray{T}, + Kinf::AbstractArray{T}, ZPinf::AbstractArray{T}, + Kstar::AbstractArray{T}, vPinftt::AbstractArray{T}, PTmp::AbstractArray{T}) where T <: AbstractFloat copy!(Pstartt, Pstar) mul!(Pstartt, transpose(ZPstar), Kinf, -1.0, 1.0) mul!(Pstartt, transpose(ZPinf), Kstar, -1.0, 1.0) end # v = y - Z*a -- basic -function get_v!(v::AbstractVector{T}, y::AbstractVecOrMat{Union{T, Missing}}, z::AbstractVecOrMat{T}, a::AbstractVector{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer} +function get_v!(v::AbstractVector{T}, y::AbstractArray{Union{T, Missing}}, z::AbstractVecOrMat{T}, a::AbstractVector{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer} copyto!(v, 1, y, iy, ny) gemv!('N', -1.0, z, a, 1.0, v) end # v = y - Z*a -- basic -- univariate -function get_v!(Y::AbstractVecOrMat{Union{T, Missing}}, Z::AbstractVecOrMat{T}, a::AbstractVector{T}, i::U) where {T <: AbstractFloat, U <: Integer} +function get_v!(Y::AbstractVector{Union{T, Missing}}, Z::AbstractVecOrMat{T}, a::AbstractVector{T}, i::U) where {T <: AbstractFloat, U <: Integer} v = Y[i] @inbounds @simd for j = 1:length(a) v -= Z[i, j]*a[j] @@ -315,37 +315,37 @@ function get_v!(Y::AbstractVecOrMat{Union{T, Missing}}, Z::AbstractVecOrMat{T}, end # v = y - a[z] -- Z selection matrix -function get_v!(v::AbstractVector{T}, y::AbstractVecOrMat{Union{T, Missing}}, z::AbstractVector{U}, a::AbstractVector{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer} +function get_v!(v::AbstractVector{T}, y::AbstractArray{Union{T, Missing}}, z::AbstractVector{U}, a::AbstractVector{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer} copyto!(v, 1, y, iy, ny) az = view(a,z) v .= v .- az end # v = y - a[z] -- Z selection matrix -- univariate -function get_v!(y::AbstractVecOrMat{Union{T, Missing}}, z::AbstractVector{U}, a::AbstractVector{T}, i::U) where {T <: AbstractFloat, U <: Integer} +function get_v!(y::AbstractVector{Union{T, Missing}}, z::AbstractVector{U}, a::AbstractVector{T}, i::U) where {T <: AbstractFloat, U <: Integer} return y[i] - a[z[i]] end # v = y - Z*a -- missing observations -function get_v!(v::AbstractVector{T}, y::AbstractVecOrMat{Union{T, Missing}}, z::AbstractMatrix{T}, a::AbstractVector{T}, t::U, pattern::Vector{U}) where {T <: AbstractFloat, U <: Integer} +function get_v!(v::AbstractVector{T}, y::AbstractMatrix{Union{T, Missing}}, z::AbstractArray{T}, a::AbstractVector{T}, t::U, pattern::Vector{U}) where {T <: AbstractFloat, U <: Integer} v .= view(y, pattern, t) gemv!('N', -1.0, z, a, 1.0, v) end # v = y - a[z] -- Z selection matrix and missing variables -function get_v!(v::AbstractVector{T}, y::AbstractVecOrMat{Union{T, Missing}}, z::AbstractVector{U}, a::AbstractVector{T}, t::U, pattern::Vector{Int64}) where {T <: AbstractFloat, U <: Integer} +function get_v!(v::AbstractVector{T}, y::AbstractMatrix{Union{T, Missing}}, z::AbstractVector{U}, a::AbstractVector{T}, t::U, pattern::Vector{Int64}) where {T <: AbstractFloat, U <: Integer} v .= view(y, pattern, t) .- view(a, z) end # v = y - c - Z*a -- basic -function get_v!(v::AbstractArray{T}, y::AbstractVecOrMat{Union{T, Missing}}, c::AbstractVector{T}, z::AbstractArray{T}, a::AbstractArray{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer} +function get_v!(v::AbstractArray{T}, y::AbstractArray{Union{T, Missing}}, c::AbstractVector{T}, z::AbstractArray{T}, a::AbstractArray{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer} copyto!(v, 1, y, iy, ny) v .-= c gemm!('N', 'N', -1.0, z, a, 1.0, v) end # v = y - c - Z*a -- basic -- univariate -function get_v!(Y::AbstractVecOrMat{Union{T, Missing}}, c::AbstractVector{T}, Z::AbstractVecOrMat{T}, a::AbstractVector{T}, i::U) where {T <: AbstractFloat, U <: Integer} +function get_v!(Y::AbstractVector{Union{T, Missing}}, c::AbstractVector{T}, Z::AbstractVecOrMat{T}, a::AbstractVector{T}, i::U) where {T <: AbstractFloat, U <: Integer} v = Y[i] - c[i] @inbounds @simd for j = 1:length(a) v -= Z[i, j]*a[j] @@ -354,25 +354,25 @@ function get_v!(Y::AbstractVecOrMat{Union{T, Missing}}, c::AbstractVector{T}, Z: end # v = y - c - a[z] -- Z selection matrix -function get_v!(v::AbstractArray{T}, y::AbstractVecOrMat{Union{T, Missing}}, c::AbstractArray{T}, z::AbstractVector{U}, a::AbstractArray{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer} +function get_v!(v::AbstractVector{T}, y::AbstractArray{Union{T, Missing}}, c::AbstractVector{T}, z::AbstractVector{U}, a::AbstractArray{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer} copyto!(v, 1, y, iy, ny) az = view(a,z) v .-= c .+ az end # v = y - c - a[z] -- Z selection matrix -- univariate -function get_v!(y::AbstractVecOrMat{Union{T, Missing}}, c::AbstractVector{T}, z::AbstractVector{U}, a::AbstractVector{T}, i::U) where {T <: AbstractFloat, U <: Integer} +function get_v!(y::AbstractVector{Union{T, Missing}}, c::AbstractVector{T}, z::AbstractVector{U}, a::AbstractVector{T}, i::U) where {T <: AbstractFloat, U <: Integer} return y[i] - c[i] - a[z[i]] end # v = y - c - Z*a -- missing observations -function get_v!(v::AbstractArray{T}, y::AbstractVecOrMat{Union{T, Missing}}, c::AbstractArray{T}, z::AbstractArray{T}, a::AbstractArray{T}, t::U, pattern::Vector{U}) where {T <: AbstractFloat, U <: Integer} +function get_v!(v::AbstractVector{T}, y::AbstractMatrix{Union{T, Missing}}, c::AbstractArray{T}, z::AbstractArray{T}, a::AbstractArray{T}, t::U, pattern::Vector{U}) where {T <: AbstractFloat, U <: Integer} v .= view(y, pattern, t) .- view(c, pattern) gemm!('N', 'N', -1.0, z, a, 1.0, v) end # v = y - c - a[z] -- Z selection matrix and missing variables -function get_v!(v::AbstractArray{T}, y::AbstractVecOrMat{Union{T, Missing}}, c::AbstractArray{T}, z::AbstractVector{U}, a::AbstractArray{T}, t::U, pattern::Vector{Int64}) where {T <: AbstractFloat, U <: Integer} +function get_v!(v::AbstractVector{T}, y::AbstractMatrix{Union{T, Missing}}, c::AbstractArray{T}, z::AbstractVector{U}, a::AbstractArray{T}, t::U, pattern::Vector{Int64}) where {T <: AbstractFloat, U <: Integer} v .= view(y, pattern, t) .- view(c, pattern) .- view(a, z) end @@ -386,11 +386,11 @@ end # Valpha_t = Pstar_t - Pstar_t*N0_{t-1}*Pstar_t # -(Pinf_t*N1_{t-1}*Pstar_t)' -Pinf_t*N1_{t-1}*Pstar_t - # -Pinf_t*N2_{t-1}*Pinf_t (DK 5.30) -function get_Valpha!(Valpha::AbstractMatrix{T}, - Pstar::AbstractMatrix{T}, - Pinf::AbstractMatrix{T}, N0::AbstractMatrix{T}, - N1::AbstractMatrix{T}, N2::AbstractMatrix{T}, - Tmp::AbstractMatrix{T}) where T <: AbstractFloat +function get_Valpha!(Valpha::AbstractArray{T}, + Pstar::AbstractArray{T}, + Pinf::AbstractArray{T}, N0::AbstractArray{T}, + N1::AbstractArray{T}, N2::AbstractArray{T}, + Tmp::AbstractArray{T}) where T <: AbstractFloat copy!(Valpha, Pstar) mul!(Tmp, Pstar, N0) mul!(Tmp, Pinf, N1, 1.0, 1.0) @@ -415,7 +415,7 @@ function get_Veta!(Veta, Q, R, N, RQ, tmp) mul!(Veta, transpose(RQ), tmp, -1.0, 1.0) end -function get_vZsmall(Zsmall::AbstractMatrix{T}, iZsmall::AbstractVector{U}, Z::AbstractArray{T}, pattern::AbstractVector{U}, n::U, ny::U) where {T <: AbstractFloat, U <: Integer} +function get_vZsmall(Zsmall::AbstractArray{T}, iZsmall::AbstractVector{U}, Z::AbstractArray{T}, pattern::AbstractVector{U}, n::U, ny::U) where {T <: AbstractFloat, U <: Integer} vZsmall = view(Zsmall, 1:n, :) if n == ny copyto!(vZsmall, Z) @@ -425,7 +425,7 @@ function get_vZsmall(Zsmall::AbstractMatrix{T}, iZsmall::AbstractVector{U}, Z::A return vZsmall end -function get_vZsmall(Zsmall::AbstractMatrix{T}, iZsmall::AbstractVector{U}, z::AbstractArray{U}, pattern::AbstractVector{U}, n::U, ny::U) where {T <: AbstractFloat, U <: Integer} +function get_vZsmall(Zsmall::AbstractArray{T}, iZsmall::AbstractVector{U}, z::AbstractArray{U}, pattern::AbstractVector{U}, n::U, ny::U) where {T <: AbstractFloat, U <: Integer} vZsmall = view(iZsmall, 1:n) if n == ny copyto!(vZsmall, z) @@ -435,7 +435,7 @@ function get_vZsmall(Zsmall::AbstractMatrix{T}, iZsmall::AbstractVector{U}, z::A return vZsmall end -function get_vZsmall(Zsmall::AbstractMatrix{T}, iZsmall::AbstractVector{U}, Z::AbstractArray{T}, pattern::AbstractVector{U}, n::U, ny::U, t::U) where {T <: AbstractFloat, U <: Integer} +function get_vZsmall(Zsmall::AbstractArray{T}, iZsmall::AbstractVector{U}, Z::AbstractArray{T}, pattern::AbstractVector{U}, n::U, ny::U, t::U) where {T <: AbstractFloat, U <: Integer} changeZ = ndims(Z) > 2 vZ = changeZ ? view(Z, :, :, t) : view(Z, :, :) vZsmall = view(Zsmall, 1:n, :) @@ -447,7 +447,7 @@ function get_vZsmall(Zsmall::AbstractMatrix{T}, iZsmall::AbstractVector{U}, Z::A return vZsmall end -function get_vZsmall(Zsmall::AbstractMatrix{T}, iZsmall::AbstractVector{U}, z::AbstractArray{U}, pattern::AbstractVector{U}, n::U, ny::U, t::U) where {T <: AbstractFloat, U <: Integer} +function get_vZsmall(Zsmall::AbstractArray{T}, iZsmall::AbstractVector{U}, z::AbstractArray{U}, pattern::AbstractVector{U}, n::U, ny::U, t::U) where {T <: AbstractFloat, U <: Integer} changeZ = ndims(z) > 1 vz = changeZ ? view(z, :, t) : view(z, :,) vZsmall = view(iZsmall, 1:n) @@ -460,7 +460,7 @@ function get_vZsmall(Zsmall::AbstractMatrix{T}, iZsmall::AbstractVector{U}, z::A end # a = T(a + K'*v) -function update_a!(a::AbstractVector{U}, K::AbstractMatrix{U}, v::AbstractVector{U}, a1::Vector{U}, T::AbstractMatrix{U}) where U <: AbstractFloat +function update_a!(a::AbstractVector{U}, K::AbstractArray{U}, v::AbstractVector{U}, a1::Vector{U}, T::AbstractArray{U}) where U <: AbstractFloat copy!(a1, a) gemv!('T', 1.0, K, v, 1.0, a1) gemv!('N', 1.0, T, a1, 0.0, a) @@ -475,7 +475,7 @@ function update_a!(a1::AbstractArray{U}, a::AbstractArray, d::AbstractArray{U}, end # a = d + T*att -function update_a!(a1::AbstractVector{U}, d::AbstractVector{U}, T::AbstractMatrix{U}, att::AbstractVector{U}) where U <: AbstractFloat +function update_a!(a1::AbstractVector{U}, d::AbstractVector{U}, T::AbstractArray{U}, att::AbstractVector{U}) where U <: AbstractFloat copy!(a1, d) mul!(a1, T, att, 1.0, 1.0) end @@ -533,20 +533,20 @@ end # F1 = inv(Finf) # N1_{t-1} = Z'*F1*Z + L0'*N1_t*L0 + L1'*N0_t*L0 -function update_N1!(N1::AbstractMatrix{T}, Z::AbstractMatrix{T}, - iFZ::AbstractMatrix{T}, L0::AbstractMatrix{T}, - N1_1::AbstractMatrix{T}, L1::AbstractMatrix{T}, - N0_1::AbstractMatrix{T}, PTmp::AbstractMatrix{T}) where T <: AbstractFloat +function update_N1!(N1::AbstractArray{T}, Z::AbstractArray{T}, + iFZ::AbstractArray{T}, L0::AbstractArray{T}, + N1_1::AbstractArray{T}, L1::AbstractArray{T}, + N0_1::AbstractArray{T}, PTmp::AbstractArray{T}) where T <: AbstractFloat mul!(N1, transpose(Z), iFZ) mul!(PTmp, transpose(L0), N1_1) mul!(PTmp, transpose(L1), N0_1, 1.0, 1.0) mul!(N1, PTmp, L0, 1.0, 1.0) end -function update_N1!(N1::AbstractMatrix{T}, z::AbstractVector{U}, - iFZ::AbstractMatrix{T}, L0::AbstractMatrix{T}, - N1_1::AbstractMatrix{T}, L1::AbstractMatrix{T}, - N0_1::AbstractMatrix{T}, PTmp::AbstractMatrix{T}) where {T <: AbstractFloat, U <: Integer} +function update_N1!(N1::AbstractArray{T}, z::AbstractVector{U}, + iFZ::AbstractArray{T}, L0::AbstractArray{T}, + N1_1::AbstractArray{T}, L1::AbstractArray{T}, + N0_1::AbstractArray{T}, PTmp::AbstractArray{T}) where {T <: AbstractFloat, U <: Integer} vN1 = view(N1, z, :) vN1 .= iFZ mul!(PTmp, transpose(L0), N1_1) @@ -557,11 +557,11 @@ end # F2 = -inv(Finf)*Fstar*inv(Finv) # N2_{t-1} = Z'*F2*Z + L0'*N2_t*L0 + L0'*N1_t*L1 # + L1'*N1_t*L0 + L1'*N0_t*L1 -function update_N2!(N2::AbstractMatrix{T}, iFZ::AbstractMatrix{T}, - Fstar::AbstractMatrix{T}, L0::AbstractMatrix{T}, - N2_1::AbstractMatrix{T}, N1_1::AbstractMatrix{T}, - L1::AbstractMatrix{T}, N0_1::AbstractMatrix{T}, - Tmp1::AbstractMatrix{T}, Tmp2::AbstractMatrix{T}) where T <: AbstractFloat +function update_N2!(N2::AbstractArray{T}, iFZ::AbstractArray{T}, + Fstar::AbstractArray{T}, L0::AbstractArray{T}, + N2_1::AbstractArray{T}, N1_1::AbstractArray{T}, + L1::AbstractArray{T}, N0_1::AbstractArray{T}, + Tmp1::AbstractArray{T}, Tmp2::AbstractArray{T}) where T <: AbstractFloat mul!(Tmp1, Fstar, iFZ) mul!(N2, transpose(iFZ), Tmp1, -1.0, 0.0) mul!(Tmp2, transpose(L0), N2_1) @@ -587,14 +587,14 @@ function filtered_P!(P1::AbstractArray{U}, P::AbstractArray{U}, K::AbstractArray end # P = T*Ptt*T' + QQ -function update_P!(P::AbstractMatrix{U}, T::AbstractMatrix{U}, Ptt::AbstractMatrix{U}, QQ::AbstractMatrix{U}, Ptmp::AbstractMatrix{U}) where U <: AbstractFloat +function update_P!(P::AbstractArray{U}, T::AbstractArray{U}, Ptt::AbstractArray{U}, QQ::AbstractArray{U}, Ptmp::AbstractArray{U}) where U <: AbstractFloat mul!(Ptmp, Ptt, transpose(T)) copy!(P, QQ) mul!(P, T, Ptmp, 1.0, 1.0) end # Pinf = T*Pinftt*T' -function update_P!(P::AbstractMatrix{U}, T::AbstractMatrix{U}, Ptt::AbstractMatrix{U}, Ptmp::AbstractMatrix{U}) where U <: AbstractFloat +function update_P!(P::AbstractArray{U}, T::AbstractArray{U}, Ptt::AbstractArray{U}, Ptmp::AbstractArray{U}) where U <: AbstractFloat mul!(Ptmp, Ptt, transpose(T)) mul!(P, T, Ptmp) end @@ -632,9 +632,9 @@ function update_r!(r1::AbstractVector{T}, z::AbstractVector{U}, iFv::AbstractVec end # r1_{t-1} = Z_t'*iFinf_t*v_t + L0_t'r1_t + L1_t'*0_t (DK 5.21) -function update_r!(r1::AbstractVector{T}, Z::AbstractMatrix{T}, - iFv::AbstractVector{T}, L0::AbstractMatrix{T}, - r1_1::AbstractVector{T}, L1::AbstractMatrix{T}, +function update_r!(r1::AbstractVector{T}, Z::AbstractArray{T}, + iFv::AbstractVector{T}, L0::AbstractArray{T}, + r1_1::AbstractVector{T}, L1::AbstractArray{T}, r0_1::AbstractVector{T}) where T <: AbstractFloat mul!(r1, transpose(Z), iFv) mul!(r1, transpose(L0), r1_1, 1.0, 1.0) @@ -643,8 +643,8 @@ end # r1_{t-1} = z_t'*iFinf_t*v_t + L0_t'r1_t + L1_t'*r0_t (DK 5.21) function update_r!(r1::AbstractVector{T}, z::AbstractVector{U}, - iFv::AbstractVector{T}, L0::AbstractMatrix{T}, - r1_1::AbstractVector{T}, L1::AbstractMatrix{T}, + iFv::AbstractVector{T}, L0::AbstractArray{T}, + r1_1::AbstractVector{T}, L1::AbstractArray{T}, r0_1::AbstractVector{T}) where {T <: AbstractFloat, U <: Integer} vr1 = view(r1, z) vr1 .= iFv diff --git a/src/kalman_likelihood.jl b/src/kalman_likelihood.jl index e1cda1fb49681f2f3b4c686bd4f99eb085126032..7b3b0ef5aeed2cf05dce8142ff3396148909b93b 100644 --- a/src/kalman_likelihood.jl +++ b/src/kalman_likelihood.jl @@ -1,54 +1,54 @@ struct KalmanLikelihoodWs{T, U} <: KalmanWs{T, U} csmall::Vector{T} - Zsmall::Matrix{T} + Zsmall::Array{T} # necessary for Z selecting vector with missing variables iZsmall::Vector{U} - RQ::Matrix{T} - QQ::Matrix{T} + RQ::Array{T} + QQ::Array{T} v::Vector{T} - F::Matrix{T} - cholF::Matrix{T} - cholH::Matrix{T} - LTcholH::Matrix{T} + F::Array{T} + cholF::Array{T} + cholH::Array{T} + LTcholH::Array{T} iFv::Vector{T} a1::Vector{T} - K::Matrix{T} - ZP::Matrix{T} + K::Array{T} + ZP::Array{T} iFZ::SubArray{T} - PTmp::Matrix{T} - oldP::Matrix{T} + PTmp::Array{T} + oldP::Array{T} lik::Vector{T} cholHset::Bool ystar::Vector{T} - Zstar::Matrix{T} - Hstar::Matrix{T} + Zstar::Array{T} + Hstar::Array{T} tmp_ns::Vector{T} PZi::Vector{T} kalman_tol::T function KalmanLikelihoodWs{T, U}(ny::U, ns::U, np::U, nobs::U) where {T <: AbstractFloat, U <: Integer} csmall = Vector{T}(undef, ny) - Zsmall = Matrix{T}(undef, ny, ns) + Zsmall = Array{T}(undef, ny, ns) iZsmall = Vector{U}(undef, ny) - RQ = Matrix{T}(undef, ns, np) - QQ = Matrix{T}(undef, ns, ns) - F = Matrix{T}(undef, ny, ny) - cholF = Matrix{T}(undef, ny, ny) - cholH = Matrix{T}(undef, ny, ny) - LTcholH = Matrix{T}(undef, ny, ny) + RQ = Array{T}(undef, ns, np) + QQ = Array{T}(undef, ns, ns) + F = Array{T}(undef, ny, ny) + cholF = Array{T}(undef, ny, ny) + cholH = Array{T}(undef, ny, ny) + LTcholH = Array{T}(undef, ny, ny) v = Vector{T}(undef, ny) iFv = Vector{T}(undef, ny) a1 = Vector{T}(undef, ns) - K = Matrix{T}(undef, ny, ns) - PTmp = Matrix{T}(undef, ns, ns) - oldP = Matrix{T}(undef, ns, ns) - ZP = Matrix{T}(undef, ny, ns) + K = Array{T}(undef, ny, ns) + PTmp = Array{T}(undef, ns, ns) + oldP = Array{T}(undef, ns, ns) + ZP = Array{T}(undef, ny, ns) iFZ = view(PTmp,1:ny,:) lik = Vector{T}(undef, nobs) cholHset = false ystar = Vector{T}(undef, ny) - Zstar = Matrix{T}(undef, ny, ns) - Hstar = Matrix{T}(undef, ny, ny) + Zstar = Array{T}(undef, ny, ns) + Hstar = Array{T}(undef, ny, ny) tmp_ns = Vector{T}(undef, ns) PZi = Vector{T}(undef, ns) kalman_tol = 1e-12 @@ -77,7 +77,7 @@ function kalman_likelihood(Y::AbstractArray{Union{U, Missing}}, ny = size(Y,1) nobs = last - start + 1 # QQ = R*Q*R' - get_QQ!(ws.QQ, R, Q, ws.RQ) + @inbounds get_QQ!(ws.QQ, R, Q, ws.RQ) cholHset = false t = start iy = (start - 1)*ny + 1 @@ -109,8 +109,8 @@ function kalman_likelihood(Y::AbstractArray{Union{U, Missing}}, t += 1 end lik_cst = (nobs - presample)*ny*log(2*pi) - vlik = view(ws.lik, start + presample:last) - return -0.5*(lik_cst + sum(vlik)) + @inbounds vlik = view(ws.lik, start + presample:last) + return @inbounds -0.5*(lik_cst + sum(vlik)) end function kalman_likelihood(Y::AbstractArray{Union{U, Missing}}, @@ -129,7 +129,7 @@ function kalman_likelihood(Y::AbstractArray{Union{U, Missing}}, ny = size(Y,1) nobs = last - start + 1 # QQ = R*Q*R' - get_QQ!(ws.QQ, R, Q, ws.RQ) + @inbounds get_QQ!(ws.QQ, R, Q, ws.RQ) l2pi = log(2*pi) t = start iy = (start - 1)*ny + 1 @@ -174,18 +174,18 @@ function kalman_likelihood(Y::AbstractArray{Union{U, Missing}}, end t += 1 end - vlik = view(ws.lik, start + presample:last) - return -0.5*sum(vlik) + @inbounds vlik = view(ws.lik, start + presample:last) + return @inbounds -0.5*sum(vlik) end -function kalman_likelihood_monitored(Y::Matrix{U}, +function kalman_likelihood_monitored(Y::AbstractArray{Union{U, Missing}}, Z::AbstractArray{W}, - H::Matrix{U}, - T::Matrix{U}, - R::Matrix{U}, - Q::Matrix{U}, + H::AbstractArray{U}, + T::AbstractArray{U}, + R::AbstractArray{U}, + Q::AbstractArray{U}, a::Vector{U}, - P::Matrix{U}, + P::AbstractArray{U}, start::V, last::V, presample::V, @@ -194,11 +194,11 @@ function kalman_likelihood_monitored(Y::Matrix{U}, nobs = last - start + 1 ns = size(T,1) # QQ = R*Q*R' - get_QQ!(ws.QQ, R, Q, ws.RQ) + @inbounds get_QQ!(ws.QQ, R, Q, ws.RQ) t = start iy = (start - 1)*ny + 1 steady = false - copy!(ws.oldP, P) + @inbounds copy!(ws.oldP, P) @inbounds while t <= last # v = Y[:,t] - Z*a get_v!(ws.v, Y, Z, a, iy, ny) @@ -231,7 +231,7 @@ function kalman_likelihood_monitored(Y::Matrix{U}, # P = T*(P - K'*Z*P)*T'+ QQ update_P!(P, T, ws.QQ, ws.K, ws.ZP, ws.PTmp) ws.oldP .-= P - if norm(ws.oldP) < ns*eps() + if norm(ws.oldP) < ns*eps()*5e10 steady = true else copy!(ws.oldP, P) @@ -240,8 +240,8 @@ function kalman_likelihood_monitored(Y::Matrix{U}, t += 1 end lik_cst = (nobs - presample)*ny*log(2*pi) - vlik = view(ws.lik, start + presample:last) - return -0.5*(lik_cst + sum(vlik)) + @inbounds vlik = view(ws.lik, start + presample:last) + return @inbounds -0.5*(lik_cst + sum(vlik)) end function kalman_likelihood_monitored(Y::AbstractArray{Union{U, Missing}}, @@ -261,12 +261,12 @@ function kalman_likelihood_monitored(Y::AbstractArray{Union{U, Missing}}, nobs = last - start + 1 ns = size(T,1) # QQ = R*Q*R' - get_QQ!(ws.QQ, R, Q, ws.RQ) + @inbounds get_QQ!(ws.QQ, R, Q, ws.RQ) l2pi = log(2*pi) t = start iy = (start - 1)*ny + 1 steady = false - copy!(ws.oldP, P) + @inbounds copy!(ws.oldP, P) cholHset = false @inbounds while t <= last pattern = data_pattern[t] @@ -310,7 +310,7 @@ function kalman_likelihood_monitored(Y::AbstractArray{Union{U, Missing}}, # P = T*(P - K'*Z*P)*T'+ QQ update_P!(P, T, ws.QQ, vK, vZP, ws.PTmp) ws.oldP .-= P - if norm(ws.oldP) < ns*eps() + if norm(ws.oldP) < ns*eps()*5e5 steady = true else copy!(ws.oldP, P) @@ -318,63 +318,63 @@ function kalman_likelihood_monitored(Y::AbstractArray{Union{U, Missing}}, end t += 1 end - vlik = view(ws.lik, start + presample:last) - return -0.5*sum(vlik) + @inbounds vlik = view(ws.lik, start + presample:last) + return @inbounds -0.5*sum(vlik) end struct FastKalmanLikelihoodWs{T, U} <: KalmanWs{T, U} csmall::Vector{T} - Zsmall::Matrix{T} + Zsmall::Array{T} iZsmall::Vector{U} - QQ::Matrix{T} + QQ::Array{T} v::Vector{T} - F::Matrix{T} - cholF::Matrix{T} + F::Array{T} + cholF::Array{T} iFv::Vector{T} a1::Vector{T} - K::Matrix{T} - RQ::Matrix{T} - ZP::Matrix{T} - M::Matrix{T} - W::Matrix{T} - ZW::Matrix{T} - ZWM::Matrix{T} - iFZWM::Matrix{T} - TW::Matrix{T} - iFZW::Matrix{T} - KtiFZW::Matrix{T} + K::Array{T} + RQ::Array{T} + ZP::Array{T} + M::Array{T} + W::Array{T} + ZW::Array{T} + ZWM::Array{T} + iFZWM::Array{T} + TW::Array{T} + iFZW::Array{T} + KtiFZW::Array{T} ystar::Vector{T} - Zstar::Matrix{T} - Hstar::Matrix{T} + Zstar::Array{T} + Hstar::Array{T} PZi::Vector{T} lik::Vector{T} kalman_tol::T function FastKalmanLikelihoodWs{T, U}(ny::U, ns::U, np::U, nobs::U) where {T <: AbstractFloat, U <: Integer} csmall = Vector{T}(undef, ny) - Zsmall = Matrix{T}(undef, ny, ns) + Zsmall = Array{T}(undef, ny, ns) iZsmall = Vector{U}(undef, ny) - QQ = Matrix{T}(undef, ns, ns) - RQ = Matrix{T}(undef, ns, np) - F = Matrix{T}(undef, ny, ny) - cholF = Matrix{T}(undef, ny, ny) + QQ = Array{T}(undef, ns, ns) + RQ = Array{T}(undef, ns, np) + F = Array{T}(undef, ny, ny) + cholF = Array{T}(undef, ny, ny) v = Vector{T}(undef, ny) iFv = Vector{T}(undef, ny) a1 = Vector{T}(undef, ns) - K = Matrix{T}(undef, ny, ns) - M = Matrix{T}(undef, ny, ny) - W = Matrix{T}(undef, ns, ny) - ZP = Matrix{T}(undef, ny, ns) - ZW = Matrix{T}(undef, ny, ny) - ZWM = Matrix{T}(undef, ny, ny) - iFZWM = Matrix{T}(undef, ny, ny) - TW = Matrix{T}(undef, ns, ny) - iFZW = Matrix{T}(undef, ny, ny) - KtiFZW = Matrix{T}(undef, ns, ny) + K = Array{T}(undef, ny, ns) + M = Array{T}(undef, ny, ny) + W = Array{T}(undef, ns, ny) + ZP = Array{T}(undef, ny, ns) + ZW = Array{T}(undef, ny, ny) + ZWM = Array{T}(undef, ny, ny) + iFZWM = Array{T}(undef, ny, ny) + TW = Array{T}(undef, ns, ny) + iFZW = Array{T}(undef, ny, ny) + KtiFZW = Array{T}(undef, ns, ny) ystar = Vector{T}(undef, ny) - Zstar = Matrix{T}(undef, ny, ns) - Hstar = Matrix{T}(undef, ny, ny) + Zstar = Array{T}(undef, ny, ns) + Hstar = Array{T}(undef, ny, ny) PZi = Vector{T}(undef, ns) lik = Vector{T}(undef, nobs) kalman_tol = 1e-12 @@ -388,14 +388,14 @@ FastKalmanLikelihoodWs(ny, ns, np, nobs) = FastKalmanLikelihoodWs{Float64, Int64 """ K doesn't represent the same matrix as above """ -function fast_kalman_likelihood(Y::Matrix{U}, +function fast_kalman_likelihood(Y::Array{U}, Z::AbstractArray{W}, - H::Matrix{U}, - T::Matrix{U}, - R::Matrix{U}, - Q::Matrix{U}, + H::Array{U}, + T::Array{U}, + R::Array{U}, + Q::Array{U}, a::Vector{U}, - P::Matrix{U}, + P::Array{U}, start::V, last::V, presample::V, @@ -441,14 +441,14 @@ function fast_kalman_likelihood(Y::Matrix{U}, return -0.5*(lik_cst + sum(vlik)) end -function fast_kalman_likelihood(Y::Matrix{U}, +function fast_kalman_likelihood(Y::Array{U}, Z::AbstractArray{W}, - H::Matrix{U}, - T::Matrix{U}, - R::Matrix{U}, - Q::Matrix{U}, + H::Array{U}, + T::Array{U}, + R::Array{U}, + Q::Array{U}, a::Vector{U}, - P::Matrix{U}, + P::Array{U}, start::V, last::V, presample::V, @@ -512,57 +512,57 @@ end struct DiffuseKalmanLikelihoodWs{T, U} <: KalmanWs{T, U} csmall::Vector{T} - Zsmall::Matrix{T} + Zsmall::Array{T} iZsmall::Vector{U} - QQ::Matrix{T} - RQ::Matrix{T} + QQ::Array{T} + RQ::Array{T} v::Vector{T} - F::Matrix{T} - iF::Matrix{T} + F::Array{T} + iF::Array{T} iFv::Vector{T} a1::Vector{T} - cholF::Matrix{T} - ZP::Matrix{T} - Fstar::Matrix{T} - ZPstar::Matrix{T} - K::Matrix{T} - iFZ::Matrix{T} - Kstar::Matrix{T} - PTmp::Matrix{T} + cholF::Array{T} + ZP::Array{T} + Fstar::Array{T} + ZPstar::Array{T} + K::Array{T} + iFZ::Array{T} + Kstar::Array{T} + PTmp::Array{T} uKinf::Vector{T} uKstar::Vector{T} Kinf_Finf::Vector{T} ystar::Vector{T} - Zstar::Matrix{T} - Hstar::Matrix{T} + Zstar::Array{T} + Hstar::Array{T} PZi::Vector{T} lik::Vector{T} kalman_tol::T function DiffuseKalmanLikelihoodWs{T, U}(ny::U, ns::U, np::U, nobs::U) where {T <: AbstractFloat, U <: Integer} csmall = Vector{T}(undef, ny) - Zsmall = Matrix{T}(undef, ny, ns) + Zsmall = Array{T}(undef, ny, ns) iZsmall = Vector{U}(undef, ny) - QQ = Matrix{T}(undef, ns, ns) - RQ = Matrix{T}(undef, ns, np) + QQ = Array{T}(undef, ns, ns) + RQ = Array{T}(undef, ns, np) v = Vector{T}(undef, ny) - F = Matrix{T}(undef, ny, ny) - iF = Matrix{T}(undef, ny,ny ) + F = Array{T}(undef, ny, ny) + iF = Array{T}(undef, ny,ny ) iFv = Vector{T}(undef, ny) a1 = Vector{T}(undef, ns) - cholF = Matrix{T}(undef, ny, ny) - ZP = Matrix{T}(undef, ny, ns) - Fstar = Matrix{T}(undef, ny, ny) - ZPstar = Matrix{T}(undef, ny, ns) - K = Matrix{T}(undef, ny, ns) - iFZ = Matrix{T}(undef, ny, ns) - Kstar = Matrix{T}(undef, ny, ns) - PTmp = Matrix{T}(undef, ns, ns) + cholF = Array{T}(undef, ny, ny) + ZP = Array{T}(undef, ny, ns) + Fstar = Array{T}(undef, ny, ny) + ZPstar = Array{T}(undef, ny, ns) + K = Array{T}(undef, ny, ns) + iFZ = Array{T}(undef, ny, ns) + Kstar = Array{T}(undef, ny, ns) + PTmp = Array{T}(undef, ns, ns) uKinf = Vector{T}(undef, ns) uKstar = Vector{T}(undef, ns) Kinf_Finf = Vector{T}(undef, ns) ystar = Vector{T}(undef, ny) - Zstar = Matrix{T}(undef, ny, ns) - Hstar = Matrix{T}(undef, ny, ny) + Zstar = Array{T}(undef, ny, ns) + Hstar = Array{T}(undef, ny, ny) PZi = Vector{T}(undef, ns) lik = zeros(T, nobs) kalman_tol = 1e-12 @@ -574,14 +574,14 @@ end DiffuseKalmanLikelihoodWs(ny, ns, np, nobs) = DiffuseKalmanLikelihoodWs{Float64, Int64}(ny, ns, np, nobs) -function diffuse_kalman_likelihood_init!(Y::Matrix{U}, +function diffuse_kalman_likelihood_init!(Y::Array{U}, Z::AbstractArray{W}, - H::Matrix{U}, - T::Matrix{U}, - QQ::Matrix{U}, + H::Array{U}, + T::Array{U}, + QQ::Array{U}, a::Vector{U}, - Pinf::Matrix{U}, - Pstar::Matrix{U}, + Pinf::Array{U}, + Pstar::Array{U}, start::V, last::V, tol::U, @@ -636,14 +636,14 @@ function diffuse_kalman_likelihood_init!(Y::Matrix{U}, t end -function diffuse_kalman_likelihood_init!(Y::Matrix{U}, +function diffuse_kalman_likelihood_init!(Y::Array{U}, Z::AbstractArray{W}, - H::Matrix{U}, - T::Matrix{U}, - QQ::Matrix{U}, + H::Array{U}, + T::Array{U}, + QQ::Array{U}, a::Vector{U}, - Pinf::Matrix{U}, - Pstar::Matrix{U}, + Pinf::Array{U}, + Pstar::Array{U}, start::V, last::V, tol::U, @@ -714,15 +714,15 @@ function diffuse_kalman_likelihood_init!(Y::Matrix{U}, t end -function diffuse_kalman_likelihood(Y::Matrix{U}, +function diffuse_kalman_likelihood(Y::Array{U}, Z::AbstractArray{W}, - H::Matrix{U}, - T::Matrix{U}, - R::Matrix{U}, - Q::Matrix{U}, + H::Array{U}, + T::Array{U}, + R::Array{U}, + Q::Array{U}, a::Vector{U}, - Pinf::Matrix{U}, - Pstar::Matrix{U}, + Pinf::Array{U}, + Pstar::Array{U}, start::V, last::V, presample::V, @@ -740,15 +740,15 @@ function diffuse_kalman_likelihood(Y::Matrix{U}, return -0.5*(lik_cst + sum(vlik)) end -function diffuse_kalman_likelihood(Y::Matrix{U}, +function diffuse_kalman_likelihood(Y::Array{U}, Z::AbstractArray{W}, - H::Matrix{U}, - T::Matrix{U}, - R::Matrix{U}, - Q::Matrix{U}, + H::Array{U}, + T::Array{U}, + R::Array{U}, + Q::Array{U}, a::Vector{U}, - Pinf::Matrix{U}, - Pstar::Matrix{U}, + Pinf::Array{U}, + Pstar::Array{U}, start::V, last::V, presample::V,