Possibly can improve constraint set error

[ofSingularMind]

What we get now is

Factorization constraint set at node UniSGP_AMF_29 is not a valid constraint set. Please check your model definition and constraint specification. (Constraint set: BitSetTuples.BoundedBitSetTuple(Bool[1 0 0 0 0; 0 1 1 1 0; 0 1 1 0 0; 0 1 0 1 0; 0 0 0 0 1]))

Maybe we can include a helper string to help users read this constraint set string. For example (acknowledging that I myself don't know how it is read, so just guessing) :

Factorization constraint set at node UniSGP_AMF_29 is not a valid constraint set. Please check your model definition and constraint specification. (Constraint set: BitSetTuples.BoundedBitSetTuple(Bool[1 0 0 0 0; 0 1 1 1 0; 0 1 1 0 0; 0 1 0 1 0; 0 0 0 0 1])).

Hint: The constraint set can be read as follows... There should not be a... (one in a certain position in each group?)... etc.

:)

[wouterwln]

The person picking this up: Ask me about it and I can point you in the right direction :)

[abpolym]

@ofSingularMind Could you create a MWE (minimum working example) which throws this error?
For anyone who wants to solve this issue (maybe this weeks development session?) here is the function that throws this error:

https://github.com/ReactiveBayes/GraphPPL.jl/blob/main/src/plugins/variational_constraints/variational_constraints_engine.jl#L622-L638

[ofSingularMind]

@abpolym Possibly not a "minimum" (not)working example, but the following code modified from RxInfer Examples, Kalman Filtering does produce the error.

using RxInfer, Distributions, StableRNGs, Plots

function generate_data(f, n; seed = 123, x_i_min = -20.0, w_i_min = 20.0, noise = 20.0, real_x_τ = 0.1, real_w_τ = 1.0)

    rng = StableRNG(seed)

    real_x = Vector{Float64}(undef, n)
    real_w = Vector{Float64}(undef, n)
    real_y = Vector{Float64}(undef, n)

    for i in 1:n
        real_x[i] = rand(rng, Normal(x_i_min, sqrt(1.0 / real_x_τ)))
        real_w[i] = rand(rng, Normal(w_i_min, sqrt(1.0 / real_w_τ)))
        real_y[i] = rand(rng, Normal(f(real_x[i], real_w[i]), sqrt(noise)))

        x_i_min = real_x[i]
        w_i_min = real_w[i]
    end
    
    return real_x, real_w, real_y
end

n = 250
real_x, real_w, real_y = generate_data(+, n);

@model function identification_problem(f, y, m_x_0, τ_x_0, a_x, b_x, m_w_0, τ_w_0, a_w, b_w, a_y, b_y)
    
    x0 ~ Normal(mean = m_x_0, precision = τ_x_0)
    τ_x ~ Gamma(shape = a_x, rate = b_x)
    w0 ~ Normal(mean = m_w_0, precision = τ_w_0)
    τ_w ~ Gamma(shape = a_w, rate = b_w)
    τ_y ~ Gamma(shape = a_y, rate = b_y)
    
    x_i_min = x0
    w_i_min = w0

    local x
    local w
    local s
    
    for i in 1:length(y)
        x[i] ~ Normal(mean = x_i_min, precision = τ_x)
        w[i] ~ Normal(mean = w_i_min, precision = τ_w)
        s[i] := f(x[i], w[i])
        y[i] ~ Normal(mean = s[i], precision = τ_y)
        
        x_i_min = x[i]
        w_i_min = w[i]
    end
    
end

constraints = @constraints begin 
    q(x0, w0, w, τ_x, τ_w, τ_y, s) = q(x0, w, w0, s)q(τ_w)q(τ_x)q(τ_y)
end

m_x_0, τ_x_0 = -20.0, 1.0
m_w_0, τ_w_0 = 20.0, 1.0

# We set relatively strong priors for random walk noise components
# and sort of vague prior for the noise of the observations
a_x, b_x = 0.01, 0.01var(real_x)
a_w, b_w = 0.01, 0.01var(real_w)
a_y, b_y = 1.0, 1.0

# We set relatively strong priors for messages
xinit = map(r -> NormalMeanPrecision(r, τ_x_0), reverse(range(-60, -20, length = n)))
winit = map(r -> NormalMeanPrecision(r, τ_w_0), range(20, 60, length = n))


init = @initialization begin
    μ(x) = xinit
    μ(w) = winit
    q(τ_x) = GammaShapeRate(a_x, b_x)
    q(τ_w) = GammaShapeRate(a_w, b_w)
    q(τ_y) = GammaShapeRate(a_y, b_y)
end

result = infer(
    model = identification_problem(f=+, m_x_0=m_x_0, τ_x_0=τ_x_0, a_x=a_x, b_x=b_x, m_w_0=m_w_0, τ_w_0=τ_w_0, a_w=a_w, b_w=b_w, a_y=a_y, b_y=b_y),
    data  = (y = real_y,), 
    options = (limit_stack_depth = 500, ), 
    constraints = constraints, 
    initialization = init,
    iterations = 50
)

The solution is provide the following constraint, which includes x:

constraints = @constraints begin 
    q(x0, w0, x, w, τ_x, τ_w, τ_y, s) = q(x, x0, w, w0, s)q(τ_w)q(τ_x)q(τ_y)
end