Broken wiki example

[willtebbutt]

The following example from the wiki

mf(vi, sampler; x=[1.5, 2.0]) = begin
  s = Turing.assume(sampler,
                    InverseGamma(2, 3),
                    Turing.VarName(vi, [:c_s, :s], ""),
                    vi)
  m = Turing.assume(sampler,
                    Normal(0,sqrt(s)),
                    Turing.VarName(vi, [:c_m, :m], ""),
                    vi)
  for i = 1:2
    Turing.observe(sampler,
                   Normal(m, sqrt(s)),
                   x[i],
                   vi)
  end
  vi
end
mf(Turing.VarInfo(), nothing)

throws the following error:

ERROR: MethodError: no method matching sqrt(::Tuple{Float64,Float64})
Closest candidates are:
  sqrt(::Float64) at math.jl:425
  sqrt(::Float32) at math.jl:426
  sqrt(::Float16) at math.jl:950
  ...
Stacktrace:
 [1] #mf#1(::Array{Float64,1}, ::Function, ::Turing.VarReplay.VarInfo, ::Void) at ./REPL[3]:6
 [2] mf(::Turing.VarReplay.VarInfo, ::Void) at ./REPL[3]:1

The culprit appears to be that s is a tuple rather than a Float64. What is the appropriate work around? This example worked a couple of weeks ago; what changes to the internals have occurred to cause this problem?

[xukai92]

In the last PR I changed the API of assume, observe a bit. Basically now they also return log-probability: I evaluate the log-probability inside the function and return. Coming with this change, I set a local variable _lp inside the model function.

The reason to do this is that ReverseDiff.jl somehow is very slow with the type is unsafe and using this local variable to accumulate probability somehow makes Turing.jl work fast with ReverseDiff.jl. Though I didn't figure out why this really happened, I guess it's something related with type stability (see related issues #400, #331 (comment))

So now this program should be changed to below

mf(vi, sampler; x=[1.5, 2.0]) = begin
  _lp = zero(Real)

  s, __lp = Turing.assume(sampler,
                    InverseGamma(2, 3),
                    Turing.VarName(vi, [:c_s, :s], ""),
                    vi)
  _lp += __lp
  m, __lp = Turing.assume(sampler,
                    Normal(0,sqrt(s)),
                    Turing.VarName(vi, [:c_m, :m], ""),
                    vi)
  _lp += __lp
  for i = 1:2
    _lp += Turing.observe(sampler,
                   Normal(m, sqrt(s)),
                   x[i],
                   vi)
  end
  vi.logp = _lp
  vi
end
mf(Turing.VarInfo(), nothing)

Tested in commit 43017ad as the current master seems broken?

[willtebbutt]

Makes sense, thanks for the quick response.