`UndefVarError` in simple Poisson-Gamma model

[slwu89]

Hi RxInfer team. I'm reporting a UndefVarError bug here that I found as a result of the MWE I posted in this question on the Julia Discourse (https://discourse.julialang.org/t/rxinfer-jl-beginner-question-re-undefreferror/126660), where I was encouraged to report this to the team.

I am on MacOS, Julia v1.11.3, and software versions:

  [31c24e10] Distributions v0.25.117
  [86711068] RxInfer v4.2.0

Here is the MWE:

using Distributions

"""
Follows 3.2.1 in the paper. Simulate data from a clinical trial enrollment process.
Note β is a rate parameter, but Distributions uses scale; we take 1/β in the function.
Note that the paper considers the "time step" to be a day.
"""
function pg_homogeneous(α, β, u, T)
    N = length(α)
    @assert N == length(β) == length(u)
    λ = rand.(Gamma.(α, 1 ./ β)) # shape, rate parameterization
    outmat = zeros(Int, N, T)
    for i in 1:N
        rate = [k < u[i] ? 0 : λ[i] for k in 1:T]
        outmat[i, :] = rand.(Poisson.(rate))
    end
    return cumsum(outmat, dims=2) # return cumulative number of pts
end

# example with time-homogeneous rates
α = 1/(1.2^2)
β = α/0.02
N = 200
T = 400
u = rand(1:29*4, N)
recruitment_target = 1000

out_homogeneous = pg_homogeneous(fill(α, N), fill(β, N), u, T)

final_date = findfirst(sum(out_homogeneous, dims=1) .>= recruitment_target)[2]
out_homogeneous = out_homogeneous[:, 1:final_date]

using RxInfer

# The inference model for the homogeneous PG model.
#     This assumes that Gamma(α,β) is the common distribution of the 
#     rate for all centers (the simulation function in the paper generalizes that somewhat
#     but all the math in the paper makes this simplifying assumption).
@model function pg_homogeneous_rxinfer(data, T, N, u)
    # these are quite terrible priors for the Gamma(α,β)
    # parameters, but the actual conjugate prior is quite a strange non-standard
    # distribution (see Gamma entry in https://en.wikipedia.org/wiki/Conjugate_prior#When_likelihood_function_is_a_continuous_distribution)
    # so let's use these for testing
    α ~ Uniform(0, 10)
    β ~ Uniform(1, 150)
    
    # the main stochastic model
    for i in 1:N
        λ[i] ~ Gamma(shape=α, scale=β) # shape, scale parameterization
        for k in 1:T
            if k ≥ u[i] # recruitment has started    
                data[i, k] ~ Poisson(λ[i]) # Poisson draw at each day
            end    
        end
    end
end

infer_homogeneous = infer(
    model = pg_homogeneous_rxinfer(T=final_date, N=N, u=u),
    data = (data = out_homogeneous, )
)

And the stack trace:

