System identification problem

[langestefan]

Hi,

I am trying to solve a system identification problem for an RC network using RxInfer. I have followed Bert's course BMLIP, but it's been 2 years and in the meantime the syntax has changed quite a bit :)

I have the following code, where mΦh and mTa are measured inputs and mx are measured outputs. The model is given by this code:

function transition(T, C, R, Φh, Ta, Δt)
    return (1 - Δt/(C*R))*T + (Δt/C)*[1/R, 1]'*[Ta, Φh]
end

@model function LGDS(mx, mΦh, mTa, Δt)

    # model parameters
    C ~ Gamma(7.7, 0.35)
    R ~ Gamma(1.6, 0.8)

    # initial states
    T[1] ~ Normal(0.0, 0.1)
    mx[1] ~ Normal(T[1], 0.1)

    # dynamics
    for i in eachindex(mx)
        T[i+1] ~ Normal(transition(T[i], C, R, mΦh[i], mTa[i], Δt), 0.1)
    end

    # observations
    for i in eachindex(mx)[2:end]
        mx[i] ~ Normal(T[i], 0.1)
    end        
end

model = LGDS(Δt = t_int)
conditioned = model | (mx = T_i_sample, mΦh = Φ_h_sample, mTa = ones(length(T_i_sample))*Tambient)

For which I then try to infer the parameters C and R. At the moment I am using simulated data so I know exactly what these values are.

result = infer(
    model = LGDS(Δt = t_int),
    data = (mx = T_i_sample, mΦh = Φ_h_sample, mTa = ones(length(T_i_sample))*Tambient),
)

Unfortunately this model is giving me the following error message: Half-edge has been found: T_53. To terminate half-edges 'Uninformative' node can be used.

I must be missing something obvious, could anyone help?

After solving this problem I would also like to:

• Make mΦh and mTa random variables since they are noisy observations
• Scale up the problem to more dimensions, I may need ContinuousTransition if I write it in state-space form (it will remain a system of linear ODE's)

[langestefan]

I was able to solve the previous error message by restating the model like this:

function transition(T, C, R, Φh, Ta, Δt)
    return (1 - Δt/(C*R))*T + (Δt/C)*[1/R, 1]'*[Ta, Φh]
end

@model function LGDS(mx, mΦh, mTa, Δt)
    # model parameters
    C ~ Gamma(7.7, 0.35)
    R ~ Gamma(1.6, 0.8)

    # initial state
    T_prev ~ Normal(μ = 0.0, var = 0.1)

    # dynamics
    for i in 1:length(mx)

        # transition
        T[i] ~ Normal(μ = transition(T_prev, C, R, mΦh[i], mTa[i], Δt), var = 0.1) 

        # observations
        mx[i] ~ Normal(μ = T[i], var = 0.1)
        
        T_prev = T[i]
    end
end

model = LGDS(Δt = t_int)
conditioned = model | (mx = T_i_sample, mΦh = Φ_h_sample, mTa = ones(length(T_i_sample))*Tambient)

Which is now giving the message:
Delta node `transition` requires meta specification with the `where { meta = ... }` in the `@model` macro or with the separate `@meta` specification. See documentation for the `DeltaMeta`.

[failin3]

This error is caused due to not having specified a linearization approximation technique for the transition function.
You should be able to solve that by creating a meta function such as:

LGDS_meta = @meta begin 
    transition() -> Linearization()
end

and passing it with the meta keyword to the infer function

You could use Linearization, Unscented or CVI, these options are shown in the nonlinear sensor fusion example.

[albertpod]

Thanks @failin3! Note that linearization won't work for @langestefan as he uses Gamma priors.
@mhidalgoaraya will sketch a solution that uses CVI.

However, if the reason for Gamma prior around nonlinearity is a matter of positive support, there are some hacks you may use to facilitate inference speed:

# define a function
f(x) = clamp(x, 0.0, Inf)

# in the model
X ~ Normal(mean = 7.0, variance=1e2)
x := f(X)
# rest of the code

# meta specification
LGDS_meta = @meta begin 
    f() -> Linearization()
    transition() -> Unscented()
end

[langestefan]

Thanks guys! The help is very appreciated. Aren't there loops in this model that are problematic? $T_{k} \to T_{k+1}$ with $C,R \to T_{k}$ and $C,R \to T_{k+1}$? As in the Bayesian Linear Regression Tutorial

[albertpod]

Good observation! There are loops, indeed. If you use CVI, the loops around the transition node will be "broken" via mean-field constraint (requirement of CVI) see #299 ; if you proceed with Linearization, then you'd have to initialize the messages yourself.

[mhidalgoaraya]

Hi @langestefan, as @albertpod mentioned this solution can work for your case using CVI. But, keep in mind that the stability of CVI is highly dependent on 1) the choice of hyperparameters (please have a look in the documentation for ProdCVI) and 2) the optimization scheme (i.e., gradient descent).

function transition(T, C, R, Φh, Ta, Δt)
    return (1 - Δt/(C*R))*T .+ (Δt/C)*[1/R, 1]'*[Ta, Φh]
end

@model function LGDS(mx, mΦh, mTa, Δt)
    # model parameters
    C ~ Gamma(7.7, 0.35)
    R ~ Gamma(1.6, 0.8)

    # initial state
    T_prev ~ Normal(μ = 0.0, var = 0.1)
    # dynamics
    for i in 1:length(mx)
        T[i] ~ Normal(μ = transition(T_prev, C, R, mΦh[i], mTa[i], Δt), var = 0.1) 
        mx[i] ~ Normal(μ = T[i], var = 0.1)        
        T_prev = T[i]
    end
end

my_meta_ = @meta begin
    transition() -> CVI(StableRNG(42), 500, 50, Optimisers.Descent(0.0001))
end

init = @initialization begin
    q(T) = NormalMeanVariance(0,1)
    q(T_prev) = NormalMeanVariance(0,1)
    μ(C) = GammaShapeRate(1.0, 1.0)
    μ(R) = GammaShapeRate(1.0, 1.0)
end

