Turing.jl: Freeze when calling sample

[baggepinnen]

On julia v1.0.1, turing v0.5 and master:

When running the following model, julia freezes upon call to sample

using Turing
weight = randn(10)
N = length(weight)

@model weightfilter(y) = begin
    N = length(y)
    shock = Vector{Real}(undef,N)
    noise = Vector{Real}(undef,N)
    weight = Vector{Real}(undef,N)
    for i = 1:N
        shock[i] ~ Poisson(1) 
        noise[i] ~ Normal(0, 0.2) 
    end
    weight[1] ~ Normal(75, 4)
    for n in 2:N
        weight[n] ~ Normal(weight[n-1], 1.)
        y[n] ~ weight[n] + shock[n] + noise[n]
    end
end;

iterations = 20
chain = sample(weightfilter(weight), PG(10,iterations));

The output is

[ Info:  Assume - `shock` is a parameter
[ Info:  Assume - `noise` is a parameter
[ Info:  Assume - `weight` is a parameter (ignoring `weight` found in global scope)
[ Info:  Observe - `y` is an observation

and then nothing ever happens, zero CPU usage. If I press Ctrl-C I get the error

^CERROR: InterruptException:
Stacktrace:
 [1] try_yieldto(::typeof(identity), ::Base.RefValue{Task}) at ./event.jl:196
 [2] yieldto at ./event.jl:182 [inlined]
 [3] schedule_and_wait(::Task, ::Nothing) at ./event.jl:141
 [4] schedule_and_wait at ./event.jl:139 [inlined]
 [5] consume(::Task) at /local/home/fredrikb/.julia/packages/Libtask/vTI2s/src/taskcopy.jl:116
 [6] consume at /local/home/fredrikb/.julia/dev/Turing/src/core/trace.jl:37 [inlined]
 [7] consume(::Turing.ParticleContainer{Turing.Trace}) at /local/home/fredrikb/.julia/dev/Turing/src/core/container.jl:86
 [8] step(::Function, ::Turing.Sampler{PG{Any,typeof(Turing.resampleSystematic)}}, ::Turing.VarReplay.VarInfo) at /local/home/fredrikb/.julia/dev/Turing/src/samplers/pgibbs.jl:73
 [9] macro expansion at ./util.jl:213 [inlined]
 [10] #sample#108(::Bool, ::Nothing, ::Int64, ::Function, ::Function, ::PG{Any,typeof(Turing.resampleSystematic)}) at /local/home/fredrikb/.julia/dev/Turing/src/samplers/pgibbs.jl:114
 [11] sample(::Function, ::PG{Any,typeof(Turing.resampleSystematic)}) at /local/home/fredrikb/.julia/dev/Turing/src/samplers/pgibbs.jl:91
 [12] top-level scope at none:0

[trappmartin]

Thanks for reporting. Could you try to rename the global variable weight so that there is no name clash. I suspect this is causing the problem.

[baggepinnen]

I renamed and restarted julia but the behavior is unfortunately the same.

[mohamed82008]

I can reproduce the problem.

[mohamed82008]

Culprit: the consume(t.task::Task) function on https://github.com/TuringLang/Turing.jl/blob/88a3643335ff88f763b3568553b6777316a50e67/src/core/trace.jl#L37.

[trappmartin]

Ok, I think the problem is the following line:

y[n] ~ weight[n] + shock[n] + noise[n]

as we assume that the right-hand side of ~ is always a distribution or a vector of distributions. However, in this case, it is the sum of RV values. Replacing this with:

y[n] = weight[n] + shock[n] + noise[n]

should solve the issue. I'll have a try myself.

[trappmartin]

@mohamed82008 I don't think we have a problem with the consume function. This has been thoroughly tested.

[mohamed82008]

Yup seems like I was tracking down a side effect. Your fix solves the problem on my machine.

[trappmartin]

Looks good on my machine too. @baggepinnen could you please check if this works for you and confirm?

[mohamed82008]

This issue shows the need for a useful error message in cases like this. I imagine mixing ~ and = maybe a common bug in user code so a useful error message will be good to have.

[baggepinnen]

I can confirm that the switch to = no longer freezes julia. The code does however not seem to do what I would like it to, which is most likely due to me not having a clue what I'm doing (trying this sort of tools for the first time).

Thanks!

[trappmartin]

Closing issue in favour of issue #575.

[baggepinnen]

I created a topic on discourse explaining what I was trying to do, since now I think the issue is on my side with my understanding of Turing and MCMC methods.
https://discourse.julialang.org/t/observation-is-sum-of-random-variables-in-turing/15963
Your expert help would be appreciated if you have time, thanks :)

[mohamed82008]

I am no expert, but I think the example on Discourse should work, worth reopening this issue with the new example.

using Turing
weight = 75 .+ cumsum(0.2randn(20)) # Generate testdata
N = length(weight)

@model weightfilter(y) = begin
    N = length(y)
    shock = Vector{Real}(undef,N-1) # noise vector
    state = Vector{Real}(undef,N) # The true weight state
    for i = 1:N-1
        shock[i] ~ Exponential(.5) # Ate before etc.
        # noise[i] ~ Normal(0, 0.2) # Hydration level etc.
    end
    state[1] ~ Normal(77, 4) # Initial state
    for n in 2:N
        state[n] ~ Normal(state[n-1], .2) # The state drifts normally
        stateshock = state[n] + shock[n-1] # The measurement is corrupted with exponential noise
        y[n] ~ Normal(stateshock, 0.2) # Measurement is also corrupted with gaussian noise
    end
end;

iterations = 500
chain = sample(weightfilter(weight), SMC(iterations)); 

This freezes on my machine.

[trappmartin]

Correct, this should work but the one above should not.

[yebai]

Hi @mohamed82008 @baggepinnen, for SMC / PG to work correctly, you need to replace Vector with TArray in Turing. The following example is working on my laptop with Turing#master. Please see https://github.com/TuringLang/Libtask.jl/blob/master/README.md for a reference. But Perhaps we should document in turing doc as well.

using Turing
weight = 75 .+ cumsum(0.2randn(20)) # Generate testdata
N = length(weight)

@model weightfilter(y) = begin
    N = length(y)
    shock = tzeros(N-1) # noise vector
    state =  tzeros(N) # The true weight state
    for i = 1:N-1
        shock[i] ~ Exponential(.5) # Ate before etc.
        # noise[i] ~ Normal(0, 0.2) # Hydration level etc.
    end
    state[1] ~ Normal(77, 4) # Initial state
    for n in 2:N
        state[n] ~ Normal(state[n-1], .2) # The state drifts normally
        stateshock = state[n] + shock[n-1] # The measurement is corrupted with exponential noise
        y[n] ~ Normal(stateshock, 0.2) # Measurement is also corrupted with gaussian noise
    end
end;

iterations = 500
chain = sample(weightfilter(weight), SMC(iterations)); 

[trappmartin]

Thanks @yebai, I completely missed that.

[yebai]

PS. the following code leads to more accurate SMC/PG inference results as it delays the sampling of shock. This delayed sampling trick is only not needed if an HMC is used to sample shock however.

using Turing
weight = 75 .+ cumsum(0.2randn(20)) # Generate testdata
N = length(weight)

@model weightfilter(y) = begin
    N = length(y)
    shock = Vector{Real}(undef,N-1) # noise vector
    state = Vector{Real}(undef,N) # The true weight state

   state[1] ~ Normal(77, 4) # Initial state
    for n in 2:N
       shock[i-1] ~ Exponential(.5) # Ate before etc.
        # noise[i-1] ~ Normal(0, 0.2) # Hydration level etc.
        state[n] ~ Normal(state[n-1], .2) # The state drifts normally
        stateshock = state[n] + shock[n-1] # The measurement is corrupted with exponential noise
        y[n] ~ Normal(stateshock, 0.2) # Measurement is also corrupted with gaussian noise
    end
end;

iterations = 500
chain = sample(weightfilter(weight), SMC(iterations)); 

[mohamed82008]

The example above works on my Linux machine and freezes on my Windows machine.

[yebai]

The example above works on my Linux machine and freezes on my Windows machine.

It should work on Windows, Linux and OSX. Perhaps try ] build Libtask to rebuild Libtask.

[mohamed82008]

Perhaps try ] build Libtask to rebuild Libtask.

No difference.

[cpfiffer]

FYI: TArray and tzeros are now in the public API page in the documentation. It's probably worth adding a bigger note somewhere (probably in the guide) going over the specific cases where they're needed.

[yebai]

The example above works on my Linux machine and freezes on my Windows machine.

I don't have a Windows computer to try this out; perhaps we should resolve the Windows issue in a separate issue.

[mohamed82008]

perhaps we should resolve the Windows issue in a separate issue.

@baggepinnen, you do not appear to be on a Windows machine so does the fix above resolve your issue?

[hhhh5]

I have the same issue on Linux and Turing v0.5.1, the program freezes after calling sample with zero CPU usage.

using Random, Turing

Random.seed!(345346);

OBS = Array{Float64}(undef,100,10);
OBS[:] = rand(Uniform(),length(OBS));

@model contamination(SNPs) = begin

    N,J = size(SNPs)
    K = 3

    gt = TArray(Int,(N,J))
    for n in 1:N, j in 1:J
        gt[n,j] ~ Categorical(ones(Real,K)/K)
    end

    culprit = TArray(Real,J)
    culprit ~ [Uniform()]

    μ = TArray(Real,(K,J))
    μ[1,:] ~ [Normal(0.1 ,0.03)]
    μ[2,:] ~ [Normal(0.5 ,0.03)]
    μ[3,:] ~ [Normal(0.95,0.03)]

    σ = TArray(Real,K)
    σ[1] ~ Normal(0.025,0.01)
    σ[2] ~ Normal(0.05 ,0.01)
    σ[3] ~ Normal(0.025,0.01)

    mix = TArray(Real,N)
    mix ~ [Beta(80,5)]

    for n in 1:N, j in 1:J
        μ′ = mix[n] * μ[gt[n,j],j] + (1-mix[n]) * culprit[j]
        σ′ = mix[n] * σ[gt[n,j]]
        SNPs[n,j] ~ Normal(μ′,σ′)
    end

    return mix
end;

tmp = sample(contamination(OBS), SMC(10));

The only output is

[ Info:  Assume - `gt` is a parameter
[ Info:  Assume - `culprit` is a parameter
[ Info:  Assume - `μ` is a parameter
[ Info:  Assume - `σ` is a parameter
[ Info:  Assume - `mix` is a parameter
[ Info:  Observe - `SNPs` is an observation

And the error message after aborting the computation

ERROR: InterruptException:
Stacktrace:
 [1] try_yieldto(::typeof(identity), ::Base.RefValue{Task}) at ./event.jl:196
 [2] yieldto at ./event.jl:182 [inlined]
 [3] schedule_and_wait(::Task, ::Nothing) at ./event.jl:141
 [4] schedule_and_wait at ./event.jl:139 [inlined]
 [5] consume(::Task) at /home/heissj01/.julia/packages/Libtask/rHi8r/src/taskcopy.jl:116
 [6] consume at /home/heissj01/.julia/packages/Turing/YnwiD/src/core/trace.jl:37 [inlined]
 [7] consume(::Turing.ParticleContainer{Turing.Trace}) at /home/heissj01/.julia/packages/Turing/YnwiD/src/core/container.jl:86
 [8] sample(::Function, ::SMC{Union{},typeof(Turing.resampleSystematic)}) at /home/heissj01/.julia/packages/Turing/YnwiD/src/samplers/smc.jl:80
 [9] top-level scope at none:0

I already tried many different combinations of Tarray/Matrix/Vector, with and without vectorization, but so far no success.

[cpfiffer]

@hhhh5 I believe this problem has to do with how Distributions handles multivariate draws. We've had some discussions here about how to treat this (see #476) , but I don't believe we have a good cure-all yet.

There is a quick-fix, which is to iterate through each element of a TArray and assign a single-variate distribution to it. The version you posted was hanging on culprit, μ, and mix. I've modified those distribution assignments so each element is distributed independently, so please try this instead:

using Random, Turing

Random.seed!(345346);

OBS = Array{Float64}(undef,100,10);
OBS[:] = rand(Uniform(),length(OBS));

@model contamination(SNPs) = begin
    N,J = size(SNPs)
    K = 3

    gt = TArray(Int,(N,J))
    for n in 1:N, j in 1:J
        gt[n,j] ~ Categorical(ones(Real,K)/K)
    end

    culprit = TArray(Real,J)
    for i in eachindex(culprit)
        culprit[i] ~ Uniform()
    end

    μ = TArray(Real,(K,J))
    μ[1,:] ~ MvNormal(fill(0.1, J), fill(0.03, J))
    μ[2,:] ~ MvNormal(fill(0.5, J), fill(0.03, J))
    μ[3,:] ~ MvNormal(fill(0.95, J), fill(0.03, J))

    σ = TArray(Real,K)
    σ[1] ~ Normal(0.025,0.01)
    σ[2] ~ Normal(0.05 ,0.01)
    σ[3] ~ Normal(0.025,0.01)

    mix = TArray(Real,N)
    for i in eachindex(mix)
        mix[i] ~ Beta(80,5)
    end

    for n in 1:N, j in 1:J
        μ′ = mix[n] * μ[gt[n,j],j] + (1-mix[n]) * culprit[j]
        σ′ = mix[n] * σ[gt[n,j]]
        SNPs[n,j] ~ Normal(μ′,σ′)
    end

    return mix
end;

tmp = sample(contamination(OBS), SMC(10));