ERROR: UndefRefError: access to undefined reference
Stacktrace:
  [1] getindex
    @ ./essentials.jl:917 [inlined]
  [2] recursive_getindex
    @ ~/.julia/packages/GraphPPL/xPNyo/src/resizable_array.jl:144 [inlined]
  [3] recursive_getindex
    @ ~/.julia/packages/GraphPPL/xPNyo/src/resizable_array.jl:152 [inlined]
  [4] getindex
    @ ~/.julia/packages/GraphPPL/xPNyo/src/resizable_array.jl:136 [inlined]
  [5] iterate(array::GraphPPL.ResizableArray{GraphPPL.NodeLabel, Vector{Vector{GraphPPL.NodeLabel}}, 2})
    @ GraphPPL ~/.julia/packages/GraphPPL/xPNyo/src/resizable_array.jl:180
  [6] iterate
    @ ./generator.jl:45 [inlined]
  [7] _collect(c::GraphPPL.ResizableArray{…}, itr::Base.Generator{…}, ::Base.EltypeUnknown, isz::Base.HasShape{…})
    @ Base ./array.jl:811
  [8] collect_similar(cont::GraphPPL.ResizableArray{…}, itr::Base.Generator{…})
    @ Base ./array.jl:720
  [9] map(f::Function, A::GraphPPL.ResizableArray{GraphPPL.NodeLabel, Vector{Vector{GraphPPL.NodeLabel}}, 2})
    @ Base ./abstractarray.jl:3371
 [10] getvarref(model::GraphPPL.Model{…}, container::GraphPPL.ResizableArray{…})
    @ RxInfer ~/.julia/packages/RxInfer/2Syro/src/model/plugins/reactivemp_inference.jl:265
 [11] (::RxInfer.var"#73#74"{GraphPPL.Model{…}})(v::GraphPPL.ResizableArray{GraphPPL.NodeLabel, Vector{…}, 2})
    @ RxInfer ~/.julia/packages/RxInfer/2Syro/src/model/plugins/reactivemp_inference.jl:261
 [12] map!(f::RxInfer.var"#73#74"{…}, out::Dictionaries.UnorderedDictionary{…}, d::Dictionaries.UnorderedDictionary{…})
    @ Dictionaries ~/.julia/packages/Dictionaries/tq7Jt/src/map.jl:57
 [13] map(f::Function, d::Dictionaries.UnorderedDictionary{Symbol, Any})
    @ Dictionaries ~/.julia/packages/Dictionaries/tq7Jt/src/map.jl:92
 [14] map(f::Function, vardict::GraphPPL.VarDict{Any})
    @ GraphPPL ~/.julia/packages/GraphPPL/xPNyo/src/graph_engine.jl:615
 [15] getvardict
    @ ~/.julia/packages/RxInfer/2Syro/src/model/plugins/reactivemp_inference.jl:261 [inlined]
 [16] getvardict
    @ ~/.julia/packages/RxInfer/2Syro/src/model/model.jl:40 [inlined]
 [17] batch_inference(; model::GraphPPL.ModelGenerator{…}, data::@NamedTuple{…}, initialization::Nothing, constraints::Nothing, meta::Nothing, options::Nothing, returnvars::Nothing, predictvars::Nothing, iterations::Nothing, free_energy::Bool, free_energy_diagnostics::Tuple{…}, allow_node_contraction::Bool, showprogress::Bool, callbacks::Nothing, addons::Nothing, postprocess::DefaultPostprocess, warn::Bool, catch_exception::Bool)
    @ RxInfer ~/.julia/packages/RxInfer/2Syro/src/inference/batch.jl:205
 [18] batch_inference
    @ ~/.julia/packages/RxInfer/2Syro/src/inference/batch.jl:96 [inlined]
 [19] #288
    @ ~/.julia/packages/RxInfer/2Syro/src/inference/inference.jl:532 [inlined]
 [20] with_session(f::RxInfer.var"#288#290"{…}, session::RxInfer.Session, label::Symbol)
    @ RxInfer ~/.julia/packages/RxInfer/2Syro/src/session.jl:253
 [21] #infer#287
    @ ~/.julia/packages/RxInfer/2Syro/src/inference/inference.jl:502 [inlined]
 [22] top-level scope
    @ ~/Desktop/misc/tmp.jl:58
Some type information was truncated. Use `show(err)` to see complete types.

[nsajko]

FWIW this reproduces on all versions of Julia I tried, so from v1.10 to v1.13 (nightly). I also tried disabling compilation with --compile=min, disabling inlining with --inline=no, and other compiler options, nothing helped, suggesting this is not an issue with Julia.

[bvdmitri]

Hey thanks for reporting the issue, we will try to figure out whats happening here

[ma-seefelder]

Hi, I get a highly similar error message with my hierarchical Bayesian model and can't determine what causes it. Hence, I added debugging messages at different parts of the model to verify that my rather complicated indexing works and the indexing is fine. Furthermore, the last info message after computing the likelihood is displayed.

@model function HierarchicalBayesianModel(samples, controls, μ, σ, a, b, nparameters, n_protocols, experiment_positions, protocol_positions, matched_positions, n_experiments, parameter_lookup) 
    #######################################
    # prior definitions                 ###
    #######################################
    # define priors
    local μ_control #   # Individual means of the control group
    local μ_sample #    # Individual means of the sample group
    local σ_control #   # Individual variance of the control group
    local σ_sample #    # Individual variance deviations of the sample group

    # 1st level priors
    σ_control[1] ~  Gamma(shape = a, scale = b) 
    σ_sample[1]  ~  Gamma(shape = a, scale = b)
    μ_control[1] ~  Normal(mean = μ, precision = 1.0 / σ)
    μ_sample[1]  ~  Normal(mean = μ, precision = 1.0 / σ)

    # 2nd and 3rd level priors
    experiment_idx = 1

    for idx ∈ 2:nparameters
        σ_control[idx]  ~  Gamma(shape = a, scale = b)  
        σ_sample[idx]   ~  Gamma(shape = a, scale = b)

        # 2nd level priors
        if idx ∈ protocol_positions
            # hyperprior for all experiments of protocol "protocol" and individual experiments
            μ_control[idx]  ~   Normal(mean = μ_control[1], precision = σ_control[1])
            μ_sample[idx]   ~   Normal(mean = μ_sample[1], precision = σ_sample[1]) 
        elseif idx ∈ experiment_positions
            # retrieve protocol position
            protocol_position = matched_positions[experiment_idx]
            experiment_idx += 1
            # define low-level priors
            μ_control[idx]  ~   Normal(mean = μ_control[protocol_position], precision = σ_control[protocol_position])
            μ_sample[idx]   ~   Normal(mean = μ_sample[protocol_position], precision = σ_sample[protocol_position])
        else
            throw(BoundsError("Position $idx is not defined"))
        end
    end

    #######################################
    # likelihood definitions
    #######################################
    for protocol ∈ 1:n_protocols, experiment ∈ 1:n_experiments[protocol]
        pos_parameter_vector = parameter_lookup[protocol, experiment]
        for idx in size(samples, 3) 
            @info "$(protocol), $(experiment), $(idx)"
            samples[protocol, experiment, idx]  ~ Normal(mean = μ_sample[pos_parameter_vector],  precision = σ_sample[pos_parameter_vector])
            controls[protocol, experiment, idx] ~ Normal(mean = μ_control[pos_parameter_vector], precision = σ_control[pos_parameter_vector])
        end
    end