result = infer(
    model = LGDS(Δt = t_int,mΦh = Φ_h_sample, mTa = ones(length(T_i_sample))*Tambient, ),
    meta = my_meta_,
    constraints = MeanField(),
    iterations = 10,
    initialization = init,
    data = (mx = T_i_sample,),
    free_energy = true
)```

[langestefan]

@langestefan in the meantime one potential solution can be to try different optimizers and see if the problem disappears. @Nimrais do you have a suggestion for this issue?

Yep, trying that. Seems to run into computational issues as well:

DomainError with -Inf:
Gamma: the condition θ > zero(θ) is not satisfied.

I am trying various initial parameters and optimisers to see if that leads to a solution. Perhaps we can also consider a change of variables for $1/C$, $1/R$ and $1/CR$, perhaps that would that be more stable?

[wouterwln]

Okay, so there's a numerical instability somewhere. Can you recreate the previous error and @show what matrix you're trying to invert? That helps us see if there's a numerical instability somewhere or if we are optimizing outside of the manifold of natural parameters

[langestefan]

I'm not explicitly doing any matrix inversion, these errors are coming from the free energy minimization.

┌ Warning: CVI approximation has not updated the initial state. The method did not converge. Set warn = false to supress this warning.
└ @ ReactiveMP C:\Users\lange.julia\packages\ReactiveMP\hqWuG\src\approximations\cvi.jl:239
DomainError with -Inf:
Gamma: the condition θ > zero(θ) is not satisfied.

[wouterwln]

Sure, but in the free energy minimizing procedure, somewhere a matrix is inverted. (Maybe I should've asked, can you print the matrix WE are trying to invert ;) ) You can do this by dev ing ReactiveMP and @show-ing the matrix that goes into the fastcholsky function on line 180 of cvi.jl in src/approximations (before the function is called and thus the error is thrown of course)

[langestefan]

Ah like that. Okay got it ;)

I got this result (there wasn't a fastcholsky function call on line 180 of cvi.jl)

Fisher = [1.6449340668482262 1.0; 1.0 1.0]
Fisher = [1.6449340668482262 1.0; 1.0 1.0]
Fisher = [1.6449340668482262 1.0; 1.0 1.0]
Fisher = [1.6449340668482262 1.0; 1.0 1.0]
Fisher = [3.0422459703088222e-6 4.542747829467995e-6; 4.542747829467995e-6 6.7833401487808994e-6]
Fisher = [3.0422459703088222e-6 4.542747829467995e-6; 4.542747829467995e-6 6.7833401487808994e-6]
Fisher = [3.0422459703088222e-6 4.542747829467995e-6; 4.542747829467995e-6 6.7833401487808994e-6]
Fisher = [3.0422459703088222e-6 4.542747829467995e-6; 4.542747829467995e-6 6.7833401487808994e-6]
Fisher = [3.0422459703088222e-6 4.542747829467995e-6; 4.542747829467995e-6 6.7833401487808994e-6]
Fisher = [3.0422459703088222e-6 4.542747829467995e-6; 4.542747829467995e-6 6.7833401487808994e-6]
Fisher = [3.0422459703088222e-6 4.542747829467995e-6; 4.542747829467995e-6 6.7833401487808994e-6]
Fisher = [3.0422459703088222e-6 4.542747829467995e-6; 4.542747829467995e-6 6.7833401487808994e-6]
Fisher = [3.0422459703088222e-6 4.542747829467995e-6; 4.542747829467995e-6 6.7833401487808994e-6]
Fisher = [3.0422459703088222e-6 4.542747829467995e-6; 4.542747829467995e-6 6.7833401487808994e-6]
Fisher = [9.900003115875254e-7 1.4349826661851786e-6; 1.4349826661851786e-6 2.0799753772216625e-6]
Fisher = [9.900003115875254e-7 1.4349826661851786e-6; 1.4349826661851786e-6 2.0799753772216625e-6]
Fisher = [9.00938836660723e-9 1.3047641513161003e-8; 1.3047641513161003e-8 1.8895949746581453e-8]
Fisher = [2.3212945371792685e-10 3.361169646747688e-10; 3.361169646747688e-10 4.866879758075584e-10]
Fisher = [2.2085489619305654e-11 3.1981962481804264e-11; 3.1981962481804264e-11 4.6313029134963225e-11]
Fisher = [2.2085489619305654e-11 3.1981962481804264e-11; 3.1981962481804264e-11 4.6313029134963225e-11]
Fisher = [1.4656656224688926e-13 2.1224209902940764e-13; 2.1224209902940764e-13 3.073464227436245e-13]
Fisher = [1.4656656224688926e-13 2.1224209902940764e-13; 2.1224209902940764e-13 3.073464227436245e-13]
Fisher = [2.3359658853317013e-15 3.3826955480231883e-15; 3.3826955480231883e-15 4.898457311584875e-15]
Fisher = [1.6268999519195686e-17 2.3559022640171636e-17; 2.3559022640171636e-17 3.4115653338439544e-17]
MethodError: no method matching cholesky(::Type{PositiveFactorizations.Positive}, ::Cholesky{Float64, StaticArraysCore.SMatrix{2, 2, Float64, 4}}; tol::Float64)

[HoangMHNguyen]

@langestefan Here is another way to run inference on your model with Linearization() or UnscentedTransform() approximation. In your model, I assume that you use Gamma distribution for R and C because of positiveness. Now we can transform R and C to log domain and use exp function inside your transition function as follows:

function my_transition(T, logC, logR, Φh, Ta, Δt)
    return (1 - Δt/(exp(logC)*exp(logR)))*T + (Δt/exp(logC))*[1/exp(logR), 1]'*[Ta, Φh]
end

Now we can specify Normal priors for logC and logR and the model now slightly become:

@model function LGDS(mx, mΦh, mTa, Δt)
    # model parameters
    lC ~ NormalMeanVariance(0, 50)
    lR ~ NormalMeanVariance(0, 50)

    # initial state
    T_prev ~ Normal(μ = 0.0, var = 0.1)
    # dynamics
    for i in 1:length(mx)
        # transition
        T[i] ~ Normal(μ =  my_transition(T_prev, lC, lR, mΦh[i], mTa[i], Δt), var = 0.1)

        # observations
        mx[i] ~ Normal(μ = T[i], var = 0.1)
        
        T_prev = T[i]
    end
end

my_meta = @meta begin
    # my_transition() -> Linearization() #alternative method
    my_transition() -> UT()
end

Since your model contain loop, we need to initialize messages for lC and lR before doing inference

init = @initialization begin
    μ(lC) =NormalMeanVariance(log(9), 10)
    μ(lR) =NormalMeanVariance(log(20), 10)
end

result = infer(
    model = LGDS(Δt = t_int,mΦh = Φ_h_sample, mTa = ones(length(T_i_sample))*Tambient, ),
    meta = my_meta,
    iterations = 51,
    initialization = init,
    data = (mx = T_i_sample,),
    free_energy = true
)

Note that you will expect an oscilating behavior in the BFE since we do loopy Belief Propagation here.

[langestefan]

Hi @HoangMHNguyen,

Thanks for your elaborate answers. To be honest I still don't fully grasp what is going on in the graph above, because I thought lR and lC were variables and not factors.

For example, in the work of Loeliger (2007) on FFG's it is very clear that C is an unkown variable to be estimated and not a factor:

[HoangMHNguyen]

Hi @langestefan!

Sorry for the confusion. Because you gave the graph as in the question so my answer was based on that graph. In RxInfer we indeed use Forney-style Factor graph (FFG) as in Loeliger's paper, and you are right that lC and lR are variables and on FFG they are edges, not factor nodes. So when you initialize messages for these variables, these messages will flow to adjacent nodes. Maybe @albertpod can elaborate on this a bit?

[albertpod]

@langestefan technically, lR and LC are not factors but variables and indeed must be displayed as edges (as in the paper from Loeliger). However, we don't have the functionality to draw proper FFG from GraphPPL.jl. It's on the roadmap, but there are no ETAs for that.

In proper FFG, the lR and lC will be represented by edges shared by the equality node. So, think of these initializations as outgoing messages from the equality node.

[langestefan]

I see, that answers a lot of my questions about how to interpret the visualization, thanks. I will make an FFG of the model above and share it here

[langestefan]

The full dynamics are described by:

$$ \begin{align} \mathbf{T}_{k+1} &= \mathbf{A_d}\mathbf{T}_k + \mathbf{B_d}\mathbf{u}_k + \mathbf {w}_k\\ \mathbf{x}_k &= \mathbf{C}\mathbf{T}_k + \mathbf {v}_k \end{align} $$

Which if I am correct is described by this FFG:

(I wrote 'prior' at the top and bottom left)

[langestefan]

Hi again @albertpod @HoangMHNguyen.

I have made quite some progress but I am now a bit stuck getting it to converge. I'm assuming I need to either figure out how to tweak the hyperparameters, which is quite difficult given how many there are, or I need to figure out a better way to specify my model. As you can see there are some matrix multiplications, which I expect is where the problems come from.

The model is now of this form, with 3 state variables and 2 output variables:

$$ \begin{align} \mathbf{T}_{k+1} &= \mathbf{A_d}\mathbf{T}_k + \mathbf{B_d}\mathbf{u}_k + \mathbf {w}_k\\ \mathbf{x}_k &= \mathbf{C}\mathbf{T}_k + \mathbf {v}_k \end{align} $$

(I have more than triple checked all the equations in the code, they should be correct.)

The error is always the same:

Failed to compute node bound free energy component. The result is NaN.
Use objective_diagnostics keyword argument in the inference function to suppress this error.

Stacktrace

{
	"name": "ErrorException",
	"message": "    Failed to compute node bound free energy component. The result is `NaN`. 
    Use `objective_diagnostics` keyword argument in the `inference` function to suppress this error.
    fform = MvNormalMeanCovariance, neighbors = Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(T_prev_31, out, NodeData in context  with properties name = T_prev, index = nothing with extra: {:bfe_stream = ProxyObservable(CountingReal, MapProxy(Marginal)), :rmp_variable = RandomVariable(ReactiveMP.MessageObservable{AbstractMessage}[ReactiveMP.MessageObservable{AbstractMessage}(RecentSubject(AbstractMessage, Subject{AbstractMessage, AsapScheduler, AsapScheduler}), LazyObservable(AbstractMessage)), ReactiveMP.MessageObservable{AbstractMessage}(RecentSubject(AbstractMessage, Subject{AbstractMessage, AsapScheduler, AsapScheduler}), LazyObservable(AbstractMessage))], ReactiveMP.MessageObservable{Message}[ReactiveMP.MessageObservable{Message}(RecentSubject(Message, Subject{Message, AsapScheduler, AsapScheduler}), LazyObservable(Message)), ReactiveMP.MessageObservable{Message}(RecentSubject(Message, Subject{Message, AsapScheduler, AsapScheduler}), LazyObservable(Message))], ReactiveMP.MarginalObservable(RecentSubject(Marginal, Subject{Marginal, AsapScheduler, AsapScheduler}), LazyObservable(Marginal)))}), (constvar_32_33, μ, NodeData in context  with properties name = constvar_32, index = nothing with extra: {:rmp_variable = ConstVariable(ReactiveMP.MarginalObservable(RecentSubject(Marginal, Subject{Marginal, AsapScheduler, AsapScheduler}), LazyObservable(Marginal)), ReactiveMP.MessageObservable{AbstractMessage}(RecentSubject(AbstractMessage, Subject{AbstractMessage, AsapScheduler, AsapScheduler}), LazyObservable(AbstractMessage)), [293.15, 293.15, 293.15], 1)}), (constvar_34_35, Σ, NodeData in context  with properties name = constvar_34, index = nothing with extra: {:rmp_variable = ConstVariable(ReactiveMP.MarginalObservable(RecentSubject(Marginal, Subject{Marginal, AsapScheduler, AsapScheduler}), LazyObservable(Marginal)), ReactiveMP.MessageObservable{AbstractMessage}(RecentSubject(AbstractMessage, Subject{AbstractMessage, AsapScheduler, AsapScheduler}), LazyObservable(AbstractMessage)), [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0], 1)})]
",
	"stack": "    Failed to compute node bound free energy component. The result is `NaN`. 
    Use `objective_diagnostics` keyword argument in the `inference` function to suppress this error.
    fform = MvNormalMeanCovariance, neighbors = Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(T_prev_31, out, NodeData in context  with properties name = T_prev, index = nothing with extra: {:bfe_stream = ProxyObservable(CountingReal, MapProxy(Marginal)), :rmp_variable = RandomVariable(ReactiveMP.MessageObservable{AbstractMessage}[ReactiveMP.MessageObservable{AbstractMessage}(RecentSubject(AbstractMessage, Subject{AbstractMessage, AsapScheduler, AsapScheduler}), LazyObservable(AbstractMessage)), ReactiveMP.MessageObservable{AbstractMessage}(RecentSubject(AbstractMessage, Subject{AbstractMessage, AsapScheduler, AsapScheduler}), LazyObservable(AbstractMessage))], ReactiveMP.MessageObservable{Message}[ReactiveMP.MessageObservable{Message}(RecentSubject(Message, Subject{Message, AsapScheduler, AsapScheduler}), LazyObservable(Message)), ReactiveMP.MessageObservable{Message}(RecentSubject(Message, Subject{Message, AsapScheduler, AsapScheduler}), LazyObservable(Message))], ReactiveMP.MarginalObservable(RecentSubject(Marginal, Subject{Marginal, AsapScheduler, AsapScheduler}), LazyObservable(Marginal)))}), (constvar_32_33, μ, NodeData in context  with properties name = constvar_32, index = nothing with extra: {:rmp_variable = ConstVariable(ReactiveMP.MarginalObservable(RecentSubject(Marginal, Subject{Marginal, AsapScheduler, AsapScheduler}), LazyObservable(Marginal)), ReactiveMP.MessageObservable{AbstractMessage}(RecentSubject(AbstractMessage, Subject{AbstractMessage, AsapScheduler, AsapScheduler}), LazyObservable(AbstractMessage)), [293.15, 293.15, 293.15], 1)}), (constvar_34_35, Σ, NodeData in context  with properties name = constvar_34, index = nothing with extra: {:rmp_variable = ConstVariable(ReactiveMP.MarginalObservable(RecentSubject(Marginal, Subject{Marginal, AsapScheduler, AsapScheduler}), LazyObservable(Marginal)), ReactiveMP.MessageObservable{AbstractMessage}(RecentSubject(AbstractMessage, Subject{AbstractMessage, AsapScheduler, AsapScheduler}), LazyObservable(AbstractMessage)), [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0], 1)})]


Stacktrace:
    [1] try_yieldto(undo::typeof(Base.ensure_rescheduled))
      @ Base .\\task.jl:931
    [2] wait()
      @ Base .\\task.jl:995
    [3] wait(c::Condition; first::Bool)
      @ Base .\\condition.jl:130
    [4] wait
      @ .\\condition.jl:125 [inlined]
    [5] limitstack(callback::RxInfer.var\"#7#8\"{Rocket.RecentSubjectInstance{AbstractMessage, Subject{AbstractMessage, AsapScheduler, AsapScheduler}}, DefferedMessage{Tuple{ReactiveMP.MessageObservable{Message}}, Tuple{ReactiveMP.MarginalObservable}, ReactiveMP.MessageMapping{DeltaFn{typeof(transition)}, Tuple{Val{:in}, Int64}, Marginalisation, Val{(:in,)}, Val{(:ins,)}, DeltaMeta{Unscented{Float64, Float64, Float64, Nothing}, Nothing}, Nothing, DeltaFnNode{typeof(transition), FixedArguments.FixedCallable{typeof(transition), Tuple{FixedArguments.FixedArgument{FixedArguments.FixedPosition{5}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{6}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{7}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{8}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{9}, ConstVariable}}, typeof(ReactiveMP.__unpack_latest_static)}, 4, Tuple{FixedArguments.FixedArgument{FixedArguments.FixedPosition{5}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{6}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{7}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{8}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{9}, ConstVariable}}, ReactiveMP.FactorNodeLocalClusters{Tuple{ReactiveMP.FactorNodeLocalMarginal, ReactiveMP.FactorNodeLocalMarginal}, Nothing}}}}}, instance::RxInfer.LimitStackScheduler)
      @ RxInfer C:\\Users\\lange\\.julia\\packages\\RxInfer\\AszMu\\src\\rocket.jl:64
    [6] scheduled_next!
      @ C:\\Users\\lange\\.julia\\packages\\RxInfer\\AszMu\\src\\rocket.jl:86 [inlined]
    [7] actor_on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:255 [inlined]
    [8] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:203 [inlined]
    [9] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\scheduled.jl:12 [inlined]
   [10] next!(actor::Rocket.ScheduledActor{AbstractMessage, RxInfer.LimitStackScheduler, Rocket.RecentSubjectInstance{AbstractMessage, Subject{AbstractMessage, AsapScheduler, AsapScheduler}}}, data::DefferedMessage{Tuple{ReactiveMP.MessageObservable{Message}}, Tuple{ReactiveMP.MarginalObservable}, ReactiveMP.MessageMapping{DeltaFn{typeof(transition)}, Tuple{Val{:in}, Int64}, Marginalisation, Val{(:in,)}, Val{(:ins,)}, DeltaMeta{Unscented{Float64, Float64, Float64, Nothing}, Nothing}, Nothing, DeltaFnNode{typeof(transition), FixedArguments.FixedCallable{typeof(transition), Tuple{FixedArguments.FixedArgument{FixedArguments.FixedPosition{5}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{6}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{7}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{8}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{9}, ConstVariable}}, typeof(ReactiveMP.__unpack_latest_static)}, 4, Tuple{FixedArguments.FixedArgument{FixedArguments.FixedPosition{5}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{6}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{7}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{8}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{9}, ConstVariable}}, ReactiveMP.FactorNodeLocalClusters{Tuple{ReactiveMP.FactorNodeLocalMarginal, ReactiveMP.FactorNodeLocalMarginal}, Nothing}}}})
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:206
   [11] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\operators\\map.jl:62 [inlined]
   [12] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:206 [inlined]
--- the last 2 lines are repeated 1 more time ---
   [15] next_received!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined.jl:101 [inlined]
   [16] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined.jl:68 [inlined]
   [17] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:206 [inlined]
--- the last 3 lines are repeated 1 more time ---
   [21] next_received!(wrapper::Rocket.CombineLatestUpdatesActorWrapper{Tuple{ReactiveMP.MarginalObservable}, Rocket.CombineLatestInnerActor{Tuple{ReactiveMP.MarginalObservable}, Rocket.CombineLatestActorWrapper{Rocket.MStorage2{Tuple{ReactiveMP.MessageObservable{Message}}, Tuple{ReactiveMP.MarginalObservable}}, Rocket.CombineLatestInnerActor{Tuple{Tuple{ReactiveMP.MessageObservable{Message}}, Tuple{ReactiveMP.MarginalObservable}}, Rocket.CombineLatestActorWrapper{Rocket.MStorage2{Tuple{Tuple{ReactiveMP.MessageObservable{Message}}, Tuple{ReactiveMP.MarginalObservable}}, NTuple{5, AbstractMessage}}, Rocket.MapActor{Tuple{Tuple{Tuple{ReactiveMP.MessageObservable{Message}}, Tuple{ReactiveMP.MarginalObservable}}, NTuple{5, AbstractMessage}}, Rocket.MapActor{Tuple{Tuple{ReactiveMP.MessageObservable{Message}}, Tuple{ReactiveMP.MarginalObservable}}, Rocket.ScheduledActor{AbstractMessage, RxInfer.LimitStackScheduler, Rocket.RecentSubjectInstance{AbstractMessage, Subject{AbstractMessage, AsapScheduler, AsapScheduler}}}, ReactiveMP.var\"#507#509\"{ReactiveMP.MessageMapping{DeltaFn{typeof(transition)}, Tuple{Val{:in}, Int64}, Marginalisation, Val{(:in,)}, Val{(:ins,)}, DeltaMeta{Unscented{Float64, Float64, Float64, Nothing}, Nothing}, Nothing, DeltaFnNode{typeof(transition), FixedArguments.FixedCallable{typeof(transition), Tuple{FixedArguments.FixedArgument{FixedArguments.FixedPosition{5}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{6}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{7}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{8}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{9}, ConstVariable}}, typeof(ReactiveMP.__unpack_latest_static)}, 4, Tuple{FixedArguments.FixedArgument{FixedArguments.FixedPosition{5}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{6}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{7}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{8}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{9}, ConstVariable}}, ReactiveMP.FactorNodeLocalClusters{Tuple{ReactiveMP.FactorNodeLocalMarginal, ReactiveMP.FactorNodeLocalMarginal}, Nothing}}}}}, typeof(first)}, PushEach, Rocket.UInt8UpdatesStatus}}, PushNew, Rocket.UInt8UpdatesStatus}}, PushNew, Rocket.UInt8UpdatesStatus, typeof(identity), Nothing}, data::Marginal{JointNormal{MvNormalMeanCovariance{Float64, Vector{Float64}, Matrix{Float64}}, NTuple{4, Tuple{Int64}}}, Nothing}, index::Int64)
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined_updates.jl:72
   [22] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined_updates.jl:34 [inlined]
   [23] scheduled_next!(actor::Rocket.CombineLatestUpdatesInnerActor{Marginal, Rocket.CombineLatestUpdatesActorWrapper{Tuple{ReactiveMP.MarginalObservable}, Rocket.CombineLatestInnerActor{Tuple{ReactiveMP.MarginalObservable}, Rocket.CombineLatestActorWrapper{Rocket.MStorage2{Tuple{ReactiveMP.MessageObservable{Message}}, Tuple{ReactiveMP.MarginalObservable}}, Rocket.CombineLatestInnerActor{Tuple{Tuple{ReactiveMP.MessageObservable{Message}}, Tuple{ReactiveMP.MarginalObservable}}, Rocket.CombineLatestActorWrapper{Rocket.MStorage2{Tuple{Tuple{ReactiveMP.MessageObservable{Message}}, Tuple{ReactiveMP.MarginalObservable}}, NTuple{5, AbstractMessage}}, Rocket.MapActor{Tuple{Tuple{Tuple{ReactiveMP.MessageObservable{Message}}, Tuple{ReactiveMP.MarginalObservable}}, NTuple{5, AbstractMessage}}, Rocket.MapActor{Tuple{Tuple{ReactiveMP.MessageObservable{Message}}, Tuple{ReactiveMP.MarginalObservable}}, Rocket.ScheduledActor{AbstractMessage, RxInfer.LimitStackScheduler, Rocket.RecentSubjectInstance{AbstractMessage, Subject{AbstractMessage, AsapScheduler, AsapScheduler}}}, ReactiveMP.var\"#507#509\"{ReactiveMP.MessageMapping{DeltaFn{typeof(transition)}, Tuple{Val{:in}, Int64}, Marginalisation, Val{(:in,)}, Val{(:ins,)}, DeltaMeta{Unscented{Float64, Float64, Float64, Nothing}, Nothing}, Nothing, DeltaFnNode{typeof(transition), FixedArguments.FixedCallable{typeof(transition), Tuple{FixedArguments.FixedArgument{FixedArguments.FixedPosition{5}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{6}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{7}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{8}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{9}, ConstVariable}}, typeof(ReactiveMP.__unpack_latest_static)}, 4, Tuple{FixedArguments.FixedArgument{FixedArguments.FixedPosition{5}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{6}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{7}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{8}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{9}, ConstVariable}}, ReactiveMP.FactorNodeLocalClusters{Tuple{ReactiveMP.FactorNodeLocalMarginal, ReactiveMP.FactorNodeLocalMarginal}, Nothing}}}}}, typeof(first)}, PushEach, Rocket.UInt8UpdatesStatus}}, PushNew, Rocket.UInt8UpdatesStatus}}, PushNew, Rocket.UInt8UpdatesStatus, typeof(identity), Nothing}}, value::Marginal{JointNormal{MvNormalMeanCovariance{Float64, Vector{Float64}, Matrix{Float64}}, NTuple{4, Tuple{Int64}}}, Nothing}, ::AsapScheduler)
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\schedulers\\asap.jl:23
   [24] on_next!(subject::Subject{Marginal, AsapScheduler, AsapScheduler}, data::Marginal{JointNormal{MvNormalMeanCovariance{Float64, Vector{Float64}, Matrix{Float64}}, NTuple{4, Tuple{Int64}}}, Nothing})
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\subjects\\subject.jl:62
   [25] actor_on_next!(::BaseActorTrait{Marginal}, actor::Subject{Marginal, AsapScheduler, AsapScheduler}, data::Marginal{JointNormal{MvNormalMeanCovariance{Float64, Vector{Float64}, Matrix{Float64}}, NTuple{4, Tuple{Int64}}}, Nothing})
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:250
   [26] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:202 [inlined]
   [27] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\subjects\\recent.jl:62 [inlined]
   [28] actor_on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:250 [inlined]
   [29] next!(actor::Rocket.RecentSubjectInstance{Marginal, Subject{Marginal, AsapScheduler, AsapScheduler}}, data::Marginal{JointNormal{MvNormalMeanCovariance{Float64, Vector{Float64}, Matrix{Float64}}, NTuple{4, Tuple{Int64}}}, Nothing})
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:202
   [30] next_received!(wrapper::Rocket.CombineLatestUpdatesActorWrapper{Tuple{ProxyObservable{Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.ManyOfObservable{Rocket.CombineLatestUpdatesObservable{NTuple{4, ReactiveMP.MessageObservable{Message}}, NTuple{4, ReactiveMP.MessageObservable{Message}}, PushNew, typeof(identity), Nothing}}}, Rocket.CombineLatestObservable{Tuple{Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.ManyOfObservable{Rocket.CombineLatestUpdatesObservable{NTuple{4, ReactiveMP.MessageObservable{Message}}, NTuple{4, ReactiveMP.MessageObservable{Message}}, PushNew, typeof(identity), Nothing}}}, NTuple{5, AbstractMessage}}, Tuple{Rocket.CombineLatestUpdatesObservable{Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.ManyOfObservable{Rocket.CombineLatestUpdatesObservable{NTuple{4, ReactiveMP.MessageObservable{Message}}, NTuple{4, ReactiveMP.MessageObservable{Message}}, PushNew, typeof(identity), Nothing}}}, Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.ManyOfObservable{Rocket.CombineLatestUpdatesObservable{NTuple{4, ReactiveMP.MessageObservable{Message}}, NTuple{4, ReactiveMP.MessageObservable{Message}}, PushNew, typeof(identity), Nothing}}}, PushNew, typeof(identity), Nothing}, Rocket.CombineLatestObservable{NTuple{5, AbstractMessage}, NTuple{5, ReactiveMP.MessageObservable{AbstractMessage}}, PushNew}}, PushEach}, Rocket.MapProxy{Tuple{Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.ManyOfObservable{Rocket.CombineLatestUpdatesObservable{NTuple{4, ReactiveMP.MessageObservable{Message}}, NTuple{4, ReactiveMP.MessageObservable{Message}}, PushNew, typeof(identity), Nothing}}}, NTuple{5, AbstractMessage}}, typeof(first)}}, ProxyObservable{Nothing, Rocket.CombineLatestObservable{Tuple{Nothing, NTuple{5, AbstractMessage}}, Tuple{SingleObservable{Nothing, AsapScheduler}, Rocket.CombineLatestObservable{NTuple{5, AbstractMessage}, NTuple{5, ReactiveMP.MessageObservable{AbstractMessage}}, PushNew}}, PushEach}, Rocket.MapProxy{Tuple{Nothing, NTuple{5, AbstractMessage}}, typeof(first)}}}, Rocket.RecentSubjectInstance{Marginal, Subject{Marginal, AsapScheduler, AsapScheduler}}, PushNew, Rocket.UInt8UpdatesStatus, ReactiveMP.MarginalMapping{DeltaFn{typeof(transition)}, Val{:ins}, Val{(:out, :ins)}, Nothing, DeltaMeta{Unscented{Float64, Float64, Float64, Nothing}, Nothing}, DeltaFnNode{typeof(transition), FixedArguments.FixedCallable{typeof(transition), Tuple{FixedArguments.FixedArgument{FixedArguments.FixedPosition{5}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{6}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{7}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{8}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{9}, ConstVariable}}, typeof(ReactiveMP.__unpack_latest_static)}, 4, Tuple{FixedArguments.FixedArgument{FixedArguments.FixedPosition{5}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{6}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{7}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{8}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{9}, ConstVariable}}, ReactiveMP.FactorNodeLocalClusters{Tuple{ReactiveMP.FactorNodeLocalMarginal, ReactiveMP.FactorNodeLocalMarginal}, Nothing}}}, typeof(ReactiveMP.reset_vstatus)}, data::Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.ManyOfObservable{Rocket.CombineLatestUpdatesObservable{NTuple{4, ReactiveMP.MessageObservable{Message}}, NTuple{4, ReactiveMP.MessageObservable{Message}}, PushNew, typeof(identity), Nothing}}}, index::Int64)
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined_updates.jl:72
   [31] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined_updates.jl:34 [inlined]
   [32] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:206 [inlined]
   [33] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\operators\\map.jl:62 [inlined]
   [34] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:206 [inlined]
   [35] next_received!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined.jl:101 [inlined]
   [36] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined.jl:68 [inlined]
   [37] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:206 [inlined]
   [38] next_received!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined_updates.jl:72 [inlined]
   [39] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined_updates.jl:34 [inlined]
   [40] scheduled_next!(actor::Rocket.CombineLatestUpdatesInnerActor{Message, Rocket.CombineLatestUpdatesActorWrapper{Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.ManyOfObservable{Rocket.CombineLatestUpdatesObservable{NTuple{4, ReactiveMP.MessageObservable{Message}}, NTuple{4, ReactiveMP.MessageObservable{Message}}, PushNew, typeof(identity), Nothing}}}, Rocket.CombineLatestInnerActor{Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.ManyOfObservable{Rocket.CombineLatestUpdatesObservable{NTuple{4, ReactiveMP.MessageObservable{Message}}, NTuple{4, ReactiveMP.MessageObservable{Message}}, PushNew, typeof(identity), Nothing}}}, Rocket.CombineLatestActorWrapper{Rocket.MStorage2{Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.ManyOfObservable{Rocket.CombineLatestUpdatesObservable{NTuple{4, ReactiveMP.MessageObservable{Message}}, NTuple{4, ReactiveMP.MessageObservable{Message}}, PushNew, typeof(identity), Nothing}}}, NTuple{5, AbstractMessage}}, Rocket.MapActor{Tuple{Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.ManyOfObservable{Rocket.CombineLatestUpdatesObservable{NTuple{4, ReactiveMP.MessageObservable{Message}}, NTuple{4, ReactiveMP.MessageObservable{Message}}, PushNew, typeof(identity), Nothing}}}, NTuple{5, AbstractMessage}}, Rocket.CombineLatestUpdatesInnerActor{Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.ManyOfObservable{Rocket.CombineLatestUpdatesObservable{NTuple{4, ReactiveMP.MessageObservable{Message}}, NTuple{4, ReactiveMP.MessageObservable{Message}}, PushNew, typeof(identity), Nothing}}}, Rocket.CombineLatestUpdatesActorWrapper{Tuple{ProxyObservable{Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.ManyOfObservable{Rocket.CombineLatestUpdatesObservable{NTuple{4, ReactiveMP.MessageObservable{Message}}, NTuple{4, ReactiveMP.MessageObservable{Message}}, PushNew, typeof(identity), Nothing}}}, Rocket.CombineLatestObservable{Tuple{Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.ManyOfObservable{Rocket.CombineLatestUpdatesObservable{NTuple{4, ReactiveMP.MessageObservable{Message}}, NTuple{4, ReactiveMP.MessageObservable{Message}}, PushNew, typeof(identity), Nothing}}}, NTuple{5, AbstractMessage}}, Tuple{Rocket.CombineLatestUpdatesObservable{Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.ManyOfObservable{Rocket.CombineLatestUpdatesObservable{NTuple{4, ReactiveMP.MessageObservable{Message}}, NTuple{4, ReactiveMP.MessageObservable{Message}}, PushNew, typeof(identity), Nothing}}}, Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.ManyOfObservable{Rocket.CombineLatestUpdatesObservable{NTuple{4, ReactiveMP.MessageObservable{Message}}, NTuple{4, ReactiveMP.MessageObservable{Message}}, PushNew, typeof(identity), Nothing}}}, PushNew, typeof(identity), Nothing}, Rocket.CombineLatestObservable{NTuple{5, AbstractMessage}, NTuple{5, ReactiveMP.MessageObservable{AbstractMessage}}, PushNew}}, PushEach}, Rocket.MapProxy{Tuple{Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.ManyOfObservable{Rocket.CombineLatestUpdatesObservable{NTuple{4, ReactiveMP.MessageObservable{Message}}, NTuple{4, ReactiveMP.MessageObservable{Message}}, PushNew, typeof(identity), Nothing}}}, NTuple{5, AbstractMessage}}, typeof(first)}}, ProxyObservable{Nothing, Rocket.CombineLatestObservable{Tuple{Nothing, NTuple{5, AbstractMessage}}, Tuple{SingleObservable{Nothing, AsapScheduler}, Rocket.CombineLatestObservable{NTuple{5, AbstractMessage}, NTuple{5, ReactiveMP.MessageObservable{AbstractMessage}}, PushNew}}, PushEach}, Rocket.MapProxy{Tuple{Nothing, NTuple{5, AbstractMessage}}, typeof(first)}}}, Rocket.RecentSubjectInstance{Marginal, Subject{Marginal, AsapScheduler, AsapScheduler}}, PushNew, Rocket.UInt8UpdatesStatus, ReactiveMP.MarginalMapping{DeltaFn{typeof(transition)}, Val{:ins}, Val{(:out, :ins)}, Nothing, DeltaMeta{Unscented{Float64, Float64, Float64, Nothing}, Nothing}, DeltaFnNode{typeof(transition), FixedArguments.FixedCallable{typeof(transition), Tuple{FixedArguments.FixedArgument{FixedArguments.FixedPosition{5}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{6}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{7}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{8}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{9}, ConstVariable}}, typeof(ReactiveMP.__unpack_latest_static)}, 4, Tuple{FixedArguments.FixedArgument{FixedArguments.FixedPosition{5}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{6}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{7}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{8}, ConstVariable}, FixedArguments.FixedArgument{FixedArguments.FixedPosition{9}, ConstVariable}}, ReactiveMP.FactorNodeLocalClusters{Tuple{ReactiveMP.FactorNodeLocalMarginal, ReactiveMP.FactorNodeLocalMarginal}, Nothing}}}, typeof(ReactiveMP.reset_vstatus)}}, typeof(first)}, PushEach, Rocket.UInt8UpdatesStatus}}, PushNew, Rocket.UInt8UpdatesStatus, typeof(identity), Nothing}}, value::Message{MvNormalMeanCovariance{Float64, Vector{Float64}, Matrix{Float64}}, Nothing}, ::AsapScheduler)
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\schedulers\\asap.jl:23
   [41] on_next!(subject::Subject{Message, AsapScheduler, AsapScheduler}, data::Message{MvNormalMeanCovariance{Float64, Vector{Float64}, Matrix{Float64}}, Nothing})
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\subjects\\subject.jl:62
   [42] actor_on_next!(::BaseActorTrait{Message}, actor::Subject{Message, AsapScheduler, AsapScheduler}, data::Message{MvNormalMeanCovariance{Float64, Vector{Float64}, Matrix{Float64}}, Nothing})
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:250
   [43] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:202 [inlined]
   [44] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\subjects\\recent.jl:62 [inlined]
   [45] actor_on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:250 [inlined]
   [46] next!(actor::Rocket.RecentSubjectInstance{Message, Subject{Message, AsapScheduler, AsapScheduler}}, data::Message{MvNormalMeanCovariance{Float64, Vector{Float64}, Matrix{Float64}}, Nothing})
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:202
   [47] next_received!(wrapper::Rocket.CollectLatestObservableWrapper{AbstractMessage, Rocket.RecentSubjectInstance{Message, Subject{Message, AsapScheduler, AsapScheduler}}, Vector{AbstractMessage}, BitVector, Vector{Teardown}, ReactiveMP.var\"#15#17\"{GenericProd, UnspecifiedFormConstraint}, typeof(ReactiveMP.reset_vstatus)}, data::DefferedMessage{Tuple{ReactiveMP.MessageObservable{Message}}, Tuple{ReactiveMP.MarginalObservable}, ReactiveMP.MessageMapping{MvNormalMeanPrecision, Val{:μ}, Marginalisation, Val{(:out,)}, Val{(:Λ,)}, Nothing, Nothing, Nothing}}, index::CartesianIndex{1})
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\collected.jl:104
   [48] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\collected.jl:93 [inlined]
   [49] scheduled_next!(actor::Rocket.CollectLatestObservableInnerActor{AbstractMessage, CartesianIndex{1}, Rocket.CollectLatestObservableWrapper{AbstractMessage, Rocket.RecentSubjectInstance{Message, Subject{Message, AsapScheduler, AsapScheduler}}, Vector{AbstractMessage}, BitVector, Vector{Teardown}, ReactiveMP.var\"#15#17\"{GenericProd, UnspecifiedFormConstraint}, typeof(ReactiveMP.reset_vstatus)}}, value::DefferedMessage{Tuple{ReactiveMP.MessageObservable{Message}}, Tuple{ReactiveMP.MarginalObservable}, ReactiveMP.MessageMapping{MvNormalMeanPrecision, Val{:μ}, Marginalisation, Val{(:out,)}, Val{(:Λ,)}, Nothing, Nothing, Nothing}}, ::AsapScheduler)
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\schedulers\\asap.jl:23
   [50] on_next!(subject::Subject{AbstractMessage, AsapScheduler, AsapScheduler}, data::DefferedMessage{Tuple{ReactiveMP.MessageObservable{Message}}, Tuple{ReactiveMP.MarginalObservable}, ReactiveMP.MessageMapping{MvNormalMeanPrecision, Val{:μ}, Marginalisation, Val{(:out,)}, Val{(:Λ,)}, Nothing, Nothing, Nothing}})
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\subjects\\subject.jl:62
   [51] actor_on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:250 [inlined]
   [52] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:202 [inlined]
   [53] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\subjects\\recent.jl:62 [inlined]
   [54] #7
      @ C:\\Users\\lange\\.julia\\packages\\RxInfer\\AszMu\\src\\rocket.jl:86 [inlined]
   [55] limitstack(callback::RxInfer.var\"#7#8\"{Rocket.RecentSubjectInstance{AbstractMessage, Subject{AbstractMessage, AsapScheduler, AsapScheduler}}, DefferedMessage{Tuple{ReactiveMP.MessageObservable{Message}}, Tuple{ReactiveMP.MarginalObservable}, ReactiveMP.MessageMapping{MvNormalMeanPrecision, Val{:μ}, Marginalisation, Val{(:out,)}, Val{(:Λ,)}, Nothing, Nothing, Nothing}}}, instance::RxInfer.LimitStackScheduler)
      @ RxInfer C:\\Users\\lange\\.julia\\packages\\RxInfer\\AszMu\\src\\rocket.jl:52
   [56] scheduled_next!
      @ C:\\Users\\lange\\.julia\\packages\\RxInfer\\AszMu\\src\\rocket.jl:86 [inlined]
   [57] actor_on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:255 [inlined]
   [58] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:203 [inlined]
   [59] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\scheduled.jl:12 [inlined]
--- the last 28 lines are repeated 1 more time ---
--- the last 78 lines are repeated 48 more times ---
 [3832] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:206 [inlined]
 [3833] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\operators\\map.jl:62 [inlined]
 [3834] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:206 [inlined]
 [3835] next_received!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined.jl:101 [inlined]
 [3836] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined.jl:68 [inlined]
 [3837] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:206 [inlined]
 [3838] next_received!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined_updates.jl:72 [inlined]
 [3839] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined_updates.jl:34 [inlined]
 [3840] scheduled_next!(actor::Rocket.CombineLatestUpdatesInnerActor{Message, Rocket.CombineLatestUpdatesActorWrapper{Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.MessageObservable{AbstractMessage}}, Rocket.CombineLatestInnerActor{Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.MessageObservable{AbstractMessage}}, Rocket.CombineLatestActorWrapper{Rocket.MStorage2{Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.MessageObservable{AbstractMessage}}, Nothing}, Rocket.MapActor{Tuple{Tuple{ReactiveMP.MessageObservable{Message}, ReactiveMP.MessageObservable{AbstractMessage}}, Nothing}, Rocket.ScheduledActor{AbstractMessage, RxInfer.LimitStackScheduler, Rocket.RecentSubjectInstance{AbstractMessage, Subject{AbstractMessage, AsapScheduler, AsapScheduler}}}, ReactiveMP.var\"#178#181\"{ReactiveMP.MessageMapping{typeof(*), Val{:in}, Marginalisation, Val{(:out, :A)}, Nothing, MatrixCorrectionTools.ReplaceZeroDiagonalEntries{TinyHugeNumbers.TinyNumber}, Nothing, Nothing}}}, PushNew, Rocket.UInt8UpdatesStatus}}, PushNew, Rocket.UInt8UpdatesStatus, typeof(identity), Nothing}}, value::Message{MvNormalMeanPrecision{Float64, Vector{Float64}, Matrix{Float64}}, Nothing}, ::AsapScheduler)
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\schedulers\\asap.jl:23
 [3841] on_next!(subject::Subject{Message, AsapScheduler, AsapScheduler}, data::Message{MvNormalMeanPrecision{Float64, Vector{Float64}, Matrix{Float64}}, Nothing})
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\subjects\\subject.jl:62
 [3842] actor_on_next!(::BaseActorTrait{Message}, actor::Subject{Message, AsapScheduler, AsapScheduler}, data::Message{MvNormalMeanPrecision{Float64, Vector{Float64}, Matrix{Float64}}, Nothing})
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:250
 [3843] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:202 [inlined]
 [3844] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\subjects\\recent.jl:62 [inlined]
 [3845] actor_on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:250 [inlined]
 [3846] next!(actor::Rocket.RecentSubjectInstance{Message, Subject{Message, AsapScheduler, AsapScheduler}}, data::Message{MvNormalMeanPrecision{Float64, Vector{Float64}, Matrix{Float64}}, Nothing})
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:202
 [3847] next_received!(wrapper::Rocket.CollectLatestObservableWrapper{AbstractMessage, Rocket.RecentSubjectInstance{Message, Subject{Message, AsapScheduler, AsapScheduler}}, Vector{AbstractMessage}, BitVector, Vector{Teardown}, ReactiveMP.var\"#15#17\"{GenericProd, UnspecifiedFormConstraint}, typeof(ReactiveMP.reset_vstatus)}, data::DefferedMessage{Nothing, Tuple{ReactiveMP.MarginalObservable, ReactiveMP.MarginalObservable}, ReactiveMP.MessageMapping{MvNormalMeanPrecision, Val{:μ}, Marginalisation, Nothing, Val{(:out, :Λ)}, Nothing, Nothing, Nothing}}, index::CartesianIndex{1})
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\collected.jl:104
 [3848] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\collected.jl:93 [inlined]
 [3849] scheduled_next!(actor::Rocket.CollectLatestObservableInnerActor{AbstractMessage, CartesianIndex{1}, Rocket.CollectLatestObservableWrapper{AbstractMessage, Rocket.RecentSubjectInstance{Message, Subject{Message, AsapScheduler, AsapScheduler}}, Vector{AbstractMessage}, BitVector, Vector{Teardown}, ReactiveMP.var\"#15#17\"{GenericProd, UnspecifiedFormConstraint}, typeof(ReactiveMP.reset_vstatus)}}, value::DefferedMessage{Nothing, Tuple{ReactiveMP.MarginalObservable, ReactiveMP.MarginalObservable}, ReactiveMP.MessageMapping{MvNormalMeanPrecision, Val{:μ}, Marginalisation, Nothing, Val{(:out, :Λ)}, Nothing, Nothing, Nothing}}, ::AsapScheduler)
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\schedulers\\asap.jl:23
 [3850] on_next!(subject::Subject{AbstractMessage, AsapScheduler, AsapScheduler}, data::DefferedMessage{Nothing, Tuple{ReactiveMP.MarginalObservable, ReactiveMP.MarginalObservable}, ReactiveMP.MessageMapping{MvNormalMeanPrecision, Val{:μ}, Marginalisation, Nothing, Val{(:out, :Λ)}, Nothing, Nothing, Nothing}})
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\subjects\\subject.jl:62
 [3851] actor_on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:250 [inlined]
 [3852] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:202 [inlined]
 [3853] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\subjects\\recent.jl:62 [inlined]
 [3854] #7
      @ C:\\Users\\lange\\.julia\\packages\\RxInfer\\AszMu\\src\\rocket.jl:86 [inlined]
 [3855] limitstack(callback::RxInfer.var\"#7#8\"{Rocket.RecentSubjectInstance{AbstractMessage, Subject{AbstractMessage, AsapScheduler, AsapScheduler}}, DefferedMessage{Nothing, Tuple{ReactiveMP.MarginalObservable, ReactiveMP.MarginalObservable}, ReactiveMP.MessageMapping{MvNormalMeanPrecision, Val{:μ}, Marginalisation, Nothing, Val{(:out, :Λ)}, Nothing, Nothing, Nothing}}}, instance::RxInfer.LimitStackScheduler)
      @ RxInfer C:\\Users\\lange\\.julia\\packages\\RxInfer\\AszMu\\src\\rocket.jl:52
 [3856] scheduled_next!
      @ C:\\Users\\lange\\.julia\\packages\\RxInfer\\AszMu\\src\\rocket.jl:86 [inlined]
 [3857] actor_on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:255 [inlined]
 [3858] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:203 [inlined]
 [3859] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\scheduled.jl:12 [inlined]
 [3860] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:206 [inlined]
 [3861] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\operators\\map.jl:62 [inlined]
 [3862] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:206 [inlined]
 [3863] next_received!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined.jl:101 [inlined]
 [3864] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined.jl:68 [inlined]
 [3865] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:206 [inlined]
 [3866] next_received!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined_updates.jl:72 [inlined]
 [3867] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\observable\\combined_updates.jl:34 [inlined]
 [3868] scheduled_next!(actor::Rocket.CombineLatestUpdatesInnerActor{Marginal, Rocket.CombineLatestUpdatesActorWrapper{Tuple{ReactiveMP.MarginalObservable, ReactiveMP.MarginalObservable}, Rocket.CombineLatestInnerActor{Tuple{ReactiveMP.MarginalObservable, ReactiveMP.MarginalObservable}, Rocket.CombineLatestActorWrapper{Rocket.MStorage2{Nothing, Tuple{ReactiveMP.MarginalObservable, ReactiveMP.MarginalObservable}}, Rocket.MapActor{Tuple{Nothing, Tuple{ReactiveMP.MarginalObservable, ReactiveMP.MarginalObservable}}, Rocket.ScheduledActor{AbstractMessage, RxInfer.LimitStackScheduler, Rocket.RecentSubjectInstance{AbstractMessage, Subject{AbstractMessage, AsapScheduler, AsapScheduler}}}, ReactiveMP.var\"#178#181\"{ReactiveMP.MessageMapping{MvNormalMeanPrecision, Val{:μ}, Marginalisation, Nothing, Val{(:out, :Λ)}, Nothing, Nothing, Nothing}}}, PushNew, Rocket.UInt8UpdatesStatus}}, PushNew, Rocket.UInt8UpdatesStatus, typeof(identity), Nothing}}, value::Marginal{PointMass{Vector{Float64}}, Nothing}, ::AsapScheduler)
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\schedulers\\asap.jl:23
 [3869] on_next!(subject::Subject{Marginal, AsapScheduler, AsapScheduler}, data::Marginal{PointMass{Vector{Float64}}, Nothing})
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\subjects\\subject.jl:62
 [3870] actor_on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:250 [inlined]
 [3871] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:202 [inlined]
 [3872] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\subjects\\recent.jl:62 [inlined]
 [3873] actor_on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:250 [inlined]
 [3874] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:202 [inlined]
 [3875] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\operators\\map.jl:62 [inlined]
 [3876] scheduled_next!(actor::Rocket.MapActor{Message, Rocket.RecentSubjectInstance{Marginal, Subject{Marginal, AsapScheduler, AsapScheduler}}, typeof(as_marginal)}, value::Message{PointMass{Vector{Float64}}, Nothing}, ::AsapScheduler)
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\schedulers\\asap.jl:23
 [3877] on_next!(subject::Subject{Message, AsapScheduler, AsapScheduler}, data::Message{PointMass{Vector{Float64}}, Nothing})
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\subjects\\subject.jl:62
 [3878] actor_on_next!(::BaseActorTrait{Message}, actor::Subject{Message, AsapScheduler, AsapScheduler}, data::Message{PointMass{Vector{Float64}}, Nothing})
      @ Rocket C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:250
 [3879] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:202 [inlined]
 [3880] on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\subjects\\recent.jl:62 [inlined]
 [3881] actor_on_next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:250 [inlined]
 [3882] next!
      @ C:\\Users\\lange\\.julia\\packages\\Rocket\\sUaNu\\src\\actor.jl:202 [inlined]
 [3883] update!
      @ C:\\Users\\lange\\.julia\\packages\\ReactiveMP\\hqWuG\\src\\variables\\data.jl:79 [inlined]
 [3884] #133
      @ C:\\Users\\lange\\.julia\\packages\\ReactiveMP\\hqWuG\\src\\variables\\data.jl:87 [inlined]
 [3885] foreach(f::ReactiveMP.var\"#133#134\", itr::Base.Iterators.Zip{Tuple{Vector{DataVariable{Rocket.RecentSubjectInstance{Message, Subject{Message, AsapScheduler, AsapScheduler}}, ReactiveMP.MarginalObservable}}, Vector{Vector{Float64}}}})
      @ Base .\\abstractarray.jl:3097
 [3886] update!(datavars::Vector{DataVariable{Rocket.RecentSubjectInstance{Message, Subject{Message, AsapScheduler, AsapScheduler}}, ReactiveMP.MarginalObservable}}, data::Vector{Vector{Float64}})
      @ ReactiveMP C:\\Users\\lange\\.julia\\packages\\ReactiveMP\\hqWuG\\src\\variables\\data.jl:86
 [3887] batch_inference(; model::GraphPPL.ModelGenerator{typeof(LGDS), @Kwargs{param::Dict{Symbol, Any}, mΦhp::Vector{Float64}, mΦcv::Vector{Float64}, mΦs::Vector{Float64}, mTa::Vector{Float64}, Δt::Float64}, GraphPPL.PluginsCollection{Tuple{}}, RxInfer.ReactiveMPGraphPPLBackend}, data::@NamedTuple{x::Vector{Vector{Float64}}}, initialization::RxInfer.InitSpecification, constraints::GraphPPL.Constraints{Vector{GraphPPL.FactorizationConstraint}, Vector{GraphPPL.MarginalFormConstraint}, Vector{GraphPPL.MessageFormConstraint}, Dict{Function, GraphPPL.GeneralSubModelConstraints}, Dict{GraphPPL.FactorID, GraphPPL.SpecificSubModelConstraints}}, meta::GraphPPL.MetaSpecification, options::@NamedTuple{limit_stack_depth::Int64}, returnvars::Nothing, predictvars::Nothing, iterations::Int64, free_energy::Bool, free_energy_diagnostics::Tuple{RxInfer.ObjectiveDiagnosticCheckNaNs, RxInfer.ObjectiveDiagnosticCheckInfs}, showprogress::Bool, callbacks::Nothing, addons::Nothing, postprocess::DefaultPostprocess, warn::Bool, catch_exception::Bool)
      @ RxInfer C:\\Users\\lange\\.julia\\packages\\RxInfer\\AszMu\\src\\inference\\batch.jl:300
 [3888] batch_inference
      @ C:\\Users\\lange\\.julia\\packages\\RxInfer\\AszMu\\src\\inference\\batch.jl:94 [inlined]
 [3889] #infer#230
      @ C:\\Users\\lange\\.julia\\packages\\RxInfer\\AszMu\\src\\inference\\inference.jl:306 [inlined]"
}

The second sentence:

Use objective_diagnostics keyword argument in the inference function to suppress this error.

Confuses me a bit because this option objective_diagnostics doesn't seem to exist? At least I can't find it anywhere in the code or in the docs. Also I don't think 'surpressing' the error would solve the problem?

Sometimes I get it to make a few iterations, but it always dies after less than 5 iterations. When it does make an iteration the variances are so high the results are useless.

I would be really grateful if you could give me some pointers on how to tweak the model, the hyperparameters or how to circumvent the NaN computation :)

The code I am working on now can be found here:
full_model.zip

For convenience I also put the most important parts below:

Prior parameters

n_states = length(unknowns(sys))
n_measured = length(T_sample)

IdT = diageye(n_states)
IdM = diageye(n_measured)

# initial state
μ_T0 = [T_bias for _ in 1:n_states]
Σ_T0 = IdT * 1

# thermal capacitance prior
μ_C = [10 for _ in 1:n_states]
Σ_C = IdT * 1e2

# thermal resistance prior
μ_R = [10 for _ in 1:n_states]
Σ_R = IdT * 1e2

# area prior
n_areas = 2
μ_A = [10 for _ in 1:n_areas]
Σ_A = IdM * 1e2

# process noise
ν_p = n_states
Λ_p = IdT * 1e1

# measurement noise
ν_m = n_measured
Λ_m = IdM * 1e1

# buid prior params dict
prior_params = Dict{Symbol,Any}(
    :μ_T0 => μ_T0, 
    :Σ_T0 => Σ_T0, 
    :μ_C => μ_C, 
    :Σ_C => Σ_C, 
    :μ_R => μ_R, 
    :Σ_R => Σ_R, 
    :μ_A => μ_A, 
    :Σ_A => Σ_A, 
    :ν_p => ν_p, 
    :Λ_p => Λ_p, 
    :ν_m => ν_m, 
    :Λ_m => Λ_m
)

Transition function

function transition(T, lC, lR, lA, Φ_hp, Φ_cv, Φ_s, T_a, Δt)
    """
    Transition function for the thermal model

    Args:
        params: Dict{Symbol,Any} - prior parameters for the model   [-]
        T: Array{Float64,1} - state vector [T_i, T_e, T_h]          [K]
        lC: Array{Float64,1} - log capacitance parameter vector     [J/K]
        lR: Array{Float64,1} - log resistance parameter vector      [K/W]
        lA: Array{Float64,1} - log area parameter vector            [m^2]
        Φ_hp: Float64 - thermal power injected by heat pump         [W]
        Φ_cv: Float64 - thermal power injected by boiler            [W]
        Φ_s: Float64 - solar power                                  [W/m^2]
        T_a: Float64 - ambient temperature                          [K]
    """
    # unpack parameter vectors
    lCi, lCe, lCh = lC
    lRie, lRea, lRih = lR
    lAw, lAe = lA

    # define A matrix
    A = [
       -exp(-lCi - lRih) - exp(-lCi - lRie)  exp(-lCi - lRie)                     exp(-lCi - lRih);
        exp(-lCe - lRie)                     -exp(-lCe - lRea) - exp(-lCe - lRie) 0;
        exp(-lCh - lRih)                     0                                    -exp(-lCh - lRih)
    ]

    # define B matrix
    B = [
        0                 0         0          exp(lAw - lCi)    0;
        exp(-lCe - lRea)  0         0          0                 exp(lAe - lCe);
        0                 exp(-lCh) exp(-lCh)  0                 0
    ]

    # discrete transition matrices
    A_d = IdT + A * Δt
    B_d = B * Δt

    # input vector
    u = [T_a, Φ_hp, Φ_cv, Φ_s, Φ_s]

    # state transition to T_new
    return A_d * T + B_d * u
end

# emission matrix
C = [
    1 0 0;
    0 0 1
]

Model definition

@model function LGDS(param, x, mΦhp, mΦcv, mΦs, mTa, Δt)
    """
    Build the probabilistic parameter/state estimation model for the thermal model

    Args:
        x: Array{Float64,2} - [x_i, x_h] measured interior and heater temp. [K]
        mΦhp: Array{Float64,1} - measured heat pump power                   [W]
        mΦcv: Array{Float64,1} - measured boiler power                      [W]
        mΦs: Array{Float64,1} - measured solar irradiance                   [W/m^2]
        mTa: Array{Float64,1} - measured ambient temperature                [K]
        Δt: Float64 - time step between measurements                        [s]
    """
    # model parameter priors
    lC ~ MvNormalMeanCovariance(param[:μ_C], param[:Σ_C])
    lR ~ MvNormalMeanCovariance(param[:μ_R], param[:Σ_R])
    lA ~ MvNormalMeanCovariance(param[:μ_A], param[:Σ_A])

    # process noise
    W ~ Wishart(param[:ν_p], param[:Λ_p])

    # measurement noise
    V ~ Wishart(param[:ν_m], param[:Λ_m])

    # initialize state
    T_prev ~ MvNormalMeanCovariance(param[:μ_T0], param[:Σ_T0])
    
    # dynamics
    for i in 1:size(x)[1]
        
        # deterministic transition variable
        trans[i] := transition(T_prev, lC, lR, lA, mΦhp[i], mΦcv[i], mΦs[i], mTa[i], Δt)

        # Gaussian state
        T[i] ~ MvNormalMeanPrecision(trans[i], W)

        # Gaussian observations
        x[i] ~ MvNormalMeanPrecision(C * T[i], V)

        # update state
        T_prev = T[i]        
    end
end

# meta data such as approximation techniques
meta = @meta begin
    transition() -> UT()
end

# model initialization due to loops    
init = @initialization begin

    # initial messages
    μ(lC) = MvNormalMeanPrecision(zeros(n_states), 1e-3 * diageye(n_states))
    μ(lR) = MvNormalMeanPrecision(zeros(n_states), 1e-3 * diageye(n_states))
    μ(lA) = MvNormalMeanPrecision(zeros(n_measured), 1e-3 * diageye(2))

    # initial marginals
    q(W) = Wishart(n_states, diageye(n_states))
    q(V) = Wishart(n_measured, diageye(n_measured))
end

# model factorization constraints (for variational inference)
my_constraints = @constraints begin
    q(T, trans, W, V) = q(T, trans)q(W)q(V)
end

Inference function

dataset = (x = [collect(i) for i in zip(T_sample["C_i"], T_sample["C_h"])],);

n_iter = 100

result = infer(
    model = LGDS(param = prior_params, 
                mΦhp = U_sample[:Φ_hp],
                mΦcv = U_sample[:Φ_cv],
                mΦs = U_sample[:Φ_s],
                mTa = U_sample[:T_a],
                Δt = t_int),                
    meta = meta,
    iterations = n_iter,
    constraints = my_constraints,
    initialization = init,
    data = dataset,
    free_energy = true,
    options = (limit_stack_depth = 100,),
    showprogress = true
)

[HoangMHNguyen]

Hi @langestefan! Sorry for the late reply! Here is the code for online setting:

@model function myLGDS(x, mΦhp, mΦcv, mΦs, mTa, Δt, μ_C, Σ_C, μ_R, Σ_R, μ_A, Σ_A, ν_p, Λ_p, ν_m, Λ_m,μ_T,Σ_T)
    """
    Build the probabilistic parameter/state estimation model for the thermal model

    Args:
        x: Array{Float64,2} - [x_i, x_h] measured interior and heater temp. [K]
        mΦhp: Array{Float64,1} - measured heat pump power                   [W]
        mΦcv: Array{Float64,1} - measured boiler power                      [W]
        mΦs: Array{Float64,1} - measured solar irradiance                   [W/m^2]
        mTa: Array{Float64,1} - measured ambient temperature                [K]
        Δt: Float64 - time step between measurements                        [s]
    """
    # model parameter priors
    lC ~ MvNormalMeanCovariance(μ_C, Σ_C)
    lR ~ MvNormalMeanCovariance(μ_R, Σ_R)
    lA ~ MvNormalMeanCovariance(μ_A, Σ_A)

    # process noise
    W ~ Wishart(ν_p, Λ_p)

    # measurement noise
    V ~ Wishart(ν_m, Λ_m)

    # prev state
    T_prev ~ MvNormalMeanCovariance(μ_T, Σ_T)

    # Gaussian state
    T ~ MvNormalMeanPrecision(transition(T_prev, lC, lR, lA, mΦhp, mΦcv, mΦs, mTa, Δt), W)

    # Gaussian observations
    x ~ MvNormalMeanPrecision(C * T, V)
end

autoupdate = @autoupdates begin
    μ_C, Σ_C = mean_cov(q(lC))
    μ_R, Σ_R = mean_cov(q(lR))
    μ_A, Σ_A = mean_cov(q(lA))
    ν_p, Λ_p = RxInfer.params(q(W))
    ν_m, Λ_m = RxInfer.params(q(V))
    μ_T, Σ_T = mean_cov(q(T))
end
# meta data such as approximation techniques
meta = @meta begin
    transition() -> UT()
end

# model initialization due to loops    
filter_init = @initialization begin
    q(lC) = MvNormalMeanCovariance(prior_params[:μ_C], prior_params[:Σ_C])
    q(lR) = MvNormalMeanCovariance(prior_params[:μ_R], prior_params[:Σ_R])
    q(lA) = MvNormalMeanCovariance(prior_params[:μ_A], prior_params[:Σ_A])
    q(T) =  MvNormalMeanCovariance(prior_params[:μ_T0], prior_params[:Σ_T0])
    q(W) = WishartFast(prior_params[:ν_p], inv(prior_params[:Λ_p]))
    q(V) = WishartFast(prior_params[:ν_m], inv(prior_params[:Λ_m]))
end

# model factorization constraints (for variational inference)
my_constraints = @constraints begin
    q(T, T_prev, lC,lR,lA,W, V) = q(T, T_prev,lC,lR,lA)q(W)q(V)
end

n_iter = 10
filter_result = infer(
    model = myLGDS(Δt = t_int),  
    data = (x = dataset.x,mΦhp = U_sample[:Φ_hp],
        mΦcv = U_sample[:Φ_cv],
        mΦs = U_sample[:Φ_s],
        mTa = U_sample[:T_a],),              
    meta = meta,
    iterations = n_iter,
    constraints = my_constraints,
    initialization = filter_init,
    free_energy = true,
    keephistory    = length(dataset.x),
    historyvars    = (
        lC = KeepLast(),
        lR = KeepLast(),
        lA = KeepLast()
    ),
    autoupdates = autoupdate,
    options = (limit_stack_depth = 500,),
    showprogress = true,
    autostart      = true
)

In this setting, we basically do "filtering" update, and as @albertpod has mentioned, we don't have loops now. This setting will give you a converged FE. For the initialization, I use the parameters that we have specified for priors. Also, keep in mind that lR, lC and lA are in the log domain, so remember to set the prior values of these variables in the log domain too. For example, in your code, the prior of lA is μ_A=10, but this should be log(10) so that when you take the exponential, we get A = 10.

Hope this helps!

[langestefan]

Thanks, great work. I can't seem to use WishartFast, is that part of some specific package? This is what I have installed:

[31c24e10] Distributions v0.25.108
[62312e5e] ExponentialFamily v1.4.1
[86711068] RxInfer v3.2.0

Edit: had to explicitly call ExponentialFamily.WishartFast

[HoangMHNguyen]

you can import it from ExponentialFamily: import ExponentialFamily: WishartFast

[langestefan]

Yes, that worked. Seems like it's not exported by default. Could you perhaps explain how you came up with the constraints specification? I am still struggling to understand this part. For example if I specify q(T,T_prev,lC,lR,lA,W,V) = q(T, T_prev)q(lC,lR,lA)q(W)q(V) as the factorization constraint I get:

Unresolvable factorization constraint: Cannot resolve factorization constraint

As to the problem itself, the parameter estimates are still all over the place. I think I will code a few intermediate models of lower order to understand exactly where it starts breaking down, because it's a bit of a mystery to me now why this is happening.

[HoangMHNguyen]

First, we should note that we cannot use Variational Message Passing for a deterministic function because all variables in the function are deterministically related to each other. Thus, for deterministic functions, we should specify a joint distribution for all variables. The constraint q(T, T_prev,lC,lR,lA,W,V) = q(T, T_prev)q(lC,lR,lA)q(W)q(V) means you are imposing a mean-field assumption on the transition function, meaning you're gonna use VMP around this function which is not possible. My constraint specification simply puts a joint distribution over all variables that go into the transition function. Furthermore, my constraint implies a structured mean-field at the Normal node. Particularly, if you add an aux variable to the model as below

    # Gaussian state
    aux ~ transition(T_prev, lC, lR, lA, mΦhp, mΦcv, mΦs, mTa, Δt)
    T ~ MvNormalMeanPrecision(aux, W)

then my constraint q(T, T_prev, lC,lR,lA,W, V) = q(T,T_prev,lC,lR,lA)q(W)q(V) now is equivalent to q(T, aux,T_prev, lC,lR,lA,W, V) = q(T,aux)q(T_prev,lC,lR,lA)q(W)q(V), where q(T, aux) is the structured mean-field constraint.
Alternatively, you can have a naive mean-field constraint by specifying q(T, T_prev, lC,lR,lA,W, V) = q(T)q(T_prev,lC,lR,lA)q(W)q(V), which is equivalent to q(T, aux,T_prev, lC,lR,lA,W, V) = q(T)q(aux,T_prev,lC,lR,lA)q(W)q(V).

[albertpod]

Hi @langestefan,

Sorry for not getting back to you sooner. I've managed to redefine your model so that we can run coupled CTransition nodes within it. There are many moving parts in terms of hyperparameters, so the model needs rigorous analysis, which I believe you have a better grasp of than I do. I was primarily concerned with running the inference. Please take a look at the solution.

function transition_A(lCRA)
    lC = lCRA[1:n_states]
    lR = lCRA[n_states+1:n_states*2]
    lA = lCRA[n_states*2+1:end]

    lCi, lCe, lCh = lC
    lRie, lRea, lRih = lR
    lAw, lAe = lA

    A = [
       -exp(-lCi - lRih) - exp(-lCi - lRie)  exp(-lCi - lRie)                     exp(-lCi - lRih);
        exp(-lCe - lRie)                     -exp(-lCe - lRea) - exp(-lCe - lRie) 0;
        exp(-lCh - lRih)                     0                                    -exp(-lCh - lRih)
    ]

    Δt = 0.1
    A_d = IdT + A * Δt
end

function transition_B(lCRA)
    lC = lCRA[1:n_states]
    lR = lCRA[n_states+1:n_states*2]
    lA = lCRA[n_states*2+1:end]

    lCi, lCe, lCh = lC
    lRie, lRea, lRih = lR
    lAw, lAe = lA

    B = [
        0                 0         0          exp(lAw - lCi)    0;
        exp(-lCe - lRea)  0         0          0                 exp(lAe - lCe);
        0                 exp(-lCh) exp(-lCh)  0                 0
    ]


    Δt = 0.1
    B_d = B * Δt
end

@model function LGDS(param, x, mΦhp, mΦcv, mΦs, mTa, Δt)

    lCRA ~ MvNormalMeanPrecision([param[:μ_C]; param[:μ_R]; param[:μ_A]], diageye(length(param[:μ_C]) + length(param[:μ_R]) + length(param[:μ_A])) * 1e-2)

    W ~ Wishart(param[:ν_p], param[:Λ_p])

    V ~ Wishart(param[:ν_m], param[:Λ_m])

    T_prev ~ MvNormalMeanCovariance(param[:μ_T0], inv(param[:Σ_T0]))
    
    for i in 1:size(x)[1]
        
        u[i]  ~ MvNormalMeanPrecision([mTa[i], mΦhp[i], mΦcv[i], mΦs[i], mΦs[i]], 1e2diageye(5)) # we will pass it as RV as ContinuousTransition doesn't support (yet) passing data as input
        
        uB[i] ~ ContinuousTransition(u[i], lCRA, diageye(3))  where {meta = CTMeta(transition_B)}
        Ta[i] ~ ContinuousTransition(T_prev, lCRA, diageye(3)) where {meta = CTMeta(transition_A)}

        T[i] ~ MvNormalMeanPrecision(Ta[i] + uB[i], W)

        x[i] ~ MvNormalMeanPrecision(C * T[i], V)

        T_prev = T[i]        
    end
end

init = @initialization begin

    q(lCRA) = MvNormalMeanPrecision(zeros(2n_states + n_areas), diageye(2n_states + n_areas))
    q(W) = Wishart(n_states, diageye(n_states))
    q(V) = Wishart(n_measured, diageye(n_measured))
end

constraints = @constraints begin
    q(Ta, T_prev, u, uB, lCRA, W, V, T) = q(Ta, T_prev, T, u, uB)q(lCRA)q(V)q(W)
end

dataset = (x = [collect(i) for i in zip(T_sample["C_i"], T_sample["C_h"])],);


n_iter = 2

result = infer(
    model = LGDS(param = prior_params, 
                mΦhp = U_sample[:Φ_hp],
                mΦcv = U_sample[:Φ_cv],
                mΦs = U_sample[:Φ_s],
                mTa = U_sample[:T_a],
                Δt = t_int),                
    iterations = n_iter,
    constraints = constraints,
    initialization = init,
    data = dataset,
    free_energy = false,
    options = (limit_stack_depth = 100,),
    showprogress = true,
    returnvars = KeepLast()
)

There are a few things that must be mentioned here:

1. Apparently, CTransition doesn't handle nonlinearities inside the matrices well. This is due to the fact that CTransition uses linearization (via ForwardDiff) under the hood (no unscented), so you can immediately encounter explosive gradients (hence only 2 iterations work with the original transitions). If I remove exponentiation in transition_A and transition_B, the inference runs longer (up to 500 iterations), but that's not the correct physics of the system, so the inference results are off.
2. As you can see, we no longer have to initialize messages, only marginals, which is a good aspect, but your priors and initial marginals still matter.
3. There was something else, but I forgot (will add when/if I recall).

If you have questions regarding the model, initializations, or constraint specifications, please ask!

All in all, it's a cool problem. Good luck!

[langestefan]

@failin3 thanks, I tried something similar. For my problem the values of alpha and kappa that work are much smaller (in the range 1e-2). I tried a grid approach, but I'm still not able to get more than 5 or so iterations without running into NaN's.

Maybe there is something else I can try?

I also noticed something odd. I inspected the posterior for the state T, and it looks somewhat reasonable, as in close to the sampled values. But when I run the discrete model with a point estimate of the posterior (the MAP or mean) it looks nothing like the posterior. And that's strange because the discrete model is used to compute the mean of the next sample so it should match the posterior.

[albertpod]

@langestefan it will be interesting to try out the inference of our new method that is supposed to be released early next month.

There's also an opportunity to design your own node that will implement stable and robust approximation for the messages and products, but this ofc is more involved and I usually don't recommend this route to outside collaborators, but seems you @langestefan can try to go this path :)

[langestefan]

@albertpod sounds interesting, can you share anything about this new method?

I would like to try designing my own node, but I'm not sure I will have enough time for that within my project.

[albertpod]

The method will be based on this package https://github.com/ReactiveBayes/ExponentialFamilyProjection.jl. There'll be more docs once we integrate it to RxInfer.

[albertpod]

@langestefan perhaps you can also take a look at the use-case of DifferentialEquations.jl package with RxInfer outline here.