Calling mcmc() from within a function

[arzwa]

Hi,

I'm very interested in using Mamba for a current project. I have enjoyed the tutorial and some examples so far. However I have noticed some rather unexpected errors I believe when calling mcmc() from within another function (also briefly mentioned in this closed issue https://github.com/brian-j-smith/Mamba.jl/issues/105#issuecomment-384496370). For example, calling the blr() function below gives the error ERROR: MethodError: no method matching (::getfield(Main, Symbol("##127#128")))(::Model) The applicable method may be too new: running in world age 25145, while current world is 25151.

Example function:

function blr()
    # Model specification
    model = Model(
        y  = Stochastic(1, (μ, σ²) -> MvNormal(μ, √σ²), false),
        μ  = Logical(1, (X, β) -> X * β, false),
        β  = Stochastic(1, () -> MvNormal(2, √1000)),
        σ² = Stochastic(() -> InverseGamma(0.001, 0.001))
    )

    # Sampling scheme
    scheme1 = [NUTS(:β), Slice(:σ², 3.0)]
    setsamplers!(model, scheme1)

    # Data
    line = Dict{Symbol, Any}(:x => [1, 2, 3, 4, 5], :y => [1, 3, 3, 3, 5])
    line[:X] = [ones(5) line[:x]]

    # Initial values
    inits = [Dict{Symbol, Any}(:y => line[:y], :β => rand(Normal(0, 1), 2),
        :σ² => rand(Gamma(1, 1))) for i in 1:3]

    # MCMC
    sim1 = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)
    return sim1
end 

The following code runs however without any issues:

function blr()
    # Model specification
    model = Model(
        y  = Stochastic(1, (μ, σ²) -> MvNormal(μ, √σ²), false),
        μ  = Logical(1, (X, β) -> X * β, false),
        β  = Stochastic(1, () -> MvNormal(2, √1000)),
        σ² = Stochastic(() -> InverseGamma(0.001, 0.001))
    )

    # Sampling scheme
    scheme1 = [NUTS(:β), Slice(:σ², 3.0)]
    setsamplers!(model, scheme1)

    # Data
    line = Dict{Symbol, Any}(:x => [1, 2, 3, 4, 5], :y => [1, 3, 3, 3, 5])
    line[:X] = [ones(5) line[:x]]

    # Initial values
    inits = [Dict{Symbol, Any}(:y => line[:y], :β => rand(Normal(0, 1), 2),
        :σ² => rand(Gamma(1, 1))) for i in 1:3]

    # MCMC
    return model, line, inits
end 

model, line, inits = blr()
mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)

So am I missing some details here or is this a bug?

Thanks in advance, Arthur

[t-student]

I think this seems to be a bug. I have recreated with the following:

julia> versioninfo()
Julia Version 1.0.1
Commit 0d713926f8 (2018-09-29 19:05 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-8550U CPU @ 1.80GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = "C:\Users\mark\AppData\Local\atom\app-1.32.2\atom.exe" -a
  JULIA_NUM_THREADS = 4

installed package status:

(v1.0) pkg> status
    Status `C:\Users\mark\.julia\environments\v1.0\Project.toml`
  [324073ee] Astro v0.0.0 #master (https://github.com/cormullion/Astro.jl)
  [c52e3926] Atom v0.7.10
  [336ed68f] CSV v0.4.2
  [5d742f6a] CSVFiles v0.10.0
  [159f3aea] Cairo v0.5.6
  [a93c6f00] DataFrames v0.14.1
  [31c24e10] Distributions v0.16.4
  [38e38edf] GLM v1.0.1
  [c91e804a] Gadfly v1.0.0
  [09f84164] HypothesisTests v0.8.0
  [7073ff75] IJulia v1.14.1
  [e5e0dc1b] Juno v0.5.3
  [b964fa9f] LaTeXStrings v1.0.3
  [5424a776] Mamba v0.12.0
  [86f7a689] NamedArrays v0.9.1
  [47be7bcc] ORCA v0.2.0
  [f0f68f2c] PlotlyJS v0.12.0
  [91a5bcdd] Plots v0.21.0
  [438e738f] PyCall v1.18.5
  [d330b81b] PyPlot v2.6.3
  [612083be] Queryverse v0.1.0
  [6f49c342] RCall v0.13.0
  [ce6b1742] RDatasets v0.6.1
  [79098fc4] Rmath v0.5.0
  [60ddc479] StatPlots v0.8.2
  [fce5fe82] Turing v0.5.1
  [c17dfb99] WinRPM v0.4.2
  [8ba89e20] Distributed
  [37e2e46d] LinearAlgebra

Example obtained from docs (https://mambajl.readthedocs.io/en/latest/examples/seeds.html), works when run outside a function with a minor modification to fully qualify the use of I (identity matrix) i.e. changed I to LinearAlgebra.I:

using Distributed
@everywhere using Mamba
import LinearAlgebra

function testymcmc()
  ## Data
  seeds = Dict{Symbol, Any}(
    :r => [10, 23, 23, 26, 17, 5, 53, 55, 32, 46, 10, 8, 10, 8, 23, 0, 3, 22, 15,
           32, 3],
    :n => [39, 62, 81, 51, 39, 6, 74, 72, 51, 79, 13, 16, 30, 28, 45, 4, 12, 41,
           30, 51, 7],
    :x1 => [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    :x2 => [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1]
  )
  seeds[:N] = length(seeds[:r])

  ## Model Specification
  model = Model(

    r = Stochastic(1,
      (alpha0, alpha1, x1, alpha2, x2, alpha12, b, n, N) ->
        UnivariateDistribution[(
          p = invlogit(alpha0 + alpha1 * x1[i] + alpha2 * x2[i] +
                       alpha12 * x1[i] * x2[i] + b[i]);
          Binomial(n[i], p)) for i in 1:N
        ],
      false
    ),

    b = Stochastic(1,
      s2 -> Normal(0, sqrt(s2)),
      false
    ),

    alpha0 = Stochastic(
      () -> Normal(0, 1000)
    ),

    alpha1 = Stochastic(
      () -> Normal(0, 1000)
    ),

    alpha2 = Stochastic(
      () -> Normal(0, 1000)
    ),

    alpha12 = Stochastic(
      () -> Normal(0, 1000)
    ),

    s2 = Stochastic(
      () -> InverseGamma(0.001, 0.001)
    )

  )

  ## Initial Values
  inits = [
    Dict(:r => seeds[:r], :alpha0 => 0, :alpha1 => 0, :alpha2 => 0,
         :alpha12 => 0, :s2 => 0.01, :b => zeros(seeds[:N])),
    Dict(:r => seeds[:r], :alpha0 => 0, :alpha1 => 0, :alpha2 => 0,
         :alpha12 => 0, :s2 => 1, :b => zeros(seeds[:N]))
  ]

  ## Sampling Scheme
  scheme = [AMM([:alpha0, :alpha1, :alpha2, :alpha12], 0.01 * Matrix{Float64}(LinearAlgebra.I, 4, 4)),
            AMWG(:b, 0.01),
            AMWG(:s2, 0.1)]
  setsamplers!(model, scheme)

  ## MCMC Simulations
  sim = mcmc(model, seeds, inits, 12500, burnin=2500, thin=2, chains=2)
end

out = testymcmc()

The error and call stack:

MethodError: no method matching (::getfield(Main, Symbol("##123#124")))(::Model)
The applicable method may be too new: running in world age 25222, while current world is 25232.
Closest candidates are:
  #123(::Model) at none:0 (method too new to be called from this world context.)
setinits!(::ScalarStochastic, ::Model, ::Int64) at dependent.jl:163
setinits!(::Model, ::Dict{Symbol,Any}) at initialization.jl:11
setinits!(::Model, ::Array{Dict{Symbol,Any},1}) at initialization.jl:24
#mcmc#23(::Int64, ::Int64, ::Int64, ::Bool, ::Function, ::Model, ::Dict{Symbol,Any}, ::Array{Dict{Symbol,Any},1}, ::Int64) at mcmc.jl:30
(::getfield(Mamba, Symbol("#kw##mcmc")))(::NamedTuple{(:burnin, :thin, :chains),Tuple{Int64,Int64,Int64}}, ::typeof(mcmc), ::Model, ::Dict{Symbol,Any}, ::Array{Dict{Symbol,Any},1}, ::Int64) at none:0
testymcmc() at julia_tmp4.jl:76
top-level scope at none:0
include_string(::Module, ::String, ::String) at loading.jl:1005
include_string(::Module, ::String, ::String, ::Int64) at eval.jl:30
(::getfield(Atom, Symbol("##114#119")){String,Int64,String})() at eval.jl:94
withpath(::getfield(Atom, Symbol("##114#119")){String,Int64,String}, ::String) at utils.jl:30
withpath at eval.jl:46 [inlined]
#113 at eval.jl:93 [inlined]
with_logstate(::getfield(Atom, Symbol("##113#118")){String,Int64,String}, ::Base.CoreLogging.LogState) at logging.jl:397
with_logger at logging.jl:493 [inlined]
#112 at eval.jl:92 [inlined]
hideprompt(::getfield(Atom, Symbol("##112#117")){String,Int64,String}) at repl.jl:85
macro expansion at eval.jl:91 [inlined]
macro expansion at dynamic.jl:24 [inlined]
(::getfield(Atom, Symbol("##111#116")))(::Dict{String,Any}) at eval.jl:86
handlemsg(::Dict{String,Any}, ::Dict{String,Any}) at comm.jl:164
(::getfield(Atom, Symbol("##19#21")){Array{Any,1}})() at task.jl:259

Please can someone take a look?

Also see: #105

[t-student]

Further to above the problems seems to relate solely to calling the Mamba.mcmc function. I can encapsulate everything else within a function then return the elements in a Dict that can be used to call mcmc outside the function.

[bdeonovic]

It is a known problem with no current solution. Mamba needs to restructure how it handles and passes around anonymous functions.

[t-student]

Thanks Benjamin. I will take a look to see what is involved.

[bdeonovic]

Probably a big task. Other people in other packages have had similar issues: https://github.com/JuliaLang/julia/issues/21356

[reactive-relations]

Just a tip for people who have this problem and need a quick fix:

As a way to get around the problem, while the bug is not fixed, put the model definition outside the function and reference it. Also in the model definition make sure that all variables with input data are named, so that you can change the input values inside the function.

For me this is a much more flexible way to still be able to structure my program, as the model definition itself doesn't change, just the input data.