@info "Likelihood computation completed"
end

Stack trace

ERROR: UndefRefError: access to undefined reference
Stacktrace:
  [1] getindex
    @ .\essentials.jl:917 [inlined]
  [2] recursive_getindex
    @ C:\Users\manue\.julia\packages\GraphPPL\xPNyo\src\resizable_array.jl:144 [inlined]
  [3] recursive_getindex (repeats 2 times)
    @ C:\Users\manue\.julia\packages\GraphPPL\xPNyo\src\resizable_array.jl:152 [inlined]
  [4] getindex
    @ C:\Users\manue\.julia\packages\GraphPPL\xPNyo\src\resizable_array.jl:136 [inlined]
  [5] iterate(array::GraphPPL.ResizableArray{GraphPPL.NodeLabel, Vector{Vector{Vector{GraphPPL.NodeLabel}}}, 3})
    @ GraphPPL C:\Users\manue\.julia\packages\GraphPPL\xPNyo\src\resizable_array.jl:180
  [6] iterate
    @ .\generator.jl:45 [inlined]
  [7] _collect(c::GraphPPL.ResizableArray{…}, itr::Base.Generator{…}, ::Base.EltypeUnknown, isz::Base.HasShape{…})
    @ Base .\array.jl:811
  [8] collect_similar(cont::GraphPPL.ResizableArray{…}, itr::Base.Generator{…})
    @ Base .\array.jl:720
  [9] map(f::Function, A::GraphPPL.ResizableArray{GraphPPL.NodeLabel, Vector{Vector{Vector{GraphPPL.NodeLabel}}}, 3})
    @ Base .\abstractarray.jl:3371
 [10] getvarref(model::GraphPPL.Model{…}, container::GraphPPL.ResizableArray{…})
    @ RxInfer C:\Users\manue\.julia\packages\RxInfer\t5qay\src\model\plugins\reactivemp_inference.jl:277
 [11] (::RxInfer.var"#73#74"{GraphPPL.Model{…}})(v::GraphPPL.ResizableArray{GraphPPL.NodeLabel, Vector{…}, 3})
    @ RxInfer C:\Users\manue\.julia\packages\RxInfer\t5qay\src\model\plugins\reactivemp_inference.jl:273
 [12] map!(f::RxInfer.var"#73#74"{…}, out::Dictionaries.UnorderedDictionary{…}, d::Dictionaries.UnorderedDictionary{…})
    @ Dictionaries C:\Users\manue\.julia\packages\Dictionaries\tq7Jt\src\map.jl:57
 [13] map(f::Function, d::Dictionaries.UnorderedDictionary{Symbol, Any})
    @ Dictionaries C:\Users\manue\.julia\packages\Dictionaries\tq7Jt\src\map.jl:92
 [14] map(f::Function, vardict::GraphPPL.VarDict{Any})
    @ GraphPPL C:\Users\manue\.julia\packages\GraphPPL\xPNyo\src\graph_engine.jl:615
 [15] getvardict
    @ C:\Users\manue\.julia\packages\RxInfer\t5qay\src\model\plugins\reactivemp_inference.jl:273 [inlined]
 [16] getvardict
    @ C:\Users\manue\.julia\packages\RxInfer\t5qay\src\model\model.jl:40 [inlined]
 [17] batch_inference(; model::GraphPPL.ModelGenerator{…}, data::@NamedTuple{…}, initialization::RxInfer.InitSpecification, constraints::GraphPPL.Constraints{…}, meta::Nothing, options::Nothing, returnvars::RxInfer.KeepLast, predictvars::Nothing, iterations::Int64, free_energy::Bool, free_energy_diagnostics::Tuple{…}, allow_node_contraction::Bool, showprogress::Bool, callbacks::Nothing, addons::Nothing, postprocess::RxInfer.DefaultPostprocess, warn::Bool, catch_exception::Bool)
    @ RxInfer C:\Users\manue\.julia\packages\RxInfer\t5qay\src\inference\batch.jl:209
 [18] batch_inference
    @ C:\Users\manue\.julia\packages\RxInfer\t5qay\src\inference\batch.jl:96 [inlined]
 [19] #314
    @ C:\Users\manue\.julia\packages\RxInfer\t5qay\src\inference\inference.jl:532 [inlined]
 [20] with_session(f::RxInfer.var"#314#316"{…}, session::RxInfer.Session, label::Symbol)
    @ RxInfer C:\Users\manue\.julia\packages\RxInfer\t5qay\src\session.jl:253
 [21] #infer#313
    @ C:\Users\manue\.julia\packages\RxInfer\t5qay\src\inference\inference.jl:502 [inlined]
 [22] infer
    @ C:\Users\manue\.julia\packages\RxInfer\t5qay\src\inference\inference.jl:455 [inlined]
 [23] HBM(data::Main.BayesInteractomics.InteractionData{…}, idx::Int64; μ_0::Float64, σ_0::Float64, a_0::Float64, b_0::Float64)
    @ Main.BayesInteractomics c:\Users\manue\Documents\GitHub\BayesInteractomics\src\model_fitting.jl:234
 [24] HBM
    @ c:\Users\manue\Documents\GitHub\BayesInteractomics\src\model_fitting.jl:194 [inlined]
 [25] main(data::Main.BayesInteractomics.InteractionData{…}, idx::Int64, referenceID::Int64; ndraws::Int64, μ_0::Nothing, σ_0::Nothing, a_0::Nothing, b_0::Nothing, α::Float64, csv_file::String, plotHBMdists::Bool, plotlog2fc::Bool, plotregr::Bool, plotbayesrange::Bool, writecsv::Bool)
    @ Main.BayesInteractomics c:\Users\manue\Documents\GitHub\BayesInteractomics\src\model_fitting.jl:581
 [26] main(data::Main.BayesInteractomics.InteractionData{Int64}, idx::Int64, referenceID::Int64)
    @ Main.BayesInteractomics c:\Users\manue\Documents\GitHub\BayesInteractomics\src\model_fitting.jl:552
 [27] top-level scope
    @ c:\Users\manue\Documents\GitHub\BayesInteractomics\src\meta_analysis.jl:47
