Confusing error when sampling from empty model

[penelopeysm]

julia> using Turing

julia> @model function f(x)
           x ~ Normal()
       end
f (generic function with 2 methods)

julia> sample(f(1.0), NUTS(), 100)
Sampling 100%|██████████████████████████████████████████| Time: 0:00:00
ERROR: ArgumentError: reducing over an empty collection is not allowed; consider supplying `init` to the reducer
[...]

I think check_model() should probably error if the model is empty.

[mhauru]

Or maybe VarInfo should be able to handle this without erroring, return a chain with one dimension being zero? I don't see a reason why we should ban likelihood-only models, even if they are a bit silly. Might come around as a special case in some system of models someone is developing or somesuch.

[penelopeysm]

DynamicPPL is actually quite fine with it, it's getparams that errors so in that case this would be a Turing bug. Happy to move it, I opened it here because if it was check_model it would be here.

[penelopeysm]

We might run into this though

julia> m = Array{Float64}(undef, 0, 0, 0)
0×0×0 Array{Float64, 3}

julia> Chains(m, [])
ERROR: StackOverflowError:
Stacktrace:
     [1] Chains(::Array{Float64, 3}, ::Vector{Any}; kwargs::@Kwargs{})
       @ MCMCChains ~/.julia/packages/MCMCChains/RIFg4/src/chains.jl:24
     [2] Chains(::Array{Float64, 3}, ::Vector{Any})
       @ MCMCChains ~/.julia/packages/MCMCChains/RIFg4/src/chains.jl:18
--- the above 2 lines are repeated 26655 more times ---
 [53313] Chains(::Array{Float64, 3}, ::Vector{Any}; kwargs::@Kwargs{})
       @ MCMCChains ~/.julia/packages/MCMCChains/RIFg4/src/chains.jl:24

[penelopeysm]

I thought this would be a quick fix, then found that it's actually quite difficult. The offending code is here:

DynamicPPL.jl/src/varinfo.jl

Lines 854 to 856 in 3acd0a2

	function getindex_internal(vi::NTVarInfo, ::Colon)
	return reduce(vcat, map(Base.Fix2(getindex_internal, Colon()), vi.metadata))
	end

(there is a similar bit of code in Turing's src/mcmc/Inference.jl, with a similar problem to what's described here)

The existing code is type-stable:

using DynamicPPL, Distributions, Test

@model e(x) = x ~ Normal()
@model f() = x ~ Normal()
@model g() = y ~ Poisson()
@model function h()
    x ~ Normal()
    y ~ Poisson()
end

evi = VarInfo(e(1.0)); fvi = VarInfo(f()); gvi = VarInfo(g()); hvi = VarInfo(h());

# evi[:] errors -- this is the issue we want to fix
@inferred fvi[:] # Vector{Float64}
@inferred gvi[:] # Vector{Int64}
@inferred hvi[:] # Vector{Float64}, because Int64 gets promoted to Float64

This breaks type stability:

- return reduce(vcat, map(Base.Fix2(getindex_internal, Colon()), vi.metadata))
+ return reduce(vcat, map(Base.Fix2(getindex_internal, Colon()), vi.metadata), init=[])

This is type-stable but gvi[:] will now return Vector{Float64}:

- return reduce(vcat, map(Base.Fix2(getindex_internal, Colon()), vi.metadata))
+ return reduce(vcat, map(Base.Fix2(getindex_internal, Colon()), vi.metadata), init=Float64[])

Seems like an option is to declare a new method for an empty varinfo.