Visualizing Forney Factor Graphs when using RxInfer (Part 5)

Creating a separate FFGViz module

Bayesian Inference
Active Inference
TikZ
RxInfer
Julia
Author

Kobus Esterhuysen

Published

August 16, 2024

Modified

September 6, 2024

In this part (Part 5) we:

One of the distinguishing features of the RxInfer Julia package is that it uses reactive message passing (RMP). This means node updates in the factor graph happens asynchronously. A user usually tends to gloss over these details due to the increased complexity. When node updates happen in the more traditional way, i.e. synchronously (for example as with the PyMDP Python package), users tend to have the option of being aware of which node is currently being udpated. It is quite beneficial to conceptualize the workings of factor graphs in general and the visualization of the factor graph itself is often helpful. The existing visualization code is rather limited and that is the reason this project was undertaken. The ultimate goal is reach a visualization which conforms to the Constrained Forney Factor Graph (CFFG) as defined in

https://biaslab.github.io/pdf/arxiv/2306.08014.pdf

https://biaslab.github.io/pdf/arxiv/2306.02733.pdf

GraphPPL.jl exports the @model macro for model specification allowing the acceptance of a regular Julia function and then builds a Forney Factor Graph (FFG) under the hood. This graph is bipartite where nodes belong to one of two sets: one for factors and one for variables. Note that we use the terms node and link rather than vertex and edge to align with RxInfer terminology.

Node Taxonomy

  • nodes
    • facs (FFG/RxInfer factor nodes)
    • vars (FFG/RxInfer variable nodes)
      • parvars (parameter variables, e.g. \(\alpha, a, A_t, B, \gamma_t, \mu\))
      • sysvars (system variables, e.g. \(x_t, s_t, u_t\))

Node Visualization

  • a fac node is visualized by a
    • square (style fac) if it is a factor
    • small black square (style dlt) if it is a delta (\(\delta\)) factor
      • delta factors connect to their associated parvars and sysvars
  • a var node is visualized by a
    • point (style var) if connected to a single fac
    • equals-square (style eql) if connected to multiple facs
  • rather than converting the bipartite graph from RxInfer to a graph that only contains a single kind of node, we keep the bipartite graph and visualize the nodes differently

For each of the examples in this project the following steps happen:

  • A @model is created

  • The @model is conditioned on some data

    • ususally a few points so that the visualization does not become unwieldy

  • A RxInfer model is created (rxi_model)

  • A GraphPPL model is obtained (gppl_model)

  • A meta graph is obtained (meta_graph)

  • The graph is rendered by the existing implementation (GraphPlot.gplot())

  • Some vertex labels are inspected for the sake of interest

  • Setup global parameters (indicated by a leading underscore character) to influence the behavior of generate_tikz()

    • _san_node_ids is provided for node names that needs to be sanitized
    • _lup_enabled is set to true/false to enable lookup of math names
    • _lup_agc is set to true/false to enable appending of global_counter values
    • _raw_links_enabled is set to true/false to enable the drawing of raw links for reference
    • _CFFG is set to true/false to enable the drawing of beads to turn a FFG into a CFFG
    • _deltas_enabled is set to true/false to enable the drawing of delta decorations
    • _lup_dict is provided for all variable names needing math equivalents
    • _tsep specifies the time delta separation of nodes
    • _lup_dict specifies the lookup of canonical latex names for nodes
      • when _lup_agc is true, global_counter values are appended
      • when _lup_agc is false, global_counter values are not appended

  • The TikZ string is generated (generate_tikz())

    • heading of the FFG
    • GraphPPL model
    • layers for the FFG
    • layer names (if applicable)
      • cntrl (control layer)
      • paramB1 (parameter B1 layer)
      • paramB2 (parameter B2 layer)
      • stateU2 (state up 2 layer, i.e. 2 layers above state layer)
      • stateU1 (state up 1 layer, i.e. 1 layer above state layer)
      • state (state layer, i.e. where state transitions happen)
      • stateD1 (state down 1 layer, i.e. 1 layer below state layer)
      • stateD2 (state down 2 layer, i.e. 2 layers below state layer)
      • obser (observation layer)
      • prefr (preference layer, for setpoints)
      • paramA (parameter A layer)

  • The TikZ string is printed

  • The FFG is saved to a pdf file (show_tikz())

Setup notes

  • If additional Ubuntu packages are needed for a Julia package
    • sudo apt-get install package
  • If you want to run bash on the devcontainer in a local terminal; kobus@Kobuss-Mac-mini ~ %
    • docker exec -it b20ca0308c9b2698fd7c19fdac794941cd50ecc7a11e900b7e47e5018383e033 /bin/bash

Development cycle notes

  • execute & save the generated tikz/pgf code
  • copy the generated .tex file from the devcontainer to the local folder; kobus@Kobuss-Mac-mini ~ %
    • use vscode command line (View -> Terminal) to discover the folder name on the devcontainer
    • use the docker ui to copy the id of the container
    • identify the local folder name from where the .tex will be rendered
    • overlay the next command with what you discovered above
    • docker cp b20ca0308c9b2698fd7c19fdac794941cd50ecc7a11e900b7e47e5018383e033:/workspaces/2024-07-05^FFGViz2/myoutput.tex /Users/kobus/2024_LaTeX
  • in the /Users/kobus/2024_LaTeX project on the mini, the pdf will be automatically rendered when a change is detected in the copied .tex
versioninfo() ## Julia version
Julia Version 1.10.4
Commit 48d4fd48430 (2024-06-04 10:41 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 12 × Intel(R) Core(TM) i7-8700B CPU @ 3.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, skylake)
Threads: 1 default, 0 interactive, 1 GC (on 12 virtual cores)
Environment:
  JULIA_NUM_THREADS = 
