mschauer / kalman.jl Goto Github PK

Hi,

I am writing to kindly request an extra feature to make H an array of matrix-free operators, where H[i] is the observation operator of the i-th time step and its usage could be matrix-free (with only access to its forward and adjoint evaluations). Maybe \Phi (evolution) as well.

Please let me know if this is meaningful to the Kalman community (it is quite important for my personal projects involving large-scale physics-driven inverse problems) and I am willing to contribute.

Cheers,

Ziyi (Francis) Yin

Contributions, comments and feature requests are most welcome.

I will settle for an interface which combines an Observation model and a Evolution model with a Filter method, so e.g.

"""
    LinearObservation(H, R)

Observe `y = Hx + v` where ``v ~ N(0, R)``
"""
struct LinearObservation{TH, TR} <: AbstractObservation
    H::TH # dxd
    R::TR # d
end

can be combined with

struct JosephForm <: FilterMethod
end

The new interface is mostly implemented, I'll do a pull request which relocates the old "LinearHomogSystem" in a less prominent place.

Any comments?

When I try the example, at the loop step

for i in 1:N
       global p
       p = Φ*p ⊕ Gaussian(zero(x0), Q)
       p, yres, _ = Kalman.correct(Kalman.JosephForm(), p, (Gaussian(ys[i], R), H))
       push!(ps, p)
end

I get an error and Julia crashes

Illegal inttoptr
      %77 = ptrtoint double addrspace(13)* %75 to i64
Illegal inttoptr
      %scevgep4546 = ptrtoint double addrspace(13)* %scevgep45 to i64
Illegal inttoptr
      %138 = inttoptr i64 %137 to i8 addrspace(13)*

signal (6): Aborted
in expression starting at none:0
gsignal at /usr/bin/../lib/x86_64-linux-gnu/libc.so.6 (unknown line)
abort at /usr/bin/../lib/x86_64-linux-gnu/libc.so.6 (unknown line)
unknown function (ip: 0x7fb280c8dd04)
_ZN4llvm13FPPassManager13runOnFunctionERNS_8FunctionE at /usr/bin/../lib/x86_64-linux-gnu/libLLVM-8.so.1 (unknown line)
_ZN4llvm13FPPassManager11runOnModuleERNS_6ModuleE at /usr/bin/../lib/x86_64-linux-gnu/libLLVM-8.so.1 (unknown line)
_ZN4llvm6legacy15PassManagerImpl3runERNS_6ModuleE at /usr/bin/../lib/x86_64-linux-gnu/libLLVM-8.so.1 (unknown line)
unknown function (ip: 0x7fb280d77ac1)
unknown function (ip: 0x7fb280d7a2d8)
unknown function (ip: 0x7fb280d7a8cd)
unknown function (ip: 0x7fb280cb5b4a)
unknown function (ip: 0x7fb280ce6ed3)
unknown function (ip: 0x7fb280d0c84b)
jl_apply_generic at /usr/bin/../lib/x86_64-linux-gnu/libjulia.so.1 (unknown line)
top-level scope at ./none:4
unknown function (ip: 0x7fb280d3a030)
jl_toplevel_eval_in at /usr/bin/../lib/x86_64-linux-gnu/libjulia.so.1 (unknown line)
eval at ./boot.jl:331 [inlined]
repleval at /home/liam/.julia/packages/Atom/BPtI0/src/repl.jl:196
#242 at /home/liam/.julia/packages/Atom/BPtI0/src/repl.jl:226
with_logstate at ./logging.jl:398
with_logger at ./logging.jl:505 [inlined]
evalrepl at /home/liam/.julia/packages/Atom/BPtI0/src/repl.jl:214
unknown function (ip: 0x7fb280d1f75b)
unknown function (ip: 0x7fb280d1f389)
unknown function (ip: 0x7fb280d1f8f0)
unknown function (ip: 0x7fb280d209c8)
unknown function (ip: 0x7fb280d21616)
unknown function (ip: 0x7fb280d39e08)
jl_toplevel_eval_in at /usr/bin/../lib/x86_64-linux-gnu/libjulia.so.1 (unknown line)
eval at ./boot.jl:331
eval_user_input at /build/julia-98cBbp/julia-1.4.1+dfsg/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
macro expansion at /build/julia-98cBbp/julia-1.4.1+dfsg/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:118 [inlined]
#26 at ./task.jl:358
unknown function (ip: 0x7fb280d24cbb)
unknown function (ip: (nil))
Allocations: 76408531 (Pool: 76388747; Big: 19784); GC: 71

Julia has exited.
Press Enter to start a new session.

I know very little about Julia, is this a problem with this package, or with my installation of Julia?

Thanks for your great work!

Is there a way to use Kalman.jl from inside ModelinToolkit.jl? I was thinking to implement it manually in MTK, but maybe that's unnessecary.

Thanks!

After running the Kalman filter on some data, I'd like to sample from the posterior, but I get

> rand(estimates[end])
ArgumentError: matrix is not symmetric/Hermitian. This error can be avoided by calling chol(Hermitian(A)) which will ignore either the upper or lower triangle of the matrix.

Stacktrace:
 [1] chol(::Array{Float64,2}) at ./linalg/cholesky.jl:185
 [2] rand(::GaussianDistributions.Gaussian{Array{Float64,1},Array{Float64,2}}) at /home/cst-jean/.julia/v0.6/GaussianDistributions/src/GaussianDistributions.jl:68
 [3] include_string(::String, ::String) at ./loading.jl:522

I admit that I don't understand these issues as well as I should. Should I just follow their advice, and use P.μ + chol(Hermitian(P.Σ))'*randn(T, length(P.μ))? It seems to work.

Hi Moritz!
I'm not sure, but I think evolve and therefore predict! ignore the temporal arguments s and t. So if we have a system where the speed x is in km per hour but the states are measured once per minute then the estimated progress is exaggerated by 60.

Ref. https://github.com/JuliaFolds/Transducers.jl

It would speed the startup...

w ~ N(0, Q) and v ~ N(0, R) are not displayed correctly in the documentation. The tilde is missing for both. Checked this for Safari and Chrome on MacOS. In Juno, it's displayed correctly.

By Julia's ! convention, I expected predict! to modify... something, but it doesn't seem to. Shall I make a PR to take out the !?

function predict!(s, x, P, t, M::LinearHomogSystem)
    x = evolve(s, x, P, t, M)
    Ppred = M.Phi*P*M.Phi' + M.Q
    x, Ppred, M.Phi
end

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Thank you for the track function, but I feel that it doesn't belong in here. This function is useful only for a small set of circumstances and just adds clutter to most of the other users. In short, if Kalman.jl is to be the backbone for all Kalman filtering then there shouldn't be any ultra specialized functions here, just general ones. I guess.