Some type information was truncated. Use `show(err)` to see complete types.

[wouterwln]

Hey @ma-seefelder , thanks for trying out RxInfer! The problem you guys are having is a bug inRxInfer. What is going on is the following: when you write your model with the @model macro, we construct a Factor Graph under the hood. We dynamically scale the containers in which you can store random variables, which means you can use indexed statements with the ~ operator. The bug comes across when we collect the posterior distributions. RxInfer assumes you fully initialize all variables in a container, which means that your model fully materializes and inference fully runs (which is why your debug statements all hit), but in the collection of posteriors we run into an error. You can avoid this by initializing all variables in a container (not skipping indices) for now. We will work on the issue whenever we have time, but the RxInfer team is pretty swamped with other tasks for the foreseeable future. If you have any other questions (maybe we can work out a hot fix) let me know!

[ma-seefelder]

Hey @wouterwln, thank you for your fast response and the explanations. I will try to initialise the containers beforehand. I already used RxInfer.jl several times after I followed the presentation at the JuliaCon and really enjoy using it, mainly because of its speed. Thanks for the good work to the whole RxInfer team.

[wouterwln]

Posting this for completeness:

Hi! We are moving to a cleaner Epic -> Feature -> Task issue hierarchy to better organize our backlog. This issue is currently either underspecified or not tagged appropriately.

To keep this issue open, please do the following within the next 7 days (by 25-11-2025):

1. Update/Replace: Ensure the description is clear and actionable.
2. Tag Correctly:For Tasks/Features, add the correct label (e.g., feature, task) AND include a link to the Parent Epic or Feature it belongs to.
3. For Bugs, add the Bug label. (Bugs do not require a parent link.)

Issues not updated, linked, or tagged correctly by the deadline will be closed and purged.

Thank you for helping us clean up and organize our backlog!

For this specific issue, I will file the bug report since I've figured out where in RxInfer the bug occurs, feel free to disregard this, but to give you a heads up why this issue will be closed without a fix implemented