import Pkg; Pkg.instantiate(); Pkg.activate()
Pkg.add(Pkg.PackageSpec(;name="RxInfer", version="3.6.0"))
Pkg.add(Pkg.PackageSpec(;name="MetaGraphsNext"))
Pkg.add(Pkg.PackageSpec(;name="GraphPlot"))
  Activating project at `~/.julia/environments/v1.10`
   Resolving package versions...
  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`
   Resolving package versions...
  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`
   Resolving package versions...
  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`
using RxInfer, MetaGraphsNext
using GraphPlot
include("/workspaces/2024-08-16^FFGViz5/FFGViz^v2.jl")
Layer
Pkg.status()
Status `~/.julia/environments/v1.10/Project.toml`
  [a2cc645c] GraphPlot v0.6.0
  [fa8bd995] MetaGraphsNext v0.7.0
  [86711068] RxInfer v3.6.0
##
# # If your PDF contains a TikZ picture (a vector graphic created using LaTeX), 
# # you can still use the method I mentioned earlier to display it in a Julia notebook. 
# # Here’s the method again for your reference:
# import Pkg; Pkg.add("PGFPlotsX")
# using PGFPlotsX

# # Create a struct that represents a PDF file
# struct PDF
#     file::String
# end

# # Define a function to show the PDF as inline SVG
# function Base.show(io::IO, ::MIME"image/svg+xml", pdf::PDF)
#     svg = first(splitext(pdf.file)) * ".svg"
#     PGFPlotsX.convert_pdf_to_svg(pdf.file, svg)
#     write(io, read(svg))
#     try; rm(svg; force=true); catch e; end
#     return nothing
# end

# # Now you can display a PDF file like this:
# display(PDF("myoutput.pdf"))

# # Internal Error: cairo context error: invalid matrix (not invertible)<0a>
# # cairo error: invalid matrix (not invertible)
_ytop = 30 #18 ## grid y-value at the top of picture
_show_grid = false ## displays a grid with the FFG if true
_node_size = "10mm" #"15mm" ## the length of the sides of factor node boxes (in millimeters)
## _show_var_points = false ## draws the var circles if true

_xsep = 2
# _xsep = 2.5

_tsep = 6*_xsep
# _tsep = 7*_xsep
# _tsep = 7.5*_xsep
# _tsep = 8*_xsep
# _tsep = 9*_xsep

_ysep = 3

_divn = setup_divn(15)
Dict{Any, Any} with 15 entries:
  5  => [0.166667, 0.333333, 0.5, 0.666667, 0.833333]
  12 => [0.0769231, 0.153846, 0.230769, 0.307692, 0.384615, 0.461538, 0.538462,…
  8  => [0.111111, 0.222222, 0.333333, 0.444444, 0.555556, 0.666667, 0.777778, …
  1  => [0.5]
  6  => [0.142857, 0.285714, 0.428571, 0.571429, 0.714286, 0.857143]
  11 => [0.0833333, 0.166667, 0.25, 0.333333, 0.416667, 0.5, 0.583333, 0.666667…
  9  => [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
  14 => [0.0666667, 0.133333, 0.2, 0.266667, 0.333333, 0.4, 0.466667, 0.533333,…
  3  => [0.25, 0.5, 0.75]
  7  => [0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875]
  4  => [0.2, 0.4, 0.6, 0.8]
  13 => [0.0714286, 0.142857, 0.214286, 0.285714, 0.357143, 0.428571, 0.5, 0.57…
  15 => [0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375, 0.5, 0.5625, 0.625…
  2  => [0.333333, 0.666667]
  10 => [0.0909091, 0.181818, 0.272727, 0.363636, 0.454545, 0.545455, 0.636364,…

Coin-toss Model

In this example, we are going to perform an exact inference for a coin-toss model that can be represented as:

\[\begin{aligned} p(\theta) &= \mathrm{Beta}(\theta|a, b),\\ p(y_i|\theta) &= \mathrm{Bern}(y_i|\theta),\\ \end{aligned}\]

where \(y_i \in \{0, 1\}\) is a binary observation induced by Bernoulli likelihood while \(\theta\) is a Beta prior distribution on the parameter of Bernoulli. We are interested in inferring the posterior distribution of \(\theta\).

The generative model is:

\[ \begin{aligned} p(\boldsymbol{y},\theta) &= p(\boldsymbol{y} \mid \theta) \cdot p(\theta) \\ &= \prod_{i=1}^N{p(y_i \mid \theta)} \cdot p(\theta) \\ &= \prod_{i=1}^N{\mathrm{Bern}(y_i \mid \theta)} \cdot \mathrm{Beta}(\theta \mid a,b) \end{aligned} \]

See Fig 4.1 in 2023_Bagaev_RPPfSBI_Thesis. So, \(N\) \(\delta\) factors, 1 Beta factor, \(N\) Ber factors.

@model function coin_model(y, a, b)
    ## We endow θ parameter of our model with some prior
    θ ~ Beta(a, b)

    ## We assume that outcome of each coin flip is governed by the Bernoulli distribution
    for i in eachindex(y)
        y[i] ~ Bernoulli(θ)
    end
end
coin_conditioned = coin_model(a = 2.0, b = 7.0) | (y = [ true, false, true ], )
coin_rxi_model = RxInfer.create_model(coin_conditioned)
coin_gppl_model = RxInfer.getmodel(coin_rxi_model)
coin_meta_graph = coin_gppl_model.graph
Meta graph based on a Graphs.SimpleGraphs.SimpleGraph{Int64} with vertex labels of type GraphPPL.NodeLabel, vertex metadata of type GraphPPL.NodeData, edge metadata of type GraphPPL.EdgeLabel, graph metadata given by GraphPPL.Context(0, coin_model, "", nothing, {Beta = 1, Bernoulli = 3}, {}, {(Bernoulli, 2) = Bernoulli_8, (Beta, 1) = Beta_4, (Bernoulli, 1) = Bernoulli_6, (Bernoulli, 3) = Bernoulli_10}, {:constvar_2 = constvar_2, :θ = θ_1, :constvar_3 = constvar_3}, {:y = ResizableArray{GraphPPL.NodeLabel,1}(GraphPPL.NodeLabel[y_5, y_7, y_9])}, {}, {}, Base.RefValue{Any}(nothing)), and default weight 1.0
##
## fieldnames(GraphPPL.Context)
# coin_gppl_model_context = GraphPPL.getcontext(coin_gppl_model)
# coin_gppl_model_context.individual_variables
# coin_gppl_model_context.vector_variables
# coin_gppl_model_context.factor_nodes
print_meta_graph(coin_meta_graph)
============== NODES ==============
1: θ_1
2: constvar_2
3: constvar_3
4: Beta_4
5: y_5
6: Bernoulli_6
7: y_7
8: Bernoulli_8
9: y_9
10: Bernoulli_10

============== LINKS ==============
θ_1 -> Beta_4
θ_1 -> Bernoulli_6
θ_1 -> Bernoulli_8
θ_1 -> Bernoulli_10
constvar_2 -> Beta_4
constvar_3 -> Beta_4
y_5 -> Bernoulli_6
y_7 -> Bernoulli_8
y_9 -> Bernoulli_10
meta_graph = coin_meta_graph

## Shorten some labels to make graph more readable
# str_labels = [string(lab) for lab in labels(meta_graph)]
# replacements = 
#     Pair("MvNormalMeanCovariance", "Nmc"), 
#     Pair("MvNormalMeanPrecision", "Nmp"), 
#     Pair("constvar", "cv"),
#     Pair("x", "XXXXX") ## make more obvious
# short_labels = [replace(s, replacements...) for s in str_labels]

GraphPlot.gplot( ## existing plotting functionality
    meta_graph,
    layout=spring_layout,
    nodelabel=collect(labels(meta_graph)),
    ## nodelabel=short_labels,
    nodelabelsize=0.1,
    NODESIZE=0.02, ## diameter of the nodes
    NODELABELSIZE=1.5,
    # nodelabelc="white",
    nodelabelc="green",
    nodelabeldist=0.0,
    nodefillc=nothing, ## "cyan"
    edgestrokec="red",
    ## ImageSize = (800, 800) ##- does not work
)
θ_1 constvar_2 constvar_3 Beta_4 y_5 Bernoulli_6 y_7 Bernoulli_8 y_9 Bernoulli_10 
_san_node_ids = Dict("dummy" => "dummy")
_lup_enabled = true
_lup_agc = false ## append global_counter values
_raw_links_enabled = false
_CFFG = true
_deltas_enabled = true
_tsep = 3*_xsep
_lup_dict = Dict{String, String}(
    "constvar_2_3" => _lup_agc ? "a\\_2\\_3" : "a",
    "constvar_4_5" => _lup_agc ? "b\\_4\\_5" : "b",
    "Beta_6" => _lup_agc ? "\\mathcal{B}eta\\_6" : "\\mathcal{B}eta",
    "Bernoulli_8" => _lup_agc ? "\\mathcal{B}er\\_6" : "\\mathcal{B}er",
    "Bernoulli_10" => _lup_agc ? "\\mathcal{B}er\\_10" : "\\mathcal{B}er",
    "Bernoulli_12" => _lup_agc ? "\\mathcal{B}er\\_12" : "\\mathcal{B}er",
    "θ_1" => _lup_agc ? "\\theta\\_1" : "\\theta",
    "y_7" => _lup_agc ? "y_1\\_7" : "y_1",
    "y_9" => _lup_agc ? "y_2\\_9" : "y_2",
    "y_11" => _lup_agc ? "y_3\\_11" : "y_3",
);
coin_kobus_graph = create_kobus_graph(coin_gppl_model)
1: θ_1
2: constvar_2
3: constvar_3
4: Beta_4
5: y_5
6: Bernoulli_6
7: y_7
8: Bernoulli_8
9: y_9
10: Bernoulli_10
Dict{String, Dict{Symbol, Any}} with 10 entries:
  "Bernoulli_10" => Dict(:node_lup=>"\\mathcal{B}er", :node_code=>10, :node_dat…
  "constvar_2"   => Dict(:node_lup=>"constvar\\_2", :node_code=>2, :node_data=>…
  "constvar_3"   => Dict(:node_lup=>"constvar\\_3", :node_code=>3, :node_data=>…
  "y_5"          => Dict(:node_lup=>"y\\_5", :node_code=>5, :node_data=>NodeDat…
  "Beta_4"       => Dict(:node_lup=>"Beta\\_4", :node_code=>4, :node_data=>Node…
  "Bernoulli_6"  => Dict(:node_lup=>"Bernoulli\\_6", :node_code=>6, :node_data=…
  "Bernoulli_8"  => Dict(:node_lup=>"\\mathcal{B}er", :node_code=>8, :node_data…
  "θ_1"          => Dict(:node_lup=>"\\theta", :node_code=>1, :node_data=>NodeD…
  "y_7"          => Dict(:node_lup=>"y_1", :node_code=>7, :node_data=>NodeData …
  "y_9"          => Dict(:node_lup=>"y_2", :node_code=>9, :node_data=>NodeData …
layers = Dict()
kobus_graph = coin_kobus_graph
gppl_model = coin_gppl_model

obser_sysvars, obser_facs, obser_parvars = find_obser_nodes(kobus_graph, "y", gppl_model)
layers["obser"] = Layer(
    name="obser",
    anchor=3*_xsep,
    iniparvars = [],
    inifac = "",
    parvar="",
    ## sysvars = ["y_7", "y_9", "y_11"],
    ## sysvars = str_vector_vars(gppl_model, "y"),
    sysvars = obser_sysvars,
    ## facs = [["Bernoulli_8", "Bernoulli_10", "Bernoulli_12"]],
    ## facs = [str_facs(gppl_model, "Bernoulli")],
    facs = [obser_facs],
    parvars = []
    ## parvars = [obser_parvars]
    )

## paramA_inifac = str_facs(gppl_model, "Beta")[1]
paramA_inifac = find_inifac_node(kobus_graph, "θ", gppl_model)
paramA_out_var, paramA_iniparvars = find_iniparvar_nodes(kobus_graph, paramA_inifac, gppl_model)    
layers["paramA"] = Layer(
    name="paramA",
    anchor=0*_xsep,
    ## iniparvars = ["constvar_2_3", "constvar_4_5"],
    iniparvars = paramA_iniparvars,
    ## inifac = "Beta_6",
    inifac = paramA_inifac,
    ## parvar="θ_1",
    parvar=paramA_out_var,
    sysvars = [],
    facs = [],
    parvars = [],
    link_fac_side = "west")
## layers["paramA"].links = layers["obser"].facs    

layers["paramA"].links = layers["obser"].facs
obser_var_label_strs: ["y_5", "y_7", "y_9"]
obser_var_label_strs[i]: y_5
    var_neighbours[j]: Bernoulli_6
        in_neighbor_label: θ_1
        out_neighbor_label: y_5
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(y_5, out, NodeData in context  with properties name = y, index = 1), (θ_1, p, NodeData in context  with properties name = θ, index = nothing)]
obser_var_label_strs[i]: y_7
    var_neighbours[j]: Bernoulli_8
        in_neighbor_label: θ_1
        out_neighbor_label: y_7
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(y_7, out, NodeData in context  with properties name = y, index = 2), (θ_1, p, NodeData in context  with properties name = θ, index = nothing)]
obser_var_label_strs[i]: y_9
    var_neighbours[j]: Bernoulli_10
        in_neighbor_label: θ_1
        out_neighbor_label: y_9
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(y_9, out, NodeData in context  with properties name = y, index = 3), (θ_1, p, NodeData in context  with properties name = θ, index = nothing)]
    var_neighbours[1]: Beta_4
        in_neighbor_label: constvar_2
        out_neighbor_label: θ_1
inifac_out_neighbor_str: θ_1
inifac_str: Beta_4
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(θ_1, out, NodeData in context  with properties name = θ, index = nothing), (constvar_2, a, NodeData in context  with properties name = constvar, index = nothing), (constvar_3, b, NodeData in context  with properties name = constvar, index = nothing)]
        out_var: θ_1
1-element Vector{Vector{String}}:
 ["Bernoulli_6", "Bernoulli_8", "Bernoulli_10"]
coin_tikz = generate_tikz(
    heading = "Coin Toss",
    Model = coin_gppl_model,
    layers = layers) ## fac layers, i.e. not var (vars are handled with assoced layer)
## print(coin_tikz)
show_tikz(coin_tikz); ## write to file rather than show in notebook
\documentclass{standalone}
\usepackage{tikz}
\usetikzlibrary{calc}
\usepackage{amsmath}
\usepackage{bm}
\newcommand\ns{10mm} %% define node (minimum) size
\newcommand\nr{.5*\ns} %% define node radius
\newcommand\ar{.4*\ns} %% define arc radius
\newcommand\br{.19*\ns} %% define bead radius
\DeclareUnicodeCharacter{03B8}{\theta}
\DeclareUnicodeCharacter{03B3}{\gamma}
\DeclareUnicodeCharacter{209C}{\0}
\DeclareUnicodeCharacter{208B}{\0}
\DeclareUnicodeCharacter{2081}{\0}
\begin{document}
\begin{tikzpicture}[
    scale=1, 
    draw=black, 
    fac/.style={rectangle, draw=black!100, fill=black!0, minimum size=\ns, inner sep=0pt},
    dlt/.style={rectangle, draw=black!100, fill=black!100, minimum size=.2*\ns},
    var/.style={circle, draw=gray!100, fill=gray!100, minimum size=1mm, inner sep=0pt},
    bus/.style={rectangle, draw=black!100, fill=black!100, minimum width=30mm, minimum height=.5mm, inner sep=0pt},
    eql/.style={rectangle, draw=black!100, fill=black!0, minimum size=.5*\ns},
    text=black,
    ]
\node[font=\Large] at (4,29) {\textbf{Coin Toss}};
⋮
⋮
⋮
\path[use as bounding box] ([xshift=-1cm, yshift=-1cm]current bounding box.south west) rectangle ([xshift=1cm, yshift=2cm]current bounding box.north east);
\end{tikzpicture}
\end{document}

Coin Toss

Hidden Markov Model with Control (Part 1)

There are 1-hot encodings on all random variables. Using regular bold font to indicate the 1-hot encodings, and also setting \(t=1\), we have:

\[\begin{aligned} p(\boldsymbol{x}, \boldsymbol{s}) &\propto \underbrace{p(\mathbf s_0)}_{\substack{\text{Initial} \\ \text{state} \\ \text{prior}}} \prod_{k=1}^{T+1} \underbrace{p(\mathbf x_k \mid \mathbf s_k)}_{\substack{\text{Observation} \\ \text{model}}} \, \underbrace{p(\mathbf s_k \mid \mathbf s_{k-1}}_{\substack{\text{Transition} \\ \text{model}}}) \end{aligned}\]

where

\[\begin{aligned} p(\mathbf s_{0}) &= \mathcal{Cat}(\mathbf s_{0} ∣ \mathbf d) \\ p(\mathbf x_k ∣ \mathbf s_k) &= \mathcal{Cat}(\mathbf x_k ∣ \mathbf A \mathbf s_k) \\ p(\mathbf s_k ∣ \mathbf s_{k-1}) &= \mathcal{Cat}(\mathbf s_k ∣ \mathbf B \mathbf s_{k-1}) \\ \end{aligned}\]

@model function hidden_markov_model(x)
    B ~ MatrixDirichlet(ones(3, 3))
    A ~ MatrixDirichlet([10.0 1.0 1.0; 
                         1.0 10.0 1.0; 
                         1.0 1.0 10.0 ])
    d = fill(1.0/3.0, 3)
    s₀ ~ Categorical(d)
    
    sₖ₋₁ = s₀
    for k in eachindex(x)
        s[k] ~ Transition(sₖ₋₁, B)
        x[k] ~ Transition(s[k], A)
        sₖ₋₁ = s[k]
    end
end
hmm_conditioned = hidden_markov_model() | (x = [[1.0, 0.0, 0.0], [0.0, 0.0, 1.0]],)
hmm_rxi_model = RxInfer.create_model(hmm_conditioned)
hmm_gppl_model = RxInfer.getmodel(hmm_rxi_model)
hmm_meta_graph = hmm_gppl_model.graph
Meta graph based on a Graphs.SimpleGraphs.SimpleGraph{Int64} with vertex labels of type GraphPPL.NodeLabel, vertex metadata of type GraphPPL.NodeData, edge metadata of type GraphPPL.EdgeLabel, graph metadata given by GraphPPL.Context(0, hidden_markov_model, "", nothing, {MatrixDirichlet = 2, Categorical{P} where P<:Real = 1, Transition = 4}, {}, {(Transition, 4) = Transition_17, (MatrixDirichlet, 1) = MatrixDirichlet_3, (Transition, 3) = Transition_15, (Transition, 1) = Transition_11, (MatrixDirichlet, 2) = MatrixDirichlet_6, (Transition, 2) = Transition_13, (Categorical{P} where P<:Real, 1) = Categorical{P} where P<:Real_9}, {:A = A_4, :constvar_5 = constvar_5, :B = B_1, :s₀ = s₀_7, :constvar_8 = constvar_8, :constvar_2 = constvar_2}, {:s = ResizableArray{GraphPPL.NodeLabel,1}(GraphPPL.NodeLabel[s_10, s_14]), :x = ResizableArray{GraphPPL.NodeLabel,1}(GraphPPL.NodeLabel[x_12, x_16])}, {}, {}, Base.RefValue{Any}(nothing)), and default weight 1.0
##
## fieldnames(GraphPPL.Context)
# hmm_gppl_model_context = GraphPPL.getcontext(hmm_gppl_model)
# hmm_gppl_model_context.individual_variables
# hmm_gppl_model_context.vector_variables
# hmm_gppl_model_context.factor_nodes
print_meta_graph(hmm_meta_graph)
============== NODES ==============
1: B_1
2: constvar_2
3: MatrixDirichlet_3
4: A_4
5: constvar_5
6: MatrixDirichlet_6
7: s₀_7
8: constvar_8
9: Categorical{P} where P<:Real_9
10: s_10
11: Transition_11
12: x_12
13: Transition_13
14: s_14
15: Transition_15
16: x_16
17: Transition_17

============== LINKS ==============
B_1 -> MatrixDirichlet_3
B_1 -> Transition_11
B_1 -> Transition_15
constvar_2 -> MatrixDirichlet_3
A_4 -> MatrixDirichlet_6
A_4 -> Transition_13
A_4 -> Transition_17
constvar_5 -> MatrixDirichlet_6
s₀_7 -> Categorical{P} where P<:Real_9
s₀_7 -> Transition_11
constvar_8 -> Categorical{P} where P<:Real_9
s_10 -> Transition_11
s_10 -> Transition_13
s_10 -> Transition_15
x_12 -> Transition_13
s_14 -> Transition_15
s_14 -> Transition_17
x_16 -> Transition_17
meta_graph = hmm_meta_graph

## Shorten some labels to make graph more readable
# str_labels = [string(lab) for lab in labels(meta_graph)]
# replacements = 
#     Pair("MvNormalMeanCovariance", "Nmc"), 
#     Pair("MvNormalMeanPrecision", "Nmp"), 
#     Pair("constvar", "cv"),
#     Pair("x", "XXXXX") ## make more obvious
# short_labels = [replace(s, replacements...) for s in str_labels]

GraphPlot.gplot( ## existing plotting functionality
    meta_graph,
    layout=spring_layout,
    nodelabel=collect(labels(meta_graph)),
    ## nodelabel=short_labels,
    nodelabelsize=0.1,
    NODESIZE=0.02, ## diameter of the nodes
    NODELABELSIZE=1.5,
    # nodelabelc="white",
    nodelabelc="green",
    nodelabeldist=0.0,
    nodefillc=nothing, ## "cyan"
    edgestrokec="red",
    ## ImageSize = (800, 800) ##- does not work
)
B_1 constvar_2 MatrixDirichlet_3 A_4 constvar_5 MatrixDirichlet_6 s₀_7 constvar_8 Categorical{P} where P<:Real_9 s_10 Transition_11 x_12 Transition_13 s_14 Transition_15 x_16 Transition_17 
## _san_node_ids = Dict("dummy" => "dummy")
_san_node_ids = Dict("Categorical{P} where P<:Real_12" => "Categorical_12")
_lup_enabled = true
_lup_agc = false ## append global_counter values
_raw_links_enabled = false
_CFFG = true
_deltas_enabled = true
_tsep = 3*_xsep
_lup_dict = Dict{String, String}(
    "MatrixDirichlet_4" => _lup_agc ? "\\mathcal{M}Dir\\_4" : "\\mathcal{M}Dir",
    "B_1" => _lup_agc ? "\\mathbf{B}\\_1" : "\\mathbf{B}",
    "constvar_10_11" => _lup_agc ? "\\mathbf{d}\\_10\\_11" : "\\mathbf{d}",
    "Categorical_12" => _lup_agc ? "\\mathcal{C}at\\_12" : "\\mathcal{C}at",
    "Transition_14" => _lup_agc ? "\\mathcal{T}\\_14" : "\\mathcal{T}",
    "Transition_18" => _lup_agc ? "\\mathcal{T}\\_18" : "\\mathcal{T}",
    "Transition_16" => _lup_agc ? "\\mathcal{T}\\_16" : "\\mathcal{T}",
    "Transition_20" => _lup_agc ? "\\mathcal{T}\\_20" : "\\mathcal{T}",
    "MatrixDirichlet_8" => _lup_agc ? "\\mathcal{M}Dir\\_8" : "\\mathcal{M}Dir",
    "A_5" => _lup_agc ? "\\mathbf{A}\\_5" : "\\mathbf{A}",
    "constvar_2_3" => _lup_agc ? "\\mathbf{B_0}\\_2\\_3" : "\\mathbf{B_0}",
    "constvar_6_7" => _lup_agc ? "\\mathbf{A_0}\\_6\\_7" : "\\mathbf{A_0}",
    ## "s₀_9" => "s_0", ## s₀ must always be in _lup_dict
    "s₀_9" => _lup_agc ? "\\mathbf{s}_0\\_9" : "\\mathbf{s}_0",
    "s_13" => _lup_agc ? "\\mathbf{s}_1\\_13" : "\\mathbf{s}_1",
    "s_17" => _lup_agc ? "\\mathbf{s}_2\\_17" : "\\mathbf{s}_2",
    "x_15" => _lup_agc ? "\\mathbf{x}_1\\_15" : "\\mathbf{x}_1",
    "x_19" => _lup_agc ? "\\mathbf{x}_2\\_19" : "\\mathbf{x}_2",
);
hmm_kobus_graph = create_kobus_graph(hmm_gppl_model)
1: B_1
2: constvar_2
3: MatrixDirichlet_3
4: A_4
5: constvar_5
6: MatrixDirichlet_6
7: s₀_7
8: constvar_8
9: Categorical{P} where P<:Real_9
10: s_10
11: Transition_11
12: x_12
13: Transition_13
14: s_14
15: Transition_15
16: x_16
17: Transition_17
Dict{String, Dict{Symbol, Any}} with 17 entries:
  "x_12"                           => Dict(:node_lup=>"x\\_12", :node_code=>12,…
  "MatrixDirichlet_6"              => Dict(:node_lup=>"MatrixDirichlet\\_6", :n…
  "B_1"                            => Dict(:node_lup=>"\\mathbf{B}", :node_code…
  "Transition_17"                  => Dict(:node_lup=>"Transition\\_17", :node_…
  "Transition_15"                  => Dict(:node_lup=>"Transition\\_15", :node_…
  "constvar_2"                     => Dict(:node_lup=>"constvar\\_2", :node_cod…
  "Categorical{P} where P<:Real_9" => Dict(:node_lup=>"Categorical{P} where P<:…
  "Transition_11"                  => Dict(:node_lup=>"Transition\\_11", :node_…
  "constvar_8"                     => Dict(:node_lup=>"constvar\\_8", :node_cod…
  "MatrixDirichlet_3"              => Dict(:node_lup=>"MatrixDirichlet\\_3", :n…
  "s₀_7"                           => Dict(:node_lup=>"s₀\\_7", :node_code=>7, …
  "constvar_5"                     => Dict(:node_lup=>"constvar\\_5", :node_cod…
  "s_10"                           => Dict(:node_lup=>"s\\_10", :node_code=>10,…
  "Transition_13"                  => Dict(:node_lup=>"Transition\\_13", :node_…
  "x_16"                           => Dict(:node_lup=>"x\\_16", :node_code=>16,…
  "A_4"                            => Dict(:node_lup=>"A\\_4", :node_code=>4, :…
  "s_14"                           => Dict(:node_lup=>"s\\_14", :node_code=>14,…
layers = Dict()
kobus_graph = hmm_kobus_graph
gppl_model = hmm_gppl_model

## paramB1_inifac = "MatrixDirichlet_4"
paramB1_inifac = find_inifac_node(kobus_graph, "B", gppl_model)
paramB1_out_var, paramB1_iniparvars = find_iniparvar_nodes(kobus_graph, paramB1_inifac, gppl_model)
layers["paramB1"] = Layer(
    name="paramB1",
    anchor=1*_xsep,
    ## iniparvars=["constvar_2_3"],
    iniparvars=paramB1_iniparvars,
    ## inifac="MatrixDirichlet_4",
    inifac=paramB1_inifac,
    ## parvar="B_1",
    parvar=paramB1_out_var,
    sysvars = [],
    facs = [],
    parvars = [],
    link_fac_side="north", link_fac_offset=_divn[1][1])
## layers["paramB1"].links = layers["state0"].facs

state_inifac = find_inifac_node(kobus_graph, "s₀", gppl_model)
_, state_vars, state_facs = find_state_nodes(kobus_graph, state_inifac, gppl_model)
state_out_var, state_iniparvars = find_iniparvar_nodes(kobus_graph, state_inifac, gppl_model)
layers["state"] = Layer(
    name="state",
    anchor=3*_xsep,
    ## iniparvars = ["constvar_10_11"],
    iniparvars = state_iniparvars,
    inifac = state_inifac,
    parvar="",
    ## sysvars = ["s₀_9", "s_13", "s_17"],
    sysvars = state_vars, ## contains only sysvars
    ## facs = [["Transition_14", "Transition_18"]],
    facs = [state_facs],
    parvars = [])
## layers["state"].links = layers["obser"].facs

obser_sysvars, obser_facs, obser_parvars = find_obser_nodes(kobus_graph, "x", gppl_model)
layers["obser"] = Layer(
    name="obser",
    anchor=5*_xsep,
    iniparvars = [],
    inifac = "",
    parvar="",
    ## sysvars = ["x_15", "x_19"],
    sysvars = obser_sysvars,
    ## facs = [["Transition_16", "Transition_20"]],
    facs = [obser_facs],
    parvars = []
    ## parvars = [obser_parvars]
    )

paramA_inifac = find_inifac_node(kobus_graph, "A", gppl_model)
paramA_out_var, paramA_iniparvars = find_iniparvar_nodes(kobus_graph, paramA_inifac, gppl_model)
layers["paramA"] = Layer(
    name="paramA",
    anchor=1*_xsep,
    ## iniparvars=["constvar_6_7"],
    iniparvars=paramA_iniparvars,
    ## inifac="MatrixDirichlet_8",
    inifac=paramA_inifac,
    ## parvar="A_5",
    parvar=paramA_out_var,
    sysvars = [],
    facs = [],
    parvars = [],
    link_fac_side="west")
## layers["paramA"].links = layers["obser"].facs

layers["paramB1"].links = layers["state"].facs
layers["state"].links = layers["obser"].facs
layers["paramA"].links = layers["obser"].facs
    var_neighbours[1]: MatrixDirichlet_3
        in_neighbor_label: constvar_2
        out_neighbor_label: B_1
inifac_out_neighbor_str: B_1
inifac_str: MatrixDirichlet_3
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(B_1, out, NodeData in context  with properties name = B, index = nothing), (constvar_2, a, NodeData in context  with properties name = constvar, index = nothing)]
        out_var: B_1
    var_neighbours[1]: Categorical{P} where P<:Real_9
        in_neighbor_label: constvar_8
        out_neighbor_label: s₀_7
inifac_out_neighbor_str: s₀_7
inifac_str: Categorical{P} where P<:Real_9
var: (s₀_7, out, NodeData in context  with properties name = s₀, index = nothing)
=== NEXT var node is: s₀_7
var_label: s₀_7
var_index: 0
var_neighbours: Any["Categorical{P} where P<:Real_9", "Transition_11"]
    var_neighbours[1]: Categorical{P} where P<:Real_9
        in_neighbor_label: constvar_8
        out_neighbor_label: s₀_7
    var_neighbours[2]: Transition_11
        in_neighbor_label: s₀_7
        out_neighbor_label: s_10
        === NEXT fac node is: Transition_11
var: (s_10, out, NodeData in context  with properties name = s, index = 1)
=== NEXT var node is: s_10
var_label: s_10
var_index: 1
var_neighbours: Any["Transition_11", "Transition_13", "Transition_15"]
    var_neighbours[1]: Transition_11
        in_neighbor_label: s₀_7
        out_neighbor_label: s_10
    var_neighbours[2]: Transition_13
        in_neighbor_label: s_10
        out_neighbor_label: x_12
    var_neighbours[3]: Transition_15
        in_neighbor_label: s_10
        out_neighbor_label: s_14
        === NEXT fac node is: Transition_15
var: (s_14, out, NodeData in context  with properties name = s, index = 2)
=== NEXT var node is: s_14
var_label: s_14
var_index: 2
var_neighbours: Any["Transition_15", "Transition_17"]
    var_neighbours[1]: Transition_15
        in_neighbor_label: s_10
        out_neighbor_label: s_14
    var_neighbours[2]: Transition_17
        in_neighbor_label: s_14
        out_neighbor_label: x_16
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(s₀_7, out, NodeData in context  with properties name = s₀, index = nothing), (constvar_8, p, NodeData in context  with properties name = constvar, index = nothing)]
        out_var: s₀_7
obser_var_label_strs: ["x_12", "x_16"]
obser_var_label_strs[i]: x_12
    var_neighbours[j]: Transition_13
        in_neighbor_label: s_10
        out_neighbor_label: x_12
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(x_12, out, NodeData in context  with properties name = x, index = 1), (s_10, in, NodeData in context  with properties name = s, index = 1), (A_4, a, NodeData in context  with properties name = A, index = nothing)]
obser_var_label_strs[i]: x_16
    var_neighbours[j]: Transition_17
        in_neighbor_label: s_14
        out_neighbor_label: x_16
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(x_16, out, NodeData in context  with properties name = x, index = 2), (s_14, in, NodeData in context  with properties name = s, index = 2), (A_4, a, NodeData in context  with properties name = A, index = nothing)]
    var_neighbours[1]: MatrixDirichlet_6
        in_neighbor_label: constvar_5
        out_neighbor_label: A_4
inifac_out_neighbor_str: A_4
inifac_str: MatrixDirichlet_6
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(A_4, out, NodeData in context  with properties name = A, index = nothing), (constvar_5, a, NodeData in context  with properties name = constvar, index = nothing)]
        out_var: A_4
1-element Vector{Vector{String}}:
 ["Transition_13", "Transition_17"]
hmm_tikz = generate_tikz(
    heading = "Hidden Markov Model",
    Model = hmm_gppl_model,
    layers = layers) ## fac layers, i.e. not var (vars are handled with assoced layer)
## print(hmm_tikz)
show_tikz(hmm_tikz); ## write to file rather than show in notebook
\documentclass{standalone}
\usepackage{tikz}
\usetikzlibrary{calc}
\usepackage{amsmath}
\usepackage{bm}
\newcommand\ns{10mm} %% define node (minimum) size
\newcommand\nr{.5*\ns} %% define node radius
\newcommand\ar{.4*\ns} %% define arc radius
\newcommand\br{.19*\ns} %% define bead radius
\DeclareUnicodeCharacter{03B8}{\theta}
\DeclareUnicodeCharacter{03B3}{\gamma}
\DeclareUnicodeCharacter{209C}{\0}
\DeclareUnicodeCharacter{208B}{\0}
\DeclareUnicodeCharacter{2081}{\0}
\begin{document}
\begin{tikzpicture}[
    scale=1, 
    draw=black, 
    fac/.style={rectangle, draw=black!100, fill=black!0, minimum size=\ns, inner sep=0pt},
    dlt/.style={rectangle, draw=black!100, fill=black!100, minimum size=.2*\ns},
    var/.style={circle, draw=gray!100, fill=gray!100, minimum size=1mm, inner sep=0pt},
    bus/.style={rectangle, draw=black!100, fill=black!100, minimum width=30mm, minimum height=.5mm, inner sep=0pt},
    eql/.style={rectangle, draw=black!100, fill=black!0, minimum size=.5*\ns},
    text=black,
    ]
\node[font=\Large] at (4,29) {\textbf{Hidden Markov Model}};
⋮
⋮
⋮
\path[use as bounding box] ([xshift=-1cm, yshift=-1cm]current bounding box.south west) rectangle ([xshift=1cm, yshift=2cm]current bounding box.north east);
\end{tikzpicture}
\end{document}
seq: ["Transition_11"]
... doing the else ...
seq: ["Transition_15"]
... doing the else ...

Hidden Markov Model

Sales Forecasting with Time-Varying Autoregressive Models (Part 4)

Our model for the Variational Bayesian Inference of a latent autoregressive model can be represented mathematically as follows:

\[\begin{aligned} p(\gamma) &= \mathcal{Gam}(\gamma|a, b)\\ p(\boldsymbol{\theta}) &= \mathcal{N}(\boldsymbol{\theta}|\mathbf{m}_{\theta}, \mathbf{V}_{\theta})\\ p(\mathbf{s}_t \mid \mathbf{s}_{t-1}, \boldsymbol{\theta}) &= \mathcal{N}\left(\mathbf{s}_t \mid \mathbf{B}\left(\boldsymbol{\theta}\right) \mathbf{s}_{t-1},\, \mathbf{V}(\gamma)\right)\\ p(x_t \mid \mathbf{s}_t) &= \mathcal{N}(x_t \mid \mathbf{s}_t, \tau^{-1}) \end{aligned}\]

where

  • \(M^{AR}\) is the order of the autoregression process

  • \(\boldsymbol{\theta}=(\theta_1, \theta_2, ..., \theta_{M^{AR}})^{\top}\) is a vector of transition coefficients

  • \(\mathbf{s}_{t-1} = (s_{t-1}, s_{t-2}, ..., s_{t-M^{AR}})^\top\) is the previous state

  • \(\mathbf{s}_t = (s_{t}, s_{t-1}, ..., s_{t-M^{AR}+1})^\top\) is the current state

  • \(\gamma\) is the transition precision

  • \(\tau\) is the observation precision

  • \(\mathbf{v}_u = (1,0,...,0)^\top\) is an \(M^{AR}\)-dimensional unit vector selecting the first component of the vector that follows

  • \(\mathbf{V}(\gamma) = (1/ \gamma)\mathbf{v}_u\mathbf{v}_u^\top\) is the transition covariance

  • \(\mathbf{B}(\boldsymbol{\theta}) = \begin{bmatrix} \boldsymbol{\qquad\theta}^\top \\ \mathbf{I}_{M^{AR}-1} & \mathbf{0} \end{bmatrix}\) is the transition matrix

The following diagram represents an autoregressive process:

Autoregressive process
@model function lar_model(
    x, ##. data/observations 
    𝚃ᴬᴿ, ##. Uni/Multi variate 
    Mᴬᴿ, ##. AR order
    vᵤ, ##. unit vector
    τ) ##. observation precision     
        ## Priors
        γ  ~ Gamma(α = 1.0, β = 1.0) ##. for transition precision    
        if 𝚃ᴬᴿ === Multivariate
            θ  ~ MvNormal(μ = zeros(Mᴬᴿ), Λ = diageye(Mᴬᴿ)) ##.kw μ,Λ only work inside macro
            s₀ ~ MvNormal(μ = zeros(Mᴬᴿ), Λ = diageye(Mᴬᴿ)) ##.kw μ,Λ only work inside macro
        else ## Univariate
            θ  ~ Normal(μ = 0.0, γ = 1.0)
            s₀ ~ Normal(μ = 0.0, γ = 1.0)
        end
        sₜ₋₁ = s₀
        for t in eachindex(x)
            s[t] ~ AR(sₜ₋₁, θ, γ) #.Eq (2b)
            if 𝚃ᴬᴿ === Multivariate
                x[t] ~ Normal(μ = dot(vᵤ, s[t]), γ = τ) #.Eq (2c)
            else
                x[t] ~ Normal(μ = vᵤ*s[t], γ = τ) #.Eq (2c)
            end
            sₜ₋₁ = s[t]
        end
    end
𝚃ᴬᴿ = Univariate
m = 1
τ̃ = 0.001 ## assumed observation precision
lar_conditioned = lar_model(
    𝚃ᴬᴿ=𝚃ᴬᴿ, 
    Mᴬᴿ=m, 
    vᵤ=ReactiveMP.ar_unit(𝚃ᴬᴿ, m), 
    τ=τ̃
) | (x = [266.0, 145.0, 183.0],)
lar_model(𝚃ᴬᴿ = Univariate, Mᴬᴿ = 1, vᵤ = 1.0, τ = 0.001) conditioned on: 
  x = [266.0, 145.0, 183.0]
lar_rxi_model = RxInfer.create_model(lar_conditioned)
lar_gppl_model = RxInfer.getmodel(lar_rxi_model)
lar_meta_graph = lar_gppl_model.graph
Meta graph based on a Graphs.SimpleGraphs.SimpleGraph{Int64} with vertex labels of type GraphPPL.NodeLabel, vertex metadata of type GraphPPL.NodeData, edge metadata of type GraphPPL.EdgeLabel, graph metadata given by GraphPPL.Context(0, lar_model, "", nothing, {* = 3, GammaShapeRate = 1, AR = 3, NormalMeanPrecision = 5}, {}, {(GammaShapeRate, 1) = GammaShapeRate_4, (NormalMeanPrecision, 3) = NormalMeanPrecision_20, (NormalMeanPrecision, 2) = NormalMeanPrecision_12, (AR, 1) = AR_14, (*, 1) = *_17, (*, 3) = *_33, (NormalMeanPrecision, 1) = NormalMeanPrecision_8, (NormalMeanPrecision, 4) = NormalMeanPrecision_28, (AR, 3) = AR_30, (*, 2) = *_25, (NormalMeanPrecision, 5) = NormalMeanPrecision_36, (AR, 2) = AR_22}, {:γ = γ_1, :constvar_19 = constvar_19, :constvar_11 = constvar_11, :constvar_16 = constvar_16, :anonymous_var_graphppl = anonymous_var_graphppl_31, :constvar_24 = constvar_24, :constvar_35 = constvar_35, :s₀ = s₀_9, :constvar_32 = constvar_32, :constvar_6 = constvar_6, :constvar_27 = constvar_27, :constvar_7 = constvar_7, :constvar_2 = constvar_2, :θ = θ_5, :constvar_10 = constvar_10, :constvar_3 = constvar_3}, {:s = ResizableArray{GraphPPL.NodeLabel,1}(GraphPPL.NodeLabel[s_13, s_21, s_29]), :x = ResizableArray{GraphPPL.NodeLabel,1}(GraphPPL.NodeLabel[x_18, x_26, x_34])}, {}, {}, Base.RefValue{Any}(nothing)), and default weight 1.0
##
## fieldnames(GraphPPL.Context)
# lar_gppl_model_context = GraphPPL.getcontext(lar_gppl_model)
# lar_gppl_model_context.individual_variables
# lar_gppl_model_context.vector_variables
# lar_gppl_model_context.factor_nodes
print_meta_graph(lar_meta_graph)
============== NODES ==============
1: γ_1
2: constvar_2
3: constvar_3
4: GammaShapeRate_4
5: θ_5
6: constvar_6
7: constvar_7
8: NormalMeanPrecision_8
9: s₀_9
10: constvar_10
11: constvar_11
12: NormalMeanPrecision_12
13: s_13
14: AR_14
15: anonymous_var_graphppl_15
16: constvar_16
17: *_17
18: x_18
19: constvar_19
20: NormalMeanPrecision_20
21: s_21
22: AR_22
23: anonymous_var_graphppl_23
24: constvar_24
25: *_25
26: x_26
27: constvar_27
28: NormalMeanPrecision_28
29: s_29
30: AR_30
31: anonymous_var_graphppl_31
32: constvar_32
33: *_33
34: x_34
35: constvar_35
36: NormalMeanPrecision_36

============== LINKS ==============
γ_1 -> GammaShapeRate_4
γ_1 -> AR_14
γ_1 -> AR_22
γ_1 -> AR_30
constvar_2 -> GammaShapeRate_4
constvar_3 -> GammaShapeRate_4
θ_5 -> NormalMeanPrecision_8
θ_5 -> AR_14
θ_5 -> AR_22
θ_5 -> AR_30
constvar_6 -> NormalMeanPrecision_8
constvar_7 -> NormalMeanPrecision_8
s₀_9 -> NormalMeanPrecision_12
s₀_9 -> AR_14
constvar_10 -> NormalMeanPrecision_12
constvar_11 -> NormalMeanPrecision_12
s_13 -> AR_14
s_13 -> *_17
s_13 -> AR_22
anonymous_var_graphppl_15 -> *_17
anonymous_var_graphppl_15 -> NormalMeanPrecision_20
constvar_16 -> *_17
x_18 -> NormalMeanPrecision_20
constvar_19 -> NormalMeanPrecision_20
s_21 -> AR_22
s_21 -> *_25
s_21 -> AR_30
anonymous_var_graphppl_23 -> *_25
anonymous_var_graphppl_23 -> NormalMeanPrecision_28
constvar_24 -> *_25
x_26 -> NormalMeanPrecision_28
constvar_27 -> NormalMeanPrecision_28
s_29 -> AR_30
s_29 -> *_33
anonymous_var_graphppl_31 -> *_33
anonymous_var_graphppl_31 -> NormalMeanPrecision_36
constvar_32 -> *_33
x_34 -> NormalMeanPrecision_36
constvar_35 -> NormalMeanPrecision_36
meta_graph = lar_meta_graph

## Shorten some labels to make graph more readable
# str_labels = [string(lab) for lab in labels(meta_graph)]
# replacements = 
#     Pair("MvNormalMeanCovariance", "Nmc"), 
#     Pair("MvNormalMeanPrecision", "Nmp"), 
#     Pair("constvar", "cv"),
#     Pair("x", "XXXXX") ## make more obvious
# short_labels = [replace(s, replacements...) for s in str_labels]

GraphPlot.gplot( ## existing plotting functionality
    meta_graph,
    layout=spring_layout,
    nodelabel=collect(labels(meta_graph)),
    ## nodelabel=short_labels,
    nodelabelsize=0.1,
    NODESIZE=0.02, ## diameter of the nodes
    NODELABELSIZE=1.5,
    # nodelabelc="white",
    nodelabelc="green",
    nodelabeldist=0.0,
    nodefillc=nothing, ## "cyan"
    edgestrokec="red",
    ## ImageSize = (800, 800) ##- does not work
)
γ_1 constvar_2 constvar_3 GammaShapeRate_4 θ_5 constvar_6 constvar_7 NormalMeanPrecision_8 s₀_9 constvar_10 constvar_11 NormalMeanPrecision_12 s_13 AR_14 anonymous_var_graphppl_15 constvar_16 *_17 x_18 constvar_19 NormalMeanPrecision_20 s_21 AR_22 anonymous_var_graphppl_23 constvar_24 *_25 x_26 constvar_27 NormalMeanPrecision_28 s_29 AR_30 anonymous_var_graphppl_31 constvar_32 *_33 x_34 constvar_35 NormalMeanPrecision_36 
_san_node_ids = Dict("dummy" => "dummy")
_lup_enabled = true
_lup_agc = false ## append global_counter values
_raw_links_enabled = false
_CFFG = true
_deltas_enabled = true
_tsep = 2*_xsep
_lup_dict = Dict{String, String}(
    "θ_7" => _lup_agc ? "\\theta\\_7" : "\\theta",
    "constvar_8_9" => _lup_agc ? "\\mu\\_8\\_9" : "\\mu",
    "constvar_10_11" => _lup_agc ? "\\gamma\\_10\\_11" : "\\gamma",
    "NormalMeanPrecision_12" => _lup_agc ? "\\mathcal{N}\\_12" : "\\mathcal{N}",

    "constvar_2_3" => _lup_agc ? "\\alpha\\_2\\_3" : "\\alpha",
    "constvar_4_5" => _lup_agc ? "\\beta\\_4\\_5" : "\\beta",
    "GammaShapeRate_6" => _lup_agc ? "\\mathcal{G}am\\_6" : "\\mathcal{G}am",
    "γ_1" => _lup_agc ? "\\gamma\\_1" : "\\gamma",

    "constvar_14_15" => _lup_agc ? "\\mu\\_14\\_15" : "\\mu",
    "constvar_16_17" => _lup_agc ? "\\gamma\\_16\\_17" : "\\gamma",
    "NormalMeanPrecision_18" => _lup_agc ? "\\mathcal{N}\\_18" : "\\mathcal{N}",
    ## "s₀_13" => "s_0", ## s₀ must always be in _lup_dict
    "s₀_13" => _lup_agc ? "\\mathbf{s}_0\\_13" : "\\mathbf{s}_0",
    "s_19" => _lup_agc ? "\\mathbf{s}_1\\_19" : "\\mathbf{s}_1",
    "s_29" => _lup_agc ? "\\mathbf{s}_2\\_29" : "\\mathbf{s}_2",
    "s_39" => _lup_agc ? "\\mathbf{s}_3\\_39" : "\\mathbf{s}_3",

    "AR_20" => _lup_agc ? "AR\\_20" : "AR",
    "AR_30" => _lup_agc ? "AR\\_30" : "AR",
    "AR_40" => _lup_agc ? "AR\\_40" : "AR",

    "constvar_22_23" => _lup_agc ? "\\mathbf{B}\\_22\\_23" : "\\mathbf{B}",
    "constvar_32_33" => _lup_agc ? "\\mathbf{B}\\_32\\_33" : "\\mathbf{B}",
    "constvar_42_43" => _lup_agc ? "\\mathbf{B}\\_42\\_43" : "\\mathbf{B}",

    "*_24" => _lup_agc ? "*\\_24" : "*",
    "*_34" => _lup_agc ? "*\\_34" : "*",
    "*_44" => _lup_agc ? "*\\_44" : "*",
    
    "anonymous_var_graphppl_21" => _lup_agc ? "\\mathbf{v}_u\\_21" : "\\mathbf{v}_u",
    "anonymous_var_graphppl_31" => _lup_agc ? "\\mathbf{v}_u\\_31" : "\\mathbf{v}_u",
    "anonymous_var_graphppl_41" => _lup_agc ? "\\mathbf{v}_u\\_41" : "\\mathbf{v}_u",

    "constvar_26_27" => _lup_agc ? "\\tau\\_26\\_27" : "\\tau",
    "constvar_36_37" => _lup_agc ? "\\tau\\_36\\_37" : "\\tau",
    "constvar_46_47" => _lup_agc ? "\\tau\\_46\\_47" : "\\tau",

    "NormalMeanPrecision_28" => _lup_agc ? "\\mathcal{N}\\_28" : "\\mathcal{N}",
    "NormalMeanPrecision_38" => _lup_agc ? "\\mathcal{N}\\_38" : "\\mathcal{N}",
    "NormalMeanPrecision_48" => _lup_agc ? "\\mathcal{N}\\_48" : "\\mathcal{N}",

    "x_25" => _lup_agc ? "x_1\\_25" : "x_1",
    "x_35" => _lup_agc ? "x_2\\_35" : "x_2",
    "x_45" => _lup_agc ? "x_3\\_45" : "x_3",
);
lar_kobus_graph = create_kobus_graph(lar_gppl_model)
1: γ_1
2: constvar_2
3: constvar_3
4: GammaShapeRate_4
5: θ_5
6: constvar_6
7: constvar_7
8: NormalMeanPrecision_8
9: s₀_9
10: constvar_10
11: constvar_11
12: NormalMeanPrecision_12
13: s_13
14: AR_14
15: anonymous_var_graphppl_15
16: constvar_16
17: *_17
18: x_18
19: constvar_19
20: NormalMeanPrecision_20
21: s_21
22: AR_22
23: anonymous_var_graphppl_23
24: constvar_24
25: *_25
26: x_26
27: constvar_27
28: NormalMeanPrecision_28
29: s_29
30: AR_30
31: anonymous_var_graphppl_31
32: constvar_32
33: *_33
34: x_34
35: constvar_35
36: NormalMeanPrecision_36
Dict{String, Dict{Symbol, Any}} with 36 entries:
  "GammaShapeRate_4"          => Dict(:node_lup=>"GammaShapeRate\\_4", :node_co…
  "s_13"                      => Dict(:node_lup=>"s\\_13", :node_code=>13, :nod…
  "anonymous_var_graphppl_15" => Dict(:node_lup=>"anonymous\\_var\\_graphppl\\_…
  "constvar_19"               => Dict(:node_lup=>"constvar\\_19", :node_code=>1…
  "constvar_35"               => Dict(:node_lup=>"constvar\\_35", :node_code=>3…
  "s₀_9"                      => Dict(:node_lup=>"s₀\\_9", :node_code=>9, :node…
  "AR_14"                     => Dict(:node_lup=>"AR\\_14", :node_code=>14, :no…
  "NormalMeanPrecision_28"    => Dict(:node_lup=>"\\mathcal{N}", :node_code=>28…
  "constvar_24"               => Dict(:node_lup=>"constvar\\_24", :node_code=>2…
  "constvar_2"                => Dict(:node_lup=>"constvar\\_2", :node_code=>2,…
  "NormalMeanPrecision_8"     => Dict(:node_lup=>"NormalMeanPrecision\\_8", :no…
  "*_17"                      => Dict(:node_lup=>"*\\_17", :node_code=>17, :nod…
  "AR_30"                     => Dict(:node_lup=>"AR", :node_code=>30, :node_da…
  "*_33"                      => Dict(:node_lup=>"*\\_33", :node_code=>33, :nod…
  "anonymous_var_graphppl_23" => Dict(:node_lup=>"anonymous\\_var\\_graphppl\\_…
  "x_18"                      => Dict(:node_lup=>"x\\_18", :node_code=>18, :nod…
  "x_34"                      => Dict(:node_lup=>"x\\_34", :node_code=>34, :nod…
  "θ_5"                       => Dict(:node_lup=>"θ\\_5", :node_code=>5, :node_…
  "constvar_3"                => Dict(:node_lup=>"constvar\\_3", :node_code=>3,…
  ⋮                           => ⋮
layers = Dict()
kobus_graph = lar_kobus_graph
gppl_model = lar_gppl_model

paramB1_inifac = find_inifac_node(kobus_graph, "θ", gppl_model)
paramB1_out_var, paramB1_iniparvars = find_iniparvar_nodes(kobus_graph, paramB1_inifac, gppl_model)
layers["paramB1"] = Layer(
    name="paramB1",
    anchor=0*_xsep,
    ## iniparvars=["constvar_8_9", "constvar_10_11"],
    iniparvars=paramB1_iniparvars,
    ## inifac="NormalMeanPrecision_12",
    inifac=paramB1_inifac,
    ## parvar="θ_7",
    parvar=paramB1_out_var,
    sysvars = [],
    facs = [],
    parvars = [],
    link_fac_side="north", link_fac_offset=_divn[3][1])
## layers["paramB1"].links = layers["state"].facs

paramB2_inifac = find_inifac_node(kobus_graph, "γ", gppl_model)
paramB2_out_var, paramB2_iniparvars = find_iniparvar_nodes(kobus_graph, paramB2_inifac, gppl_model)
layers["paramB2"] = Layer(
    name="paramB2",
    anchor=0*_xsep,
    ## iniparvars = ["constvar_2_3", "constvar_4_5"],
    iniparvars = paramB2_iniparvars,
    ## inifac = "GammaShapeRate_6",
    inifac = paramB2_inifac,
    ## parvar = "γ_1",
    parvar = paramB2_out_var,
    sysvars = [],
    facs = [],
    parvars = [],
    link_fac_side="north", link_fac_offset=_divn[3][3])
## layers["paramB2"].links = layers["state"].facs

state_inifac = find_inifac_node(kobus_graph, "s₀", gppl_model)
_, state_vars, state_facs = find_state_nodes(kobus_graph, state_inifac, gppl_model)
state_out_var, state_iniparvars = find_iniparvar_nodes(kobus_graph, state_inifac, gppl_model)
layers["state"] = Layer(
    name="state",
    anchor=3*_xsep,
    ## iniparvars = ["constvar_14_15", "constvar_16_17"],
    iniparvars = state_iniparvars,
    ## inifac = "NormalMeanPrecision_18",
    inifac = state_inifac,
    parvar="",
    ## sysvars = ["s₀_13", "s_19", "s_29", "s_39"],
    sysvars = state_vars, ## contains only sysvars
    ## facs = [["AR_20", "AR_30", "AR_40"]],
    facs = [state_facs],
    parvars = [])
## layers["state"].links = layers["stateD1"].facs

layers["stateD1"] = Layer(
    name="stateD1",
    anchor=4*_xsep,
    iniparvars = [],
    inifac = "",
    parvar="",
    sysvars = ["anonymous_var_graphppl_21", "anonymous_var_graphppl_31", "anonymous_var_graphppl_41"],
    facs = [["*_24", "*_34", "*_44"]],
    parvars = [
        [["constvar_22_23"], ["constvar_32_33"], ["constvar_42_43"]] ],
    link_fac_side="north", link_fac_offset=_divn[1][1]
)
## layers["stateD1"].links = layers["obser"].facs

obser_sysvars, obser_facs, obser_parvars = find_obser_nodes(kobus_graph, "x", gppl_model)
layers["obser"] = Layer(
    name="obser",
    anchor=4*_xsep,
    iniparvars = [],
    inifac = "",
    parvar="",
    ## MANUAL:  ## state vars not directly connected to obser facs
    sysvars = ["x_25", "x_35", "x_45"],
    ## MANUAL:  ## state vars not directly connected to obser facs
    ## facs = [obser_facs], 
    facs = [["NormalMeanPrecision_28", "NormalMeanPrecision_38", "NormalMeanPrecision_48"]],
    ## MANUAL: sysvars = obser_sysvars, ## state vars not directly connected to obser facs
    ## parvars = [obser_parvars],
    parvars = [
        [["constvar_26_27"], ["constvar_36_37"], ["constvar_46_47"]] ])

layers["paramB1"].links = layers["state"].facs
layers["paramB2"].links = layers["state"].facs
layers["state"].links = layers["stateD1"].facs
layers["stateD1"].links = layers["obser"].facs
    var_neighbours[1]: NormalMeanPrecision_8
        in_neighbor_label: constvar_6
        out_neighbor_label: θ_5
inifac_out_neighbor_str: θ_5
inifac_str: NormalMeanPrecision_8
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(θ_5, out, NodeData in context  with properties name = θ, index = nothing), (constvar_6, μ, NodeData in context  with properties name = constvar, index = nothing), (constvar_7, τ, NodeData in context  with properties name = constvar, index = nothing)]
        out_var: θ_5
    var_neighbours[1]: GammaShapeRate_4
        in_neighbor_label: constvar_2
        out_neighbor_label: γ_1
inifac_out_neighbor_str: γ_1
inifac_str: GammaShapeRate_4
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(γ_1, out, NodeData in context  with properties name = γ, index = nothing), (constvar_2, α, NodeData in context  with properties name = constvar, index = nothing), (constvar_3, β, NodeData in context  with properties name = constvar, index = nothing)]
        out_var: γ_1
    var_neighbours[1]: NormalMeanPrecision_12
        in_neighbor_label: constvar_10
        out_neighbor_label: s₀_9
inifac_out_neighbor_str: s₀_9
inifac_str: NormalMeanPrecision_12
var: (s₀_9, out, NodeData in context  with properties name = s₀, index = nothing)
=== NEXT var node is: s₀_9
var_label: s₀_9
var_index: 0
var_neighbours: Any["NormalMeanPrecision_12", "AR_14"]
    var_neighbours[1]: NormalMeanPrecision_12
        in_neighbor_label: constvar_10
        out_neighbor_label: s₀_9
    var_neighbours[2]: AR_14
        in_neighbor_label: s₀_9
        out_neighbor_label: s_13
        === NEXT fac node is: AR_14
var: (s_13, y, NodeData in context  with properties name = s, index = 1)
=== NEXT var node is: s_13
var_label: s_13
var_index: 1
var_neighbours: Any["AR_14", "*_17", "AR_22"]
    var_neighbours[1]: AR_14
        in_neighbor_label: s₀_9
        out_neighbor_label: s_13
    var_neighbours[2]: *_17
        in_neighbor_label: constvar_16
        out_neighbor_label: anonymous_var_graphppl_15
    var_neighbours[3]: AR_22
        in_neighbor_label: s_13
        out_neighbor_label: s_21
        === NEXT fac node is: AR_22
var: (s_21, y, NodeData in context  with properties name = s, index = 2)
=== NEXT var node is: s_21
var_label: s_21
var_index: 2
var_neighbours: Any["AR_22", "*_25", "AR_30"]
    var_neighbours[1]: AR_22
        in_neighbor_label: s_13
        out_neighbor_label: s_21
    var_neighbours[2]: *_25
        in_neighbor_label: constvar_24
        out_neighbor_label: anonymous_var_graphppl_23
    var_neighbours[3]: AR_30
        in_neighbor_label: s_21
        out_neighbor_label: s_29
        === NEXT fac node is: AR_30
var: (s_29, y, NodeData in context  with properties name = s, index = 3)
=== NEXT var node is: s_29
var_label: s_29
var_index: 3
var_neighbours: Any["AR_30", "*_33"]
    var_neighbours[1]: AR_30
        in_neighbor_label: s_21
        out_neighbor_label: s_29
    var_neighbours[2]: *_33
        in_neighbor_label: constvar_32
        out_neighbor_label: anonymous_var_graphppl_31
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(s₀_9, out, NodeData in context  with properties name = s₀, index = nothing), (constvar_10, μ, NodeData in context  with properties name = constvar, index = nothing), (constvar_11, τ, NodeData in context  with properties name = constvar, index = nothing)]
        out_var: s₀_9
obser_var_label_strs: ["x_18", "x_26", "x_34"]
obser_var_label_strs[i]: x_18
    var_neighbours[j]: NormalMeanPrecision_20
        in_neighbor_label: anonymous_var_graphppl_15
        out_neighbor_label: x_18
obser_var_label_strs[i]: x_26
    var_neighbours[j]: NormalMeanPrecision_28
        in_neighbor_label: anonymous_var_graphppl_23
        out_neighbor_label: x_26
obser_var_label_strs[i]: x_34
    var_neighbours[j]: NormalMeanPrecision_36
        in_neighbor_label: anonymous_var_graphppl_31
        out_neighbor_label: x_34
1-element Vector{Vector{String}}:
 ["NormalMeanPrecision_28", "NormalMeanPrecision_38", "NormalMeanPrecision_48"]
lar_tikz = generate_tikz(
    heading = "Time-Varying Autoregressive Model (Univariate)",
    Model = lar_gppl_model,
    layers = layers) ## fac layers, i.e. not var (vars are handled with assoced layer)
## print(lar_tikz)
show_tikz(lar_tikz); ## write to file rather than show in notebook
\documentclass{standalone}
\usepackage{tikz}
\usetikzlibrary{calc}
\usepackage{amsmath}
\usepackage{bm}
\newcommand\ns{10mm} %% define node (minimum) size
\newcommand\nr{.5*\ns} %% define node radius
\newcommand\ar{.4*\ns} %% define arc radius
\newcommand\br{.19*\ns} %% define bead radius
\DeclareUnicodeCharacter{03B8}{\theta}
\DeclareUnicodeCharacter{03B3}{\gamma}
\DeclareUnicodeCharacter{209C}{\0}
\DeclareUnicodeCharacter{208B}{\0}
\DeclareUnicodeCharacter{2081}{\0}
\begin{document}
\begin{tikzpicture}[
    scale=1, 
    draw=black, 
    fac/.style={rectangle, draw=black!100, fill=black!0, minimum size=\ns, inner sep=0pt},
    dlt/.style={rectangle, draw=black!100, fill=black!100, minimum size=.2*\ns},
    var/.style={circle, draw=gray!100, fill=gray!100, minimum size=1mm, inner sep=0pt},
    bus/.style={rectangle, draw=black!100, fill=black!100, minimum width=30mm, minimum height=.5mm, inner sep=0pt},
    eql/.style={rectangle, draw=black!100, fill=black!0, minimum size=.5*\ns},
    text=black,
    ]
\node[font=\Large] at (4,29) {\textbf{Time-Varying Autoregressive Model (Univariate)}};
⋮
⋮
⋮
\path[use as bounding box] ([xshift=-1cm, yshift=-1cm]current bounding box.south west) rectangle ([xshift=1cm, yshift=2cm]current bounding box.north east);
\end{tikzpicture}
\end{document}
seq: ["AR_14"]
... doing the else ...
seq: ["AR_22"]
... doing the else ...
seq: ["AR_30"]
... doing the else ...

Time-Varying Autoregressive Model (Univariate)

Fly a Drone to a Target using Active Inference

We start by specifying a probabilistic model for the agent that describes the agent’s internal beliefs over the external dynamics of the environment.

The generative model is defined as follows:

\[\begin{aligned} p'(\mathbf{x}_:, \mathbf{s}_:, \mathbf{u}_:) = \underbrace{p(\mathbf{s}_{t-1})}_{\substack{ \text{Initial} \\ \text{state} \\ \text{prior}}} \prod_{k=t}^{t+T} \underbrace{p(\mathbf{x}_k \mid \mathbf{s}_k)}_{\substack{ \text{Observation} \\ \text{model}}} \, \underbrace{p(\mathbf{s}_k \mid \mathbf{s}_{k-1}, \mathbf{u}_k}_{\substack{\text{Transition} \\ \text{model}}}) \, \underbrace{p(\mathbf{u}_k)}_{\substack { \text{Control} \\ \text{prior}}} \ \nonumber \qquad \mathrm{(2019\_VanDeLaar,\_SAIPbMP\ eq9)} \end{aligned}\]

The generative model includes future time steps.

To infer goal-driven (i.e. purposeful) behavior, we add prior beliefs \(p^+(\mathbf{x})\) about desired future observations. This leads to an extended agent model:

\[\begin{aligned} p(\mathbf{x}_:, \mathbf{s}_:, \mathbf{u}_:) &= \frac{p'(\mathbf{x}, \mathbf{s}, \mathbf{u}) p^+(\mathbf{x})}{\int_x p'(\mathbf{x}, \mathbf{s}, \mathbf{u})p^+(\mathbf{x}) \mathrm{d}\mathbf{x}} \\ &\propto \underbrace{p(\mathbf{s}_{t-1})\prod_{k=t}^{t+T} p(x_k \mid \mathbf{s}_k)p(\mathbf{s}_k \mid \mathbf{s}_{k-1}, \mathbf{u}_k)p(\mathbf{u}_k)}_{\text{original generative model}} \, \underbrace{p^+(\mathbf{x}_k)}_{\substack { \text{extension} \\ \text{Goal} \\ \text{prior}}} \\ &\propto \underbrace{p(\mathbf{s}_{t-1})}_{\substack{ \text{Initial} \\ \text{state} \\ \text{prior}}} \prod_{k=t}^{t+T} \underbrace{p(\mathbf{x}_k \mid \mathbf{s}_k)}_{\substack{ \text{Observation} \\ \text{model}}} \, \underbrace{p(\mathbf{s}_k \mid \mathbf{s}_{k-1}, \mathbf{u}_k}_{\substack{\text{Transition} \\ \text{model}}}) \, \underbrace{p(\mathbf{u}_k)}_{\substack { \text{Control} \\ \text{prior}}} \, \underbrace{p^+(\mathbf{x}_k)}_{\substack { \text{Goal} \\ \text{prior}}} \nonumber \end{aligned}\]

The factors are defined as:

observation:

\[p(\mathbf{x}_k \mid \mathbf{s}_k) = \mathcal{N}(\mathbf{x}_k \mid \mathbf{s}_k,\,\mathbf\Theta)\] where \(\mathbf{x}_k = (\chi_{1k}, \chi_{2k}, ...)\) denotes observations of the agent after interacting with the environment.

state transition:

\[ \begin{aligned} p(\mathbf{s}_k \mid \mathbf{s}_{k-1}, \mathbf{u}_k) &= \mathcal{N}(\mathbf{s}_k \mid \mathbf{s}_{k-1} + \mathbf{u}_k, ϑ) \\ p(\mathbf{s}_{t-1}) &= \mathcal{N}(\mathbf{s}_{t-1} \mid \mathbf{m}_{t-1},\, \mathbf{V}_{t-1}) \end{aligned} \]

The current state is a linear function of the previous state and action. We have endowed the agent with an accurate model of the system dynamics.

control:

\[ \begin{aligned} p(\mathbf{u}_t) &= \prod_{k=t}^{t+T} \mathcal{N}(\mathbf{u}_k \mid \mathbf{0},\,\mathbf\Xi) \\ &= \prod_{k=t}^{t+T} \mathcal{N}(\mathbf{u}_k \mid \mathbf{0},\,ξ ⋅ \mathbf{I}) \\ &= \prod_{k=t}^{t+T} \mathcal{N}(\mathbf{u}_k \mid \mathbf{0},\,1.0) \\ &= \mathcal{N}(\mathbf{u}_t \mid \mathbf{m}_u,\, \mathbf{V}_u)\\ \end{aligned} \]

This represents the control priors.

goal/target/preference:

\[ \begin{aligned} p^+(\mathbf{x}_t) &= \prod_{k=t}^{t+(T-1)} \mathcal{N}(\mathbf{x}_k \mid \mathbf{0},\,\sigma^{huge}) \cdot \mathcal{N}(\mathbf{x}_T \mid \mathbf{x}_+, \, \mathbf\Sigma)\\ &= \prod_{k=t}^{t+(T-1)} \mathcal{N}(\mathbf{x}_k \mid \mathbf{0},\,\sigma^{huge}) \cdot \mathcal{N}(\mathbf{x}_T \mid \mathbf{x}_+, \, \sigma ⋅ \mathbf I)\\ &= \prod_{k=t}^{t+(T-1)} \mathcal{N}(\mathbf{x}_k \mid \mathbf{0},\,10^{12}) \cdot \mathcal{N}(\mathbf{x}_T \mid [0.0, 0.0, 0.0*π, 0.1], \, 10^{-4} ⋅ \mathbf I)\\ &= \mathcal{N}(\mathbf{x}_t \mid \mathbf{m}_x,\, \mathbf{V}_x) \\ \end{aligned} \]

This represents the goal/target/preference priors and encodes a belief about a preferred target/goal \(\mathbf{x}₊\).

initial state:

Setting \(t=1\), \[ p(\mathbf{s}_0) = \mathcal{N}(\mathbf{s}_0 \mid 0, 10^{12}) \]

This means we set a vague prior for the initial state.

4.5.3.1 Generative Model for the Drone

The code in the next block defines the agent’s internal beliefs over the external dynamics and its probabilistic model of the environment, which correspond accurately by directly using the functions defined above. We use the @model macro from RxInfer to define the probabilistic model and the meta block to define approximation methods for the nonlinear state-transition functions.

In the model specification we in addition to the current state of the agent we include the beliefs over its future states (up to T steps ahead):

@model function dronenav_model(x, mᵤ, Vᵤ, mₓ, Vₓ, mₛ₍ₜ₋₁₎, Vₛ₍ₜ₋₁₎, T, Rᵃ)
    ## Transition function
    g = (sₜ₋₁::AbstractVector) -> begin
        sₜ = similar(sₜ₋₁) ## Next state
        sₜ = Aᵃ(sₜ₋₁, 1.0) + sₜ₋₁
        return sₜ
    end
    
    ## Function for modeling turn/yaw control
    h = (u::AbstractVector) -> Rᵃ(u[1])
    
    Γ = _γ*diageye(4) ## Transition precision
    𝚯 = _ϑ*diageye(4) ## Observation variance
    
    ## sₜ₋₁ ~ MvNormal(mean=mₛ₍ₜ₋₁₎, cov=Vₛ₍ₜ₋₁₎)
    s₀ ~ MvNormal(mean=mₛ₍ₜ₋₁₎, cov=Vₛ₍ₜ₋₁₎)
    ## sₖ₋₁ = sₜ₋₁
    sₖ₋₁ = s₀
    
    local s

    for k in 1:T
        ## Control
        u[k] ~ MvNormal(mean=mᵤ[k], cov=Vᵤ[k])
        hIuI[k] ~ h(u[k]) where { meta=DeltaMeta(method=Unscented()) }

        ## State transition
        gIsI[k] ~ g(sₖ₋₁) where { meta=DeltaMeta(method=Unscented()) }
        ghSum[k] ~ gIsI[k] + hIuI[k]#.
        s[k] ~ MvNormal(mean=ghSum[k], precision=Γ)

        ## Likelihood of future observations
        x[k] ~ MvNormal(mean=s[k], cov=𝚯)

        ## Target/Goal prior
        x[k] ~ MvNormal(mean=mₓ[k], cov=Vₓ[k])

        sₖ₋₁ = s[k]
    end
    return (s, )
end
_Fᴱⁿᵍᴸⁱᵐⁱᵗ = 0.1
function Rᵃ(a::Real) ## turn/yaw rate
    b = [ 0.0, 0.0, 1.0, 0.0 ]
    return b*_Fᴱⁿᵍᴸⁱᵐⁱᵗ*tanh(a)
end
## Rᵃ(0.25)

_γ = 1e4 ## transition precision (system noise)
_ϑ = 1e-4 ## observation variance (observation noise)

## T =_Tᵃⁱ,
## T =100
T =3
Rᵃ=Rᵃ

mᵤ = Vector{Float64}[ [0.0] for k=1:T ] ##Set control priors
ξ = 0.1
Ξ  = fill(ξ, 1, 1) ##Control prior variance
Vᵤ = Matrix{Float64}[ Ξ for k=1:T ]
mₓ      = [zeros(4) for k=1:T]
x₊ = [0.0, 0.0, 0.0*π, 0.1] ## Target/Goal state
mₓ[end] = x₊ ##Set prior mean to reach target/goal at t=T
Vₓ      = [huge*diageye(4) for k=1:T]
σ = 1e-4
Σ       = σ*diageye(4) ##Target/Goal prior variance
Vₓ[end] = Σ ##Set prior variance to reach target/goal at t=T
s₀ = [8.0, 8.0, -0.1, 0.1] ## initial state
mₛ₍ₜ₋₁₎ = s₀
Vₛ₍ₜ₋₁₎ = tiny*diageye(4)

drone_conditioned = dronenav_model(
    mᵤ= mᵤ, 
    Vᵤ= Vᵤ, 
    mₓ= mₓ, 
    Vₓ= Vₓ,
    mₛ₍ₜ₋₁₎= mₛ₍ₜ₋₁₎,
    Vₛ₍ₜ₋₁₎= Vₛ₍ₜ₋₁₎,
    T=T, 
    Rᵃ=Rᵃ
) | (x = [ [8.099, 7.990, -0.109, 0.1],  [8.198, 7.979, -0.119, 0.1],  [8.298, 7.967, -0.129, 0.1]],)
## ) | (x = [ [8.099, 7.990, -0.109],  [8.198, 7.979, -0.119],  [8.298, 7.967, -0.129]],)
dronenav_model(mᵤ = [[0.0], [0.0], [0.0]], Vᵤ = [[0.1;;], [0.1;;], [0.1;;]], mₓ = [[0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.1]], Vₓ = [[1.0e12 0.0 0.0 0.0; 0.0 1.0e12 0.0 0.0; 0.0 0.0 1.0e12 0.0; 0.0 0.0 0.0 1.0e12], [1.0e12 0.0 0.0 0.0; 0.0 1.0e12 0.0 0.0; 0.0 0.0 1.0e12 0.0; 0.0 0.0 0.0 1.0e12], [0.0001 0.0 0.0 0.0; 0.0 0.0001 0.0 0.0; 0.0 0.0 0.0001 0.0; 0.0 0.0 0.0 0.0001]], mₛ₍ₜ₋₁₎ = [8.0, 8.0, -0.1, 0.1], Vₛ₍ₜ₋₁₎ = [1.0e-12 0.0 0.0 0.0; 0.0 1.0e-12 0.0 0.0; 0.0 0.0 1.0e-12 0.0; 0.0 0.0 0.0 1.0e-12], T = 3, Rᵃ = Rᵃ) conditioned on: 
  x = [[8.099, 7.99, -0.109, 0.1], [8.198, 7.979, -0.119, 0.1], [8.298, 7.967, -0.129, 0.1]]
drone_rxi_model = RxInfer.create_model(drone_conditioned)
drone_gppl_model = RxInfer.getmodel(drone_rxi_model)
drone_meta_graph = drone_gppl_model.graph
Meta graph based on a Graphs.SimpleGraphs.SimpleGraph{Int64} with vertex labels of type GraphPPL.NodeLabel, vertex metadata of type GraphPPL.NodeData, edge metadata of type GraphPPL.EdgeLabel, graph metadata given by GraphPPL.Context(0, dronenav_model, "", nothing, {var"#69#81"() = 3, var"#70#82"{typeof(Rᵃ)}(Rᵃ) = 3, MvNormalMeanPrecision = 3, MvNormalMeanCovariance = 10, + = 3}, {}, {(#69, 1) = #69_12, (MvNormalMeanCovariance, 2) = MvNormalMeanCovariance_8, (MvNormalMeanPrecision, 1) = MvNormalMeanPrecision_17, (MvNormalMeanCovariance, 9) = MvNormalMeanCovariance_58, (MvNormalMeanCovariance, 1) = MvNormalMeanCovariance_4, (+, 2) = +_33, (#70, 3) = #70_48, (+, 3) = +_52, (MvNormalMeanCovariance, 4) = MvNormalMeanCovariance_23, (MvNormalMeanCovariance, 8) = MvNormalMeanCovariance_46, (#69, 2) = #69_31, (MvNormalMeanCovariance, 6) = MvNormalMeanCovariance_39, (#70, 2) = #70_29, (MvNormalMeanCovariance, 5) = MvNormalMeanCovariance_27, (MvNormalMeanCovariance, 10) = MvNormalMeanCovariance_61, (+, 1) = +_14, (#70, 1) = #70_10, (MvNormalMeanCovariance, 7) = MvNormalMeanCovariance_42, (#69, 3) = #69_50, (MvNormalMeanPrecision, 2) = MvNormalMeanPrecision_36, (MvNormalMeanPrecision, 3) = MvNormalMeanPrecision_55, (MvNormalMeanCovariance, 3) = MvNormalMeanCovariance_20}, {:constvar_16 = constvar_16, :constvar_22 = constvar_22, :constvar_54 = constvar_54, :constvar_35 = constvar_35, :constvar_21 = constvar_21, :constvar_26 = constvar_26, :constvar_59 = constvar_59, :constvar_44 = constvar_44, :s₀ = s₀_1, :constvar_60 = constvar_60, :constvar_19 = constvar_19, :constvar_45 = constvar_45, :constvar_6 = constvar_6, :constvar_2 = constvar_2, :constvar_57 = constvar_57, :constvar_3 = constvar_3, :constvar_41 = constvar_41, :constvar_7 = constvar_7, :constvar_40 = constvar_40, :constvar_25 = constvar_25, :constvar_38 = constvar_38}, {:x = ResizableArray{GraphPPL.NodeLabel,1}(GraphPPL.NodeLabel[x_18, x_37, x_56]), :s = ResizableArray{GraphPPL.NodeLabel,1}(GraphPPL.NodeLabel[s_15, s_34, s_53]), :ghSum = ResizableArray{GraphPPL.NodeLabel,1}(GraphPPL.NodeLabel[ghSum_13, ghSum_32, ghSum_51]), :u = ResizableArray{GraphPPL.NodeLabel,1}(GraphPPL.NodeLabel[u_5, u_24, u_43]), :gIsI = ResizableArray{GraphPPL.NodeLabel,1}(GraphPPL.NodeLabel[gIsI_11, gIsI_30, gIsI_49]), :hIuI = ResizableArray{GraphPPL.NodeLabel,1}(GraphPPL.NodeLabel[hIuI_9, hIuI_28, hIuI_47])}, {}, {}, Base.RefValue{Any}((ResizableArray{GraphPPL.NodeLabel,1}(GraphPPL.NodeLabel[s_15, s_34, s_53]),))), and default weight 1.0
##
## fieldnames(GraphPPL.Context)
# drone_gppl_model_context = GraphPPL.getcontext(drone_gppl_model)
# drone_gppl_model_context.individual_variables
# drone_gppl_model_context.vector_variables
# drone_gppl_model_context.factor_nodes
print_meta_graph(drone_meta_graph)
============== NODES ==============
1: s₀_1
2: constvar_2
3: constvar_3
4: MvNormalMeanCovariance_4
5: u_5
6: constvar_6
7: constvar_7
8: MvNormalMeanCovariance_8
9: hIuI_9
10: #70_10
11: gIsI_11
12: #69_12
13: ghSum_13
14: +_14
15: s_15
16: constvar_16
17: MvNormalMeanPrecision_17
18: x_18
19: constvar_19
20: MvNormalMeanCovariance_20
21: constvar_21
22: constvar_22
23: MvNormalMeanCovariance_23
24: u_24
25: constvar_25
26: constvar_26
27: MvNormalMeanCovariance_27
28: hIuI_28
29: #70_29
30: gIsI_30
31: #69_31
32: ghSum_32
33: +_33
34: s_34
35: constvar_35
36: MvNormalMeanPrecision_36
37: x_37
38: constvar_38
39: MvNormalMeanCovariance_39
40: constvar_40
41: constvar_41
42: MvNormalMeanCovariance_42
43: u_43
44: constvar_44
45: constvar_45
46: MvNormalMeanCovariance_46
47: hIuI_47
48: #70_48
49: gIsI_49
50: #69_50
51: ghSum_51
52: +_52
53: s_53
54: constvar_54
55: MvNormalMeanPrecision_55
56: x_56
57: constvar_57
58: MvNormalMeanCovariance_58
59: constvar_59
60: constvar_60
61: MvNormalMeanCovariance_61

============== LINKS ==============
s₀_1 -> MvNormalMeanCovariance_4
s₀_1 -> #69_12
constvar_2 -> MvNormalMeanCovariance_4
constvar_3 -> MvNormalMeanCovariance_4
u_5 -> MvNormalMeanCovariance_8
u_5 -> #70_10
constvar_6 -> MvNormalMeanCovariance_8
constvar_7 -> MvNormalMeanCovariance_8
hIuI_9 -> #70_10
hIuI_9 -> +_14
gIsI_11 -> #69_12
gIsI_11 -> +_14
ghSum_13 -> +_14
ghSum_13 -> MvNormalMeanPrecision_17
s_15 -> MvNormalMeanPrecision_17
s_15 -> MvNormalMeanCovariance_20
s_15 -> #69_31
constvar_16 -> MvNormalMeanPrecision_17
x_18 -> MvNormalMeanCovariance_20
x_18 -> MvNormalMeanCovariance_23
constvar_19 -> MvNormalMeanCovariance_20
constvar_21 -> MvNormalMeanCovariance_23
constvar_22 -> MvNormalMeanCovariance_23
u_24 -> MvNormalMeanCovariance_27
u_24 -> #70_29
constvar_25 -> MvNormalMeanCovariance_27
constvar_26 -> MvNormalMeanCovariance_27
hIuI_28 -> #70_29
hIuI_28 -> +_33
gIsI_30 -> #69_31
gIsI_30 -> +_33
ghSum_32 -> +_33
ghSum_32 -> MvNormalMeanPrecision_36
s_34 -> MvNormalMeanPrecision_36
s_34 -> MvNormalMeanCovariance_39
s_34 -> #69_50
constvar_35 -> MvNormalMeanPrecision_36
x_37 -> MvNormalMeanCovariance_39
x_37 -> MvNormalMeanCovariance_42
constvar_38 -> MvNormalMeanCovariance_39
constvar_40 -> MvNormalMeanCovariance_42
constvar_41 -> MvNormalMeanCovariance_42
u_43 -> MvNormalMeanCovariance_46
u_43 -> #70_48
constvar_44 -> MvNormalMeanCovariance_46
constvar_45 -> MvNormalMeanCovariance_46
hIuI_47 -> #70_48
hIuI_47 -> +_52
gIsI_49 -> #69_50
gIsI_49 -> +_52
ghSum_51 -> +_52
ghSum_51 -> MvNormalMeanPrecision_55
s_53 -> MvNormalMeanPrecision_55
s_53 -> MvNormalMeanCovariance_58
constvar_54 -> MvNormalMeanPrecision_55
x_56 -> MvNormalMeanCovariance_58
x_56 -> MvNormalMeanCovariance_61
constvar_57 -> MvNormalMeanCovariance_58
constvar_59 -> MvNormalMeanCovariance_61
constvar_60 -> MvNormalMeanCovariance_61
meta_graph = drone_meta_graph

## Shorten some labels to make graph more readable
str_labels = [string(lab) for lab in labels(meta_graph)]
replacements = 
    Pair("MvNormalMeanCovariance", "Nmc"), 
    Pair("MvNormalMeanPrecision", "Nmp"), 
    Pair("constvar", "cv"),
    Pair("x", "XXXXX") ## make more obvious
short_labels = [replace(s, replacements...) for s in str_labels]

GraphPlot.gplot( ## existing plotting functionality
    meta_graph,
    layout=spring_layout,
    ## nodelabel=collect(labels(meta_graph)),
    nodelabel=short_labels,
    nodelabelsize=0.1,
    NODESIZE=0.02, ## diameter of the nodes
    NODELABELSIZE=1.5,
    # nodelabelc="white",
    nodelabelc="green",
    nodelabeldist=0.0,
    nodefillc=nothing, ## "cyan"
    edgestrokec="red",
    ## ImageSize = (800, 800) ##- does not work
)
s₀_1 cv_2 cv_3 Nmc_4 u_5 cv_6 cv_7 Nmc_8 hIuI_9 #70_10 gIsI_11 #69_12 ghSum_13 +_14 s_15 cv_16 Nmp_17 XXXXX_18 cv_19 Nmc_20 cv_21 cv_22 Nmc_23 u_24 cv_25 cv_26 Nmc_27 hIuI_28 #70_29 gIsI_30 #69_31 ghSum_32 +_33 s_34 cv_35 Nmp_36 XXXXX_37 cv_38 Nmc_39 cv_40 cv_41 Nmc_42 u_43 cv_44 cv_45 Nmc_46 hIuI_47 #70_48 gIsI_49 #69_50 ghSum_51 +_52 s_53 cv_54 Nmp_55 XXXXX_56 cv_57 Nmc_58 cv_59 cv_60 Nmc_61 
_cntrl_hash_prefix = "70"
_state_hash_prefix = "69"
"69"
_san_node_ids = Dict(
    ## "#45_14" => "HX_14",
    "#$(_cntrl_hash_prefix)_14" => "CH_14",
    ## "#45_39" => "HX_39",
    "#$(_cntrl_hash_prefix)_39" => "CH_39",
    ## "#45_64" => "HX_64",
    "#$(_cntrl_hash_prefix)_64" => "CH_64",

    "#$(_state_hash_prefix)_16" => "SH_16",
    "#$(_state_hash_prefix)_41" => "SH_41",
    "#$(_state_hash_prefix)_66" => "SH_66",
)
_lup_enabled = true
_lup_agc = false ## append global_counter values
_raw_links_enabled = false
_CFFG = false
_deltas_enabled = false
_tsep = 6*_xsep
_lup_dict = Dict{String, String}(
    "constvar_8_9" => _lup_agc ? "\\mathbf{m}_u\\_8\\_9" : "\\mathbf{m}_u", 
    "constvar_10_11" => _lup_agc ? "\\mathbf{V}_u\\_10\\_11" : "\\mathbf{V}_u", 
    "MvNormalMeanCovariance_12" => _lup_agc ? "\\mathcal{N}\\_12" : "\\mathcal{N}",
    "u_7" => _lup_agc ? "\\mathbf{u}_1\\_7" : "\\mathbf{u}_1",

    "constvar_33_34" => _lup_agc ? "\\mathbf{m}_u\\_33\\_34" : "\\mathbf{m}_u", 
    "constvar_35_36" => _lup_agc ? "\\mathbf{V}_u\\_35\\_36" : "\\mathbf{V}_u", 
    "MvNormalMeanCovariance_37" => _lup_agc ? "\\mathcal{N}\\_37" : "\\mathcal{N}",
    "u_32" => _lup_agc ? "\\mathbf{u}_2\\_32" : "\\mathbf{u}_2",

    "constvar_58_59" => _lup_agc ? "\\mathbf{m}_u\\_58\\_59" : "\\mathbf{m}_u", 
    "constvar_60_61" => _lup_agc ? "\\mathbf{V}_u\\_60\\_61" : "\\mathbf{V}_u", 
    "MvNormalMeanCovariance_62" => _lup_agc ? "\\mathcal{N}\\_62" : "\\mathcal{N}",
    "u_57" => _lup_agc ? "\\mathbf{u}_3\\_57" : "\\mathbf{u}_3",

    "hIuI_13" => _lup_agc ? "h(\\mathbf{u})\\_13" : "h(\\mathbf{u})",
    "hIuI_38" => _lup_agc ? "h(\\mathbf{u})\\_38" : "h(\\mathbf{u})",
    "hIuI_63" => _lup_agc ? "h(\\mathbf{u})\\_63" : "h(\\mathbf{u})",
    "+_18" => _lup_agc ? "+\\_18" : "+",
    "+_43" => _lup_agc ? "+\\_43" : "+",
    "+_68" => _lup_agc ? "+\\_68" : "+",

    "gIsI_15" => _lup_agc ? "g(\\mathbf{s})\\_15" : "g(\\mathbf{s})",
    "gIsI_40" => _lup_agc ? "g(\\mathbf{s})\\_40" : "g(\\mathbf{s})",
    "gIsI_65" => _lup_agc ? "g(\\mathbf{s})\\_65" : "g(\\mathbf{s})",
    "ghSum_17" => _lup_agc ? "g(\\mathbf{s}) {+} h(\\mathbf{u})\\_17" : "g(\\mathbf{s}) {+} h(\\mathbf{u})",
    "ghSum_42" => _lup_agc ? "g(\\mathbf{s}) {+} h(\\mathbf{u})\\_42" : "g(\\mathbf{s}) {+} h(\\mathbf{u})",
    "ghSum_67" => _lup_agc ? "g(\\mathbf{s}) {+} h(\\mathbf{u})\\_67" : "g(\\mathbf{s}) {+} h(\\mathbf{u})",
    # "ghSum_17" => _lup_agc ? "g {+} h" : "g {+} h", ## to compress space surrounding +; better than \\!+\\! & \\mathrel{+}
    # "ghSum_42" => _lup_agc ? "g {+} h" : "g {+} h",
    # "ghSum_67" => _lup_agc ? "g {+} h" : "g {+} h",

    "MvNormalMeanPrecision_22" => _lup_agc ? "\\mathcal{N}\\_22" : "\\mathcal{N}",
    "constvar_20_21" => _lup_agc ? "\\vartheta\\_20\\_21" : "\\vartheta", 
    "MvNormalMeanPrecision_47" => _lup_agc ? "\\mathcal{N}\\_47" : "\\mathcal{N}",
    "constvar_45_46" => _lup_agc ? "\\vartheta\\_45\\_46" : "\\vartheta", 
    "MvNormalMeanPrecision_72" => _lup_agc ? "\\mathcal{N}\\_72" : "\\mathcal{N}",
    "constvar_70_71" => _lup_agc ? "\\vartheta\\_70\\_71" : "\\vartheta", 

    "MvNormalMeanCovariance_6" => _lup_agc ? "\\mathcal{N}\\_6" : "\\mathcal{N}",
    "constvar_4_5" => _lup_agc ? "10^{12}\\_4\\_5" : "10^{12}",
    "constvar_2_3" => _lup_agc ? "\\mathbf{0}\\_2\\_3" : "\\mathbf{0}",
    "s₀_1" => _lup_agc ? "\\mathbf{s}_0\\_1" : "\\mathbf{s}_0", ## s₀ must always be in _lup_dict
    "s_19" => _lup_agc ? "\\mathbf{s}_1\\_19" : "\\mathbf{s}_1",
    "s_44" => _lup_agc ? "\\mathbf{s}_2\\_44" : "\\mathbf{s}_2",
    "s_69" => _lup_agc ? "\\mathbf{s}_3\\_69" : "\\mathbf{s}_3",
    
    "MvNormalMeanCovariance_26" => _lup_agc ? "\\mathcal{N}\\_26" : "\\mathcal{N}",
    "constvar_24_25" => _lup_agc ? "\\mathbf\\Theta\\_24\\_25" : "\\mathbf\\Theta",
    "MvNormalMeanCovariance_51" => _lup_agc ? "\\mathcal{N}\\_51" : "\\mathcal{N}",
    "constvar_49_50" => _lup_agc ? "\\mathbf\\Theta\\_49\\_50" : "\\mathbf\\Theta",
    "MvNormalMeanCovariance_76" =>  _lup_agc ? "\\mathcal{N}\\_76" : "\\mathcal{N}",
    "constvar_74_75" => _lup_agc ? "\\mathbf\\Theta\\_74\\_75" : "\\mathbf\\Theta",
    "x_23" => _lup_agc ? "\\mathbf{x}_1\\_23" : "\\mathbf{x}_1",
    "x_48" => _lup_agc ? "\\mathbf{x}_2\\_48" : "\\mathbf{x}_2",
    "x_73" => _lup_agc ? "\\mathbf{x}_3\\_73" : "\\mathbf{x}_3",

    "MvNormalMeanCovariance_31" => _lup_agc ? "\\mathcal{N}\\_31" : "\\mathcal{N}",
    "constvar_29_30" => _lup_agc ? "\\mathbf{V}_x\\_29\\_30" : "\\mathbf{V}_x", 
    "constvar_27_28" => _lup_agc ? "\\mathbf{m}_x\\_27\\_28" : "\\mathbf{m}_x", 
    "MvNormalMeanCovariance_56" => _lup_agc ? "\\mathcal{N}\\_56" : "\\mathcal{N}",
    "constvar_54_55" => _lup_agc ? "\\mathbf{V}_x\\_54\\_55" : "\\mathbf{V}_x", 
    "constvar_52_53" => _lup_agc ? "\\mathbf{m}_x\\_52\\_53" : "\\mathbf{m}_x", 
    "MvNormalMeanCovariance_81" => _lup_agc ? "\\mathcal{N}\\_81" : "\\mathcal{N}",
    "constvar_79_80" => _lup_agc ? "\\mathbf{V}_x\\_79\\_80" : "\\mathbf{V}_x", 
    "constvar_77_78" => _lup_agc ? "\\mathbf{m}_x\\_77\\_78" : "\\mathbf{m}_x", 

    # "MvNormalMeanCovariance_36" => _lup_agc ? "\\bm{\\mathcal{N}}\\_36" : "",
    # "MvNormalMeanPrecision_76" => _lup_agc ? "\\bm{\\mathcal{N}}\\_76" : "",
    # "MvNormalMeanCovariance_31" =>  _lup_agc ? "\\bm{\\mathcal{N}}\\_31" : "",
);
drone_kobus_graph = create_kobus_graph(drone_gppl_model)
1: s₀_1
2: constvar_2
3: constvar_3
4: MvNormalMeanCovariance_4
5: u_5
6: constvar_6
7: constvar_7
8: MvNormalMeanCovariance_8
9: hIuI_9
10: #70_10
11: gIsI_11
12: #69_12
13: ghSum_13
14: +_14
15: s_15
16: constvar_16
17: MvNormalMeanPrecision_17
18: x_18
19: constvar_19
20: MvNormalMeanCovariance_20
21: constvar_21
22: constvar_22
23: MvNormalMeanCovariance_23
24: u_24
25: constvar_25
26: constvar_26
27: MvNormalMeanCovariance_27
28: hIuI_28
29: #70_29
30: gIsI_30
31: #69_31
32: ghSum_32
33: +_33
34: s_34
35: constvar_35
36: MvNormalMeanPrecision_36
37: x_37
38: constvar_38
39: MvNormalMeanCovariance_39
40: constvar_40
41: constvar_41
42: MvNormalMeanCovariance_42
43: u_43
44: constvar_44
45: constvar_45
46: MvNormalMeanCovariance_46
47: hIuI_47
48: #70_48
49: gIsI_49
50: #69_50
51: ghSum_51
52: +_52
53: s_53
54: constvar_54
55: MvNormalMeanPrecision_55
56: x_56
57: constvar_57
58: MvNormalMeanCovariance_58
59: constvar_59
60: constvar_60
61: MvNormalMeanCovariance_61
Dict{String, Dict{Symbol, Any}} with 61 entries:
  "MvNormalMeanCovariance_4"  => Dict(:node_lup=>"MvNormalMeanCovariance\\_4", …
  "MvNormalMeanCovariance_8"  => Dict(:node_lup=>"MvNormalMeanCovariance\\_8", …
  "MvNormalMeanCovariance_61" => Dict(:node_lup=>"MvNormalMeanCovariance\\_61",…
  "#70_10"                    => Dict(:node_lup=>"#70\\_10", :node_code=>10, :n…
  "#70_29"                    => Dict(:node_lup=>"#70\\_29", :node_code=>29, :n…
  "+_14"                      => Dict(:node_lup=>"+\\_14", :node_code=>14, :nod…
  "s_34"                      => Dict(:node_lup=>"s\\_34", :node_code=>34, :nod…
  "constvar_38"               => Dict(:node_lup=>"constvar\\_38", :node_code=>3…
  "s_15"                      => Dict(:node_lup=>"s\\_15", :node_code=>15, :nod…
  "constvar_21"               => Dict(:node_lup=>"constvar\\_21", :node_code=>2…
  "hIuI_47"                   => Dict(:node_lup=>"hIuI\\_47", :node_code=>47, :…
  "MvNormalMeanPrecision_17"  => Dict(:node_lup=>"MvNormalMeanPrecision\\_17", …
  "+_52"                      => Dict(:node_lup=>"+\\_52", :node_code=>52, :nod…
  "x_56"                      => Dict(:node_lup=>"x\\_56", :node_code=>56, :nod…
  "constvar_57"               => Dict(:node_lup=>"constvar\\_57", :node_code=>5…
  "s_53"                      => Dict(:node_lup=>"s\\_53", :node_code=>53, :nod…
  "MvNormalMeanCovariance_23" => Dict(:node_lup=>"MvNormalMeanCovariance\\_23",…
  "ghSum_13"                  => Dict(:node_lup=>"ghSum\\_13", :node_code=>13, …
  "MvNormalMeanPrecision_36"  => Dict(:node_lup=>"MvNormalMeanPrecision\\_36", …
  ⋮                           => ⋮
layers = Dict()
kobus_graph = drone_kobus_graph
gppl_model = drone_gppl_model

cntrl_sysvars, cntrl_facs, cntrl_parvars = find_cntrl_nodes(kobus_graph, "u", gppl_model)
layers["cntrl"] = Layer(
    name="cntrl",
    anchor=3*_xsep,
    iniparvars = [],
    inifac = "",
    parvar="",
    ## sysvars = ["u_7", "u_32", "u_57"],
    sysvars = cntrl_sysvars,
    ## facs = [["MvNormalMeanCovariance_12", "MvNormalMeanCovariance_37", "MvNormalMeanCovariance_62"]],
    facs = [cntrl_facs],
    ## parvars = [
    #     [
    #         ["constvar_8_9", "constvar_10_11"], 
    #         ["constvar_33_34", "constvar_35_36"],
    #         ["constvar_58_59", "constvar_60_61"]
    #     ] ]
    parvars = [cntrl_parvars]    
    )
## layers["cntrl"].links = layers["stateU2"].facs

layers["stateU1"] = Layer(
    name="stateU1",
    anchor=3*_xsep,
    iniparvars = [],
    inifac = "",
    parvar="",
    sysvars = ["hIuI_13", "hIuI_38", "hIuI_63"],
    facs = [
        [
            sanitized("#$(_cntrl_hash_prefix)_14"), 
            sanitized("#$(_cntrl_hash_prefix)_39"), 
            sanitized("#$(_cntrl_hash_prefix)_64")
        ]
    ],    
    parvars = [],
    link_fac_side="north", link_fac_offset=_divn[1][1])
## layers2["stateU1"].links = [layers2["state"].facs[3]]

state_inifac = find_inifac_node(kobus_graph, "s₀", gppl_model)
## _, state0_sysvars, state0_facs = find_state0_nodes(kobus_graph, state0_inifac, gppl_model)
_, state_vars, state_facs = find_state_nodes(kobus_graph, state_inifac, gppl_model)
state_out_var, state_iniparvars = find_iniparvar_nodes(kobus_graph, state_inifac, gppl_model)
layers["state"] = Layer(
    name="state",
    anchor=1*_xsep,
    ## iniparvars = ["constvar_2_3", "constvar_4_5"],
    iniparvars = state_iniparvars,
    ## inifac = "MvNormalMeanCovariance_6",
    inifac = state_inifac,
    parvar="",
    sysvars = ["s₀_1", "s_19", "s_44", "s_69"], ## state0_vars only return first seq, so hand-configure
    facs = [ ## in general, facs could be interleafed with vars; first & last layers must be facs
        [sanitized("#$(_state_hash_prefix)_16"), sanitized("#$(_state_hash_prefix)_41"), sanitized("#$(_state_hash_prefix)_66")],
        ["gIsI_15", "gIsI_40", "gIsI_65"], 
        ["+_18", "+_43", "+_68"],
        ["ghSum_17", "ghSum_42", "ghSum_67"], 
        ["MvNormalMeanPrecision_22", "MvNormalMeanPrecision_47", "MvNormalMeanPrecision_72"]
    ],
    parvars = [
        [],
        [], 
        [],
        [], 
        [["constvar_20_21"], ["constvar_45_46"], ["constvar_70_71"]] ])
## layers["state"].links = layers["obser"].facs

obser_sysvars, obser_facs, obser_parvars = find_obser_nodes(kobus_graph, "x", gppl_model)
layers["obser"] = Layer(
    name="obser",
    anchor=6*_xsep,
    iniparvars = [],
    inifac = "",
    parvar="",
    ## sysvars = ["x_23", "x_48", "x_73"], 
    sysvars = obser_sysvars,
    ## facs = [["MvNormalMeanCovariance_26", "MvNormalMeanCovariance_51", "MvNormalMeanCovariance_76"]],
    facs = [obser_facs],
    ## parvars = [
    #     [["constvar_24_25"], ["constvar_49_50"], ["constvar_74_75"]] ]
    parvars = [obser_parvars]
    )
## layers["obser"].links = layers["prefr"].facs

prefr_facs, prefr_parvars = find_prefr_nodes(kobus_graph, "x", gppl_model)
layers["prefr"] = Layer(
    name="prefr",
    anchor=6*_xsep,
    iniparvars = [],
    inifac = "",
    parvar="",
    sysvars = [], 
    ## facs = [["MvNormalMeanCovariance_31", "MvNormalMeanCovariance_56", "MvNormalMeanCovariance_81"]],
    facs = [prefr_facs],
    ## parvars = [
    #     [
    #         ["constvar_27_28", "constvar_29_30"], 
    #         ["constvar_52_53", "constvar_54_55"], 
    #         ["constvar_77_78", "constvar_79_80"]
    #     ] ]
    parvars = [prefr_parvars]
    )

layers["cntrl"].links = layers["stateU1"].facs
layers["stateU1"].links = [layers["state"].facs[3]]
layers["state"].links = layers["obser"].facs
layers["obser"].links = layers["prefr"].facs
cntrl_var_label_strs: ["u_5", "u_24", "u_43"]
cntrl_var_label_strs[i]: u_5
    var_neighbours[j]: MvNormalMeanCovariance_8
        in_neighbor_label: constvar_6
        out_neighbor_label: u_5
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(u_5, out, NodeData in context  with properties name = u, index = 1), (constvar_6, μ, NodeData in context  with properties name = constvar, index = nothing), (constvar_7, Σ, NodeData in context  with properties name = constvar, index = nothing)]
    var_neighbours[j]: #70_10
        in_neighbor_label: u_5
        out_neighbor_label: hIuI_9
cntrl_var_label_strs[i]: u_24
    var_neighbours[j]: MvNormalMeanCovariance_27
        in_neighbor_label: constvar_25
        out_neighbor_label: u_24
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(u_24, out, NodeData in context  with properties name = u, index = 2), (constvar_25, μ, NodeData in context  with properties name = constvar, index = nothing), (constvar_26, Σ, NodeData in context  with properties name = constvar, index = nothing)]
    var_neighbours[j]: #70_29
        in_neighbor_label: u_24
        out_neighbor_label: hIuI_28
cntrl_var_label_strs[i]: u_43
    var_neighbours[j]: MvNormalMeanCovariance_46
        in_neighbor_label: constvar_44
        out_neighbor_label: u_43
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(u_43, out, NodeData in context  with properties name = u, index = 3), (constvar_44, μ, NodeData in context  with properties name = constvar, index = nothing), (constvar_45, Σ, NodeData in context  with properties name = constvar, index = nothing)]
    var_neighbours[j]: #70_48
        in_neighbor_label: u_43
        out_neighbor_label: hIuI_47
    var_neighbours[1]: MvNormalMeanCovariance_4
        in_neighbor_label: constvar_2
        out_neighbor_label: s₀_1
inifac_out_neighbor_str: s₀_1
inifac_str: MvNormalMeanCovariance_4
var: (s₀_1, out, NodeData in context  with properties name = s₀, index = nothing)
=== NEXT var node is: s₀_1
var_label: s₀_1
var_index: 0
var_neighbours: Any["MvNormalMeanCovariance_4", "#69_12"]
    var_neighbours[1]: MvNormalMeanCovariance_4
        in_neighbor_label: constvar_2
        out_neighbor_label: s₀_1
    var_neighbours[2]: #69_12
        in_neighbor_label: s₀_1
        out_neighbor_label: gIsI_11
        Take this var_label as next fac because we ONLY have 2 var_neighbours
        === NEXT fac node is: #69_12
var: (gIsI_11, out, NodeData in context  with properties name = gIsI, index = 1)
=== NEXT var node is: gIsI_11
var_label: gIsI_11
var_index: 1
var_neighbours: Any["#69_12", "+_14"]
    var_neighbours[1]: #69_12
        in_neighbor_label: s₀_1
        out_neighbor_label: gIsI_11
    var_neighbours[2]: +_14
        in_neighbor_label: gIsI_11
        out_neighbor_label: ghSum_13
        Take this var_label as next fac because we ONLY have 2 var_neighbours
        === NEXT fac node is: +_14
var: (ghSum_13, out, NodeData in context  with properties name = ghSum, index = 1)
=== NEXT var node is: ghSum_13
var_label: ghSum_13
var_index: 1
var_neighbours: Any["+_14", "MvNormalMeanPrecision_17"]
    var_neighbours[1]: +_14
        in_neighbor_label: gIsI_11
        out_neighbor_label: ghSum_13
    var_neighbours[2]: MvNormalMeanPrecision_17
        in_neighbor_label: ghSum_13
        out_neighbor_label: s_15
        Take this var_label as next fac because we ONLY have 2 var_neighbours
        === NEXT fac node is: MvNormalMeanPrecision_17
var: (s_15, out, NodeData in context  with properties name = s, index = 1)
=== NEXT var node is: s_15
var_label: s_15
var_index: 1
var_neighbours: Any["MvNormalMeanPrecision_17", "MvNormalMeanCovariance_20", "#69_31"]
    var_neighbours[1]: MvNormalMeanPrecision_17
        in_neighbor_label: ghSum_13
        out_neighbor_label: s_15
    var_neighbours[2]: MvNormalMeanCovariance_20
        in_neighbor_label: s_15
        out_neighbor_label: x_18
    var_neighbours[3]: #69_31
        in_neighbor_label: s_15
        out_neighbor_label: gIsI_30
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(s₀_1, out, NodeData in context  with properties name = s₀, index = nothing), (constvar_2, μ, NodeData in context  with properties name = constvar, index = nothing), (constvar_3, Σ, NodeData in context  with properties name = constvar, index = nothing)]
        out_var: s₀_1
obser_var_label_strs: ["x_18", "x_37", "x_56"]
obser_var_label_strs[i]: x_18
    var_neighbours[j]: MvNormalMeanCovariance_20
        in_neighbor_label: s_15
        out_neighbor_label: x_18
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(x_18, out, NodeData in context  with properties name = x, index = 1), (s_15, μ, NodeData in context  with properties name = s, index = 1), (constvar_19, Σ, NodeData in context  with properties name = constvar, index = nothing)]
    var_neighbours[j]: MvNormalMeanCovariance_23
        in_neighbor_label: constvar_21
        out_neighbor_label: x_18
obser_var_label_strs[i]: x_37
    var_neighbours[j]: MvNormalMeanCovariance_39
        in_neighbor_label: s_34
        out_neighbor_label: x_37
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(x_37, out, NodeData in context  with properties name = x, index = 2), (s_34, μ, NodeData in context  with properties name = s, index = 2), (constvar_38, Σ, NodeData in context  with properties name = constvar, index = nothing)]
    var_neighbours[j]: MvNormalMeanCovariance_42
        in_neighbor_label: constvar_40
        out_neighbor_label: x_37
obser_var_label_strs[i]: x_56
    var_neighbours[j]: MvNormalMeanCovariance_58
        in_neighbor_label: s_53
        out_neighbor_label: x_56
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(x_56, out, NodeData in context  with properties name = x, index = 3), (s_53, μ, NodeData in context  with properties name = s, index = 3), (constvar_57, Σ, NodeData in context  with properties name = constvar, index = nothing)]
    var_neighbours[j]: MvNormalMeanCovariance_61
        in_neighbor_label: constvar_59
        out_neighbor_label: x_56
obser_var_label_strs: ["x_18", "x_37", "x_56"]
obser_var_label_strs[i]: x_18
    var_neighbours[j]: MvNormalMeanCovariance_20
        in_neighbor_label: s_15
        out_neighbor_label: x_18
    var_neighbours[j]: MvNormalMeanCovariance_23
        in_neighbor_label: constvar_21
        out_neighbor_label: x_18
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(x_18, out, NodeData in context  with properties name = x, index = 1), (constvar_21, μ, NodeData in context  with properties name = constvar, index = nothing), (constvar_22, Σ, NodeData in context  with properties name = constvar, index = nothing)]
obser_var_label_strs[i]: x_37
    var_neighbours[j]: MvNormalMeanCovariance_39
        in_neighbor_label: s_34
        out_neighbor_label: x_37
    var_neighbours[j]: MvNormalMeanCovariance_42
        in_neighbor_label: constvar_40
        out_neighbor_label: x_37
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(x_37, out, NodeData in context  with properties name = x, index = 2), (constvar_40, μ, NodeData in context  with properties name = constvar, index = nothing), (constvar_41, Σ, NodeData in context  with properties name = constvar, index = nothing)]
obser_var_label_strs[i]: x_56
    var_neighbours[j]: MvNormalMeanCovariance_58
        in_neighbor_label: s_53
        out_neighbor_label: x_56
    var_neighbours[j]: MvNormalMeanCovariance_61
        in_neighbor_label: constvar_59
        out_neighbor_label: x_56
neighbors: Tuple{GraphPPL.NodeLabel, GraphPPL.EdgeLabel, GraphPPL.NodeData}[(x_56, out, NodeData in context  with properties name = x, index = 3), (constvar_59, μ, NodeData in context  with properties name = constvar, index = nothing), (constvar_60, Σ, NodeData in context  with properties name = constvar, index = nothing)]
1-element Vector{Vector{String}}:
 ["MvNormalMeanCovariance_23", "MvNormalMeanCovariance_42", "MvNormalMeanCovariance_61"]
drone_tikz = generate_tikz(
    heading = "Drone Flying to Target",
    Model = drone_gppl_model,
    layers = layers) ## fac layers, i.e. not var (vars are handled with assoced layer)
## print(drone_tikz)
show_tikz(drone_tikz); ## write to file rather than show in notebook
\documentclass{standalone}
\usepackage{tikz}
\usetikzlibrary{calc}
\usepackage{amsmath}
\usepackage{bm}
\newcommand\ns{10mm} %% define node (minimum) size
\newcommand\nr{.5*\ns} %% define node radius
\newcommand\ar{.4*\ns} %% define arc radius
\newcommand\br{.19*\ns} %% define bead radius
\DeclareUnicodeCharacter{03B8}{\theta}
\DeclareUnicodeCharacter{03B3}{\gamma}
\DeclareUnicodeCharacter{209C}{\0}
\DeclareUnicodeCharacter{208B}{\0}
\DeclareUnicodeCharacter{2081}{\0}
\begin{document}
\begin{tikzpicture}[
    scale=1, 
    draw=black, 
    fac/.style={rectangle, draw=black!100, fill=black!0, minimum size=\ns, inner sep=0pt},
    dlt/.style={rectangle, draw=black!100, fill=black!100, minimum size=.2*\ns},
    var/.style={circle, draw=gray!100, fill=gray!100, minimum size=1mm, inner sep=0pt},
    bus/.style={rectangle, draw=black!100, fill=black!100, minimum width=30mm, minimum height=.5mm, inner sep=0pt},
    eql/.style={rectangle, draw=black!100, fill=black!0, minimum size=.5*\ns},
    text=black,
    ]
\node[font=\Large] at (4,29) {\textbf{Drone Flying to Target}};
⋮
⋮
⋮
\path[use as bounding box] ([xshift=-1cm, yshift=-1cm]current bounding box.south west) rectangle ([xshift=1cm, yshift=2cm]current bounding box.north east);
\end{tikzpicture}
\end{document}
seq: ["SH_16", "gIsI_15", "+_18", "ghSum_17", "MvNormalMeanPrecision_22"]
seq: ["SH_41", "gIsI_40", "+_43", "ghSum_42", "MvNormalMeanPrecision_47"]
seq: ["SH_66", "gIsI_65", "+_68", "ghSum_67", "MvNormalMeanPrecision_72"]

Drone flying to target