Handling missing

[trappmartin]

The following MWE breaks in the current version of Turing.

using LinearAlgebra
using Turing

# Linear Model
Δt = 0.1
A = [1.0 Δt; 0.0 1.0]
Q = [0.25 * Δt^4 0.5 * Δt^3; 0.5 * Δt^3 Δt^2] + 1e-6 * I

H = [1.0 0.0; 0.0 1.0;]
R = 0.1 * I

@model function hmm(states, ::Type{T} = Float64) where T
    # Initialize arrays for the time-series
    observations = Matrix{T}(undef, 2, length(states))

    # Prior for first step
    states[1] ~ MvNormal([0.0, 0.5], sqrt(0.001))
    observations[:, 1] ~ MvNormal(states[1], R/Δt)

    for t in 2:length(states)
        states[t] ~ MvNormal(A * states[t-1], Δt*Q)
        observations[:, t] ~ MvNormal(states[t], R/Δt)
    end

    return states
end

# Sample a trial from the prior to be used as the ground truth
hmm_prior = hmm(Union{Missing,Vector{Float64}}[missing for _ in 1:100])
states = hmm_prior()

# Condition on the sampled trial 
hmm_conditioned = hmm(states)

# Perform inference (not working)
chain = sample(hmm_conditioned, SMC(), 1000)

With the stack trace:

ERROR: CTaskException:
Malformed dims
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] TArray{Union{Missing, Array{Float64,1}},1}(::Array{Union{Missing, Array{Float64,1}},1}) at /Users/martin/.julia/packages/Libtask/Zo6uM/src/tarray.jl:36

Any idea why this happens?

When replacing states = hmm_prior() with states = Array{Array{Float64}}(hmm_prior()) everything seems fine.

[cpfiffer]

This is more of a Libtask issue than anything else, I think. Your fix is a good shortcut for the moment, but Libtask should probably be more robust to arrays with union types (Union{Missing, Array{Float64}}, in this case). It's also an issue of states not being tightly typed when it's returned by hmm_prior, which you could probably fix by adding to the return statement:

@model function hmm(states, ::Type{T} = Float64) where T
    # Initialize arrays for the time-series
    observations = Matrix{T}(undef, 2, length(states))

    # Prior for first step
    states[1] ~ MvNormal([0.0, 0.5], sqrt(0.001))
    observations[:, 1] ~ MvNormal(states[1], R/Δt)

    for t in 2:length(states)
        states[t] ~ MvNormal(A * states[t-1], Δt*Q)
        observations[:, t] ~ MvNormal(states[t], R/Δt)
    end

    return map(identity, states)
end

Granted, this is super user-unfriendly, but I can't quite see a robust fix that wouldn't take a ton of code in Libtask to accomplish.

[devmotion]

Actually, Libtask handles Union{Missing,Vector{Float64}} element types just fine, you can create such arrays without any issues.

However, there are two problems here:

• Mainly, the TArray constructors are incorrect, and hence a call of TArray{T,1}(whatever) dispatches to the definition of TArray{T,1}(UndefInitializer) (note the missing :: in the definition!) which fails since the size of the first dimension is not specified
• It is not possible to convert arrays to TArrays by just calling TArray(myarray), instead one has to use the convert method. IMO https://github.com/TuringLang/DynamicPPL.jl/blob/5376534dfd5058bb7a448170f961a9982ac4f317/src/compiler.jl#L372 should be changed to convert.

[devmotion]

Indeed, I just checked that with these two changes the example above works.

[devmotion]

The example above should work with the latest releases of Libtask and DynamicPPL.

[trappmartin]

Awesome thanks. I’ll give it a try and close the issue if it works.