-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathdeeper_tutorial.jl
More file actions
2920 lines (2430 loc) · 115 KB
/
Copy pathdeeper_tutorial.jl
File metadata and controls
2920 lines (2430 loc) · 115 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
### A Pluto.jl notebook ###
# v1.0.3
using Markdown
using InteractiveUtils
# ╔═╡ 51fb74cd-974f-5403-8f48-8a7f9e05d363
begin
using LinearAlgebra
using Statistics
using Printf
using PlutoUI
using Plots
using Random
using ReactiveMP
using RxInfer
using StableRNGs
using OneHotArrays
default(; linewidth = 2, markersize = 5, legend = :best)
# Render numerical statements from this run rather than copying stored outputs.
fmt(x) = Printf.format(Printf.Format("%.5g"), x)
fmt_range(xs) = join(fmt.(extrema(xs)), "–")
nothing
end
# ╔═╡ 46d5396c-90ab-5503-8c84-6d353eba8e5e
md"""
# Solving Linear Systems with Message Passing
"""
# ╔═╡ 65ceb5c7-4eb5-517e-9011-f4ee7f3d46d5
begin
TableOfContents()
end
# ╔═╡ 1dcd6be4-b49d-57a5-959d-4c6bb7039b68
md"""
## A linear system as a Gaussian model
Assume that ``A`` is symmetric and positive definite. The unnormalized Gaussian
```math
p(x) \propto \exp\left(-\frac{1}{2}x^\mathsf{T}Ax+b^\mathsf{T}x\right)
```
has mean ``A^{-1}b``, so finding its mean solves ``Ax=b``. Expanding the exponent gives two kinds of factor:
- a univariate Normal distribution for each diagonal entry ``\exp(b_ix_i - \frac{1}{2} x_i^2A_{ii})``;
- an edge potential ``\exp(-A_{ij}x_i x_j)`` for each nonzero off-diagonal entry.
The `GaussianCoupling` node represents exactly that edge potential: `GaussianCoupling(xᵢ, -A[i, j])`.
"""
# ╔═╡ de5bfe1c-912a-5c5b-8b26-663e60faf1f3
@model function linear_system_model(b, A)
n = length(b)
# exp(bᵢxᵢ - Aᵢᵢxᵢ²/2), written as an ordinary normalized Gaussian.
for i in 1:n
x[i] ~ Normal(mean = b[i] / A[i, i], precision = A[i, i])
end
# One GaussianCoupling factor for every edge of the sparse matrix.
for i in 1:(n - 1), j in (i + 1):n
if !iszero(A[i, j])
x[j] ~ GaussianCoupling(x[i], -A[i, j])
end
end
end
# ╔═╡ 49438deb-5e38-5d41-85cb-1444b345d76d
md"""
The graph can contain feedback loops, so it needs an initial message before updates can circulate. We start with a broad Gaussian message. This is only a computational starting point; it is not an extra prior in the model. `KeepEach()` records the belief after every sweep so that we can watch the algorithm settle.
"""
# ╔═╡ 391c4c3c-45a4-5073-9565-0ff0c8ee0406
function solve_with_messages(A, b; iterations = 20)
@assert issymmetric(A) "A must be symmetric"
@assert isposdef(A) "A must be positive definite"
initialization = @initialization begin
μ(x) = NormalMeanVariance(0.0, 1e6)
end
result = infer(
model = linear_system_model(A = A),
data = (b = b,),
initialization = initialization,
returnvars = (x = KeepEach(),),
iterations = iterations,
options = (limit_stack_depth = 100,), # keeps large graphs from exhausting the stack
)
beliefs = result.posteriors[:x]
means = reduce(hcat, (mean.(belief) for belief in beliefs))
variances = reduce(hcat, (var.(belief) for belief in beliefs))
return (; result, beliefs, means, variances)
end
# ╔═╡ fd4429d7-2e6a-5231-afd5-acd06b12f372
md"""
## Warm-up: a chain
First consider four unknowns coupled in a chain. The corresponding factor graph has no cycles. This is the friendly case: both the means and marginal variances returned by message passing are exact once the messages have converged.
"""
# ╔═╡ 90783867-f800-5880-b4d8-fb77c0f12369
begin
A_tree = [
3.0 -0.8 0.0 0.0
-0.8 3.5 -0.6 0.0
0.0 -0.6 2.8 -0.7
0.0 0.0 -0.7 2.5
]
b_tree = [1.0, -0.5, 2.0, 0.75]
tree = solve_with_messages(A_tree, b_tree; iterations = 12)
x_exact_tree = A_tree \ b_tree
v_exact_tree = diag(inv(A_tree))
println("largest mean error: ", maximum(abs.(tree.means[:, end] - x_exact_tree)))
println("largest variance error: ", maximum(abs.(tree.variances[:, end] - v_exact_tree)))
end
# ╔═╡ 7db65f1c-5b60-592a-aa60-fbae23aaafb0
md"""
Alongside the numbers, look at the graph itself. Circles are unknowns, the small dark squares are the Gaussian self-potentials, and every `A` square is one `GaussianCoupling` factor sitting on a nonzero off-diagonal entry of ``A``. For the chain there is exactly one path between any two unknowns: no loops, so converged message passing gives exact means and marginal variances. The picture at the bottom explains how to convert an A into the respective factor graph.
"""
# ╔═╡ 3a9ba1f8-6e9b-5f99-a5d5-cb46b5b66074
begin
# A picture of the factor graph itself, drawn directly with Plots so it composes with
# other subplots: a circle per unknown, a labelled square per `GaussianCoupling` factor (one for
# every nonzero off-diagonal of A), and a small dark square for each Gaussian
# self-potential. The double loop mirrors `linear_system_model` line by line.
function factor_graph(A; positions, kwargs...)
n = size(A, 1)
sub(k) = join('₀' + d for d in reverse(digits(k)))
selfpos = [p .+ (0.16, 0.16) for p in positions]
xs = [first.(positions); first.(selfpos)]
ys = [last.(positions); last.(selfpos)]
p = plot(;
legend = false, axis = false, grid = false, ticks = false, aspect_ratio = :equal,
xlims = extrema(xs) .+ (-0.5, 0.5), ylims = extrema(ys) .+ (-0.5, 0.5), kwargs...
)
# Edges first, so that the nodes are drawn on top of them.
for i in 1:n
plot!(p, [positions[i][1], selfpos[i][1]], [positions[i][2], selfpos[i][2]]; color = :gray40, linewidth = 1.5)
end
for i in 1:(n - 1), j in (i + 1):n
if !iszero(A[i, j])
plot!(p, [positions[i][1], positions[j][1]], [positions[i][2], positions[j][2]]; color = :gray40, linewidth = 1.5)
end
end
# A labelled square per `GaussianCoupling` factor, halfway along its edge.
for i in 1:(n - 1), j in (i + 1):n
if !iszero(A[i, j])
mid = (positions[i] .+ positions[j]) ./ 2
scatter!(p, [mid[1]], [mid[2]]; marker = :square, markersize = 13, color = :gray90, markerstrokecolor = :gray40)
annotate!(p, mid[1], mid[2], text("A$(sub(i))$(sub(j))", 8, :gray10))
end
end
# A small dark square per Gaussian self-potential and a circle per unknown.
scatter!(p, first.(selfpos), last.(selfpos); marker = :square, markersize = 6, color = :gray60, markerstrokecolor = :gray40)
scatter!(p, first.(positions), last.(positions); marker = :circle, markersize = 12, color = :white, markerstrokecolor = :black)
for i in 1:n
annotate!(p, positions[i][1], positions[i][2], text("x$(sub(i))", 9, :gray10))
end
return p
end
fg = factor_graph(A_tree; positions = [(float(i), 0.0) for i in 1:4], title = "Chain: factor graph");
p_tree_graph = heatmap(
abs.(A_tree .- Diagonal(diag(A_tree))),
title = "Chain: edge strengths |Aᵢⱼ|",
xlabel = "variable j", ylabel = "variable i",
aspect_ratio = :equal, color = :blues, yflip = true, xticks = 1:4, yticks = 1:4,
);
plot(p_tree_graph, fg; layout = (1, 2), size = (900, 360))
end
# ╔═╡ 1b6cacad-1d82-5b00-8ae6-6094b824ad67
md"""
## Add one edge: a loopy graph
Closing the chain into a ring creates a feedback loop. We keep the matrix strictly diagonally dominant, a convenient sufficient condition for convergence. Watch what changes: the means still converge to `A \ b`, while the variances become loop-dependent approximations.
"""
# ╔═╡ f6a7cafc-58e2-534e-99dc-609b2bffec51
begin
A_loop = copy(A_tree)
A_loop[1, 4] = A_loop[4, 1] = -0.9
loop = solve_with_messages(A_loop, b_tree; iterations = 20)
x_exact_loop = A_loop \ b_tree
v_exact_loop = diag(inv(A_loop))
residuals = [norm(A_loop * loop.means[:, k] - b_tree) for k in axes(loop.means, 2)]
@assert all(A_loop[i, i] > sum(abs, A_loop[i, :]) - A_loop[i, i] for i in axes(A_loop, 1))
println("final residual ‖Ax-b‖: ", residuals[end])
end
# ╔═╡ 625ed0a9-7296-5590-b9c3-b654b0a069df
md"""
Here is what that one extra matrix entry does to the graph: it closes the chain into a ring. Now there are two routes between any pair of unknowns and messages can circulate around the cycle indefinitely - this is exactly why the variances come out approximate below, while the means still converge to the true solution.
"""
# ╔═╡ ebd306ca-2028-59bf-9f3f-7c03449356ee
begin
factor_graph(A_loop; positions = [(0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0)])
end
# ╔═╡ d1db6159-e47a-5204-a116-ff479c61240c
begin
let
p_convergence = plot(
residuals, yscale = :log10, marker = :circle, label = false,
xlabel = "message-passing sweep", ylabel = "‖Ax-b‖₂",
title = "Local messages converge globally",
)
p_paths = plot(
xlabel = "message-passing sweep", ylabel = "mean",
title = "Belief means settling",
)
for i in axes(loop.means, 1)
plot!(p_paths, loop.means[i, :], label = "x$i", marker = :circle)
hline!(p_paths, [x_exact_loop[i]], color = i, linestyle = :dash, label = false)
end
plot(p_convergence, p_paths; layout = (1, 2), size = (950, 360))
end
end
# ╔═╡ 5d933fec-990e-54cb-8dc9-317369035dfe
begin
let
variance_error = 100 .* (loop.variances[:, end] .- v_exact_loop) ./ v_exact_loop
p_variances = bar(
(1:4) .- 0.18, v_exact_loop,
label = "diag(inv(A))", bar_width = 0.35,
xlabel = "variable", ylabel = "marginal variance",
title = "Cycles change the variance estimate", xticks = 1:4,
)
bar!(p_variances, (1:4) .+ 0.18, loop.variances[:, end], label = "loopy message passing", bar_width = 0.35)
p_variance_error = bar(
1:4, variance_error, label = false,
xlabel = "variable", ylabel = "relative error (%)",
title = "Variance error despite exact means",
)
plot(p_variances, p_variance_error; layout = (1, 2), size = (950, 360))
end
end
# ╔═╡ 4606b4bf-17a6-5ede-bbfb-96ecd24c8130
md"""
## Turn up the coupling
Message passing is most comfortable when each diagonal term dominates the edges around it. The next experiment scales every off-diagonal coupling while keeping the diagonal fixed. As stronger couplings reduce the diagonal-dominance margin, information echoes more strongly around the ring and convergence slows.
"""
# ╔═╡ 74db928c-ba0c-5048-a023-e1ca1ebbe38f
begin
let
diagonal = Diagonal(diag(A_loop))
off_diagonal = A_loop - diagonal
coupling_scales = [0.25, 0.6, 0.9]
p_coupling = plot(
yscale = :log10, xlabel = "message-passing sweep", ylabel = "‖Ax-b‖₂",
title = "Stronger feedback needs more sweeps",
)
for scale in coupling_scales
A_scaled = diagonal + scale * off_diagonal
run = solve_with_messages(Matrix(A_scaled), b_tree; iterations = 20)
errors = [norm(A_scaled * run.means[:, k] - b_tree) for k in axes(run.means, 2)]
plot!(p_coupling, errors, marker = :circle, label = "scale = $scale")
end
p_coupling
end
end
# ╔═╡ 9760db36-0009-5ee1-8dee-720fde55562c
md"""
## From one solve to a stream of solves: heat on a grid
So far ``b`` was a fixed vector. Now let ``A`` describe a ``20 \times 20`` grid of cells (the discrete Laplacian plus a small leak ``\kappa`` to the ambient), and let the right-hand side change over time: four segments of the boundary are in contact with the outside, with alternating hot and cold inputs, and those inputs drift. Each time step is one linear system, and the graph we already built solves all of them at once, while borrowing strength across time.
The only sensors are a thermometer at every source segment, and they occasionally go offline. When a source sensor drops out, its temperature has to be inferred from its observations before and after the gap and from what the other three sources are doing. This is batch smoothing, with future observations available, rather than online filtering.
The source values are additive right-hand-side inputs, shared across five cells per segment, rather than clamped boundary temperatures. The plots retain the original example’s temperature labels; this is an illustrative steady-state model, not a calibrated thermal simulation. There is no heat storage in the field.
"""
# ╔═╡ cc7c8129-b7d3-5570-a4de-2b653d6d8a6c
begin
function grid_laplacian(m, κ)
idx(r, c) = (r - 1) * m + c
A = zeros(m * m, m * m)
for r in 1:m, c in 1:m, (dr, dc) in ((1, 0), (-1, 0), (0, 1), (0, -1))
rr, cc = r + dr, c + dc
if 1 <= rr <= m && 1 <= cc <= m
A[idx(r, c), idx(rr, cc)] = -1.0
A[idx(r, c), idx(r, c)] += 1.0
end
end
return A + κ * I
end
m, κ = 20, 0.1 # stronger leakage improves diagonal dominance; sweep counts also depend on the stopping tolerance
cell(r, c) = (r - 1) * m + c
A_grid = grid_laplacian(m, κ)
n_grid = size(A_grid, 1)
source_cells = [ # north, east, south, west
[cell(1, c) for c in 6:10], [cell(r, m) for r in 6:10], [cell(m, c) for c in 11:15], [cell(r, 1) for r in 11:15],
]
K = length(source_cells)
cell_map = fill(0.0, m, m)
for i in reduce(vcat, source_cells)
cell_map[div(i - 1, m) + 1, mod(i - 1, m) + 1] = 1.0
end
p_grid_sparsity = heatmap(
abs.(A_grid .- Diagonal(diag(A_grid))), title = "Grid: edge strengths |Aᵢⱼ|",
xlabel = "cell j", ylabel = "cell i", aspect_ratio = :equal, color = :blues, yflip = true, colorbar = false,
)
p_grid_cells = heatmap(
cell_map, title = "Source segments", aspect_ratio = :equal, yflip = true,
color = cgrad([:white, :firebrick]), clims = (0, 1), colorbar = false, framestyle = :none,
)
plot(p_grid_sparsity, p_grid_cells; layout = (1, 2), size = (900, 400))
end
# ╔═╡ 619d031c-23d2-535c-8349-bf2a7a377659
md"""
The heat-grid factor graph has 400 unknowns, but its texture is already visible in a small patch. A ``4 \times 4`` corner of the grid shows what the chain and the ring could not: every plaquette of the lattice is a loop, and neighbouring loops share edges. This is the regime where the means stay exact under convergence while the variances are genuinely approximate.
"""
# ╔═╡ fb6db12c-70cd-5110-b74d-3dffdb7fb0aa
begin
factor_graph(grid_laplacian(4, κ); positions = [(float(c), -float(r)) for r in 1:4 for c in 1:4])
end
# ╔═╡ 6cdf2273-3d5b-5703-9804-4e10f6e9f3b4
md"""
The data are generated by a daily profile plus a slow drift that follows a random walk, one per source:
```math
b_k(t) = s_k(t) + d_k(t), \qquad d_k(t) = d_k(t-1) + \varepsilon_k(t), \quad \varepsilon_k(t) \sim \mathcal{N}(0, \tau_d^{-1}).
```
The true temperature field at each step is the steady state ``x_t = A^{-1} M b_t``, where ``M`` maps the four source values to their boundary cells. The model below is not told the profile ``s_k(t)``; it sees only the source sensors, which have precision ``\tau_y``, and two of those sensors go offline for a while: entries of `Y_b` become `missing`.
"""
# ╔═╡ c97e1c18-1742-5a94-8959-b26ea0cd03fd
begin
T, period = 120, 40
τ_d, τ_y = 25.0, 25.0 # drift, source sensors
rng = StableRNG(42)
profile_level, profile_swing, profile_phase = [2.0, -2.0, 1.0, -1.0], [8.0, 7.0, 8.0, 7.0], [0.0, π / 2, π, 3π / 2]
S_known = [profile_level[k] + profile_swing[k] * sin(2π * t / period + profile_phase[k]) for t in 1:T, k in 1:K]
D_true = cumsum(vcat(zeros(1, K), randn(rng, T - 1, K) ./ sqrt(τ_d)); dims = 1)
B_true = S_known .+ D_true
full_b(bk) = (v = zeros(n_grid); for k in 1:K; v[source_cells[k]] .= bk[k]; end; v)
X_true = reduce(hcat, (A_grid \ full_b(B_true[t, :]) for t in 1:T))' # T × n
dropouts = Dict(2 => 40:70, 4 => 90:100)
Y_b = Matrix{Union{Float64, Missing}}(B_true .+ randn(rng, T, K) ./ sqrt(τ_y))
for (k, window) in dropouts
Y_b[window, k] .= missing
end
source_names = ["north", "east", "south", "west"]
p_sources = plot(layout = (2, 2), size = (950, 500), legend = :topright)
for k in 1:K
plot!(p_sources[k], 1:T, S_known[:, k], linestyle = :dash, color = :gray, label = "daily profile sₖ(t)")
plot!(p_sources[k], 1:T, B_true[:, k], color = k, label = "true bₖ(t)")
scatter!(p_sources[k], 1:T, Y_b[:, k], color = k, markersize = 2, markerstrokewidth = 0, alpha = 0.6, label = "sensor")
haskey(dropouts, k) && vspan!(p_sources[k], [first(dropouts[k]), last(dropouts[k])], color = :gray, alpha = 0.2, label = "offline")
plot!(p_sources[k], title = "source: $(source_names[k])", xlabel = k > 2 ? "time" : "", ylabel = "temperature")
end
p_sources
end
# ╔═╡ 1e712a9a-a90e-5805-8867-fe6e046ab30d
md"""
### A state-space prior for the sources
Before wiring the sources into the grid, give them a prior of their own. The four temperatures at time ``t`` form one vector ``b_t`` that follows a random walk with an unknown ``4 \times 4`` transition precision ``\Lambda \sim \mathrm{Wishart}``: the sources drift together or apart, and ``\Lambda`` learns how. Each sensor is a scalar reading of one component, ``\langle b_t, e_k \rangle``, with a learned precision ``\tau``, so a single `missing` entry drops one likelihood term instead of a whole time slice.
We write this as a *submodel*. The arguments of a `@model` are its Markov blanket, the variables it shares with the outside: `y_b` (the data), `b` (the states it creates), `Λ` and `τ` (hyperparameters owned by the caller) and the constant `dim_latent`. Everything else stays internal. A submodel is invoked with `~`: all but one interface are given by keyword and the remaining one is bound to the left-hand side, so `b ~ lss_prior(y_b = y_b, Λ = Λ, τ = τ, dim_latent = 4)` hands the whole array of states back to the caller. Factorisation constraints for the factors inside a submodel live in a `for q in lss_prior` block. First we run it on its own, wrapped in a top-level model that declares ``\Lambda`` and ``\tau``.
"""
# ╔═╡ c35676e9-9053-5900-8505-f96af3ed17fd
begin
unit(k, n) = onehot(k, 1:n) # eₖ as a one-hot vector, so dot(b[t], eₖ) picks component k
@model function lss_prior(y_b, b, Λ, τ, dim_latent)
b[1] ~ MvNormalMeanPrecision(zeros(dim_latent), 0.01 * diageye(dim_latent))
for t in 1:size(y_b, 1)
if t > 1
b[t] ~ MvNormalMeanPrecision(b[t - 1], Λ)
end
for k in 1:dim_latent
y_b[t, k] ~ NormalMeanPrecision(dot(b[t], unit(k, dim_latent)), τ)
end
end
end
@model function linear_state_space_model(y_b, dim_latent)
Λ ~ Wishart(dim_latent + 2, diageye(dim_latent))
τ ~ GammaShapeRate(1, 1)
b ~ lss_prior(y_b = y_b, Λ = Λ, τ = τ, dim_latent = dim_latent)
end
lss_constraints = @constraints begin
for q in lss_prior
q(b, Λ, τ) = q(b)q(Λ)q(τ)
end
end
lss_init = @initialization begin
q(Λ) = Wishart(K + 2, diageye(K))
q(τ) = GammaShapeRate(1.0, 1.0)
end
lss_result = infer(
model = linear_state_space_model(dim_latent = K),
data = (y_b = Y_b,),
constraints = lss_constraints,
initialization = lss_init,
returnvars = (b = KeepLast(), Λ = KeepLast(), τ = KeepLast()),
options = (limit_stack_depth = 100,),
iterations = 10,
)
B_lss_mean = reduce(hcat, mean.(lss_result.posteriors[:b]))' # T × 4
B_lss_std = reduce(hcat, (sqrt.(var(q)) for q in lss_result.posteriors[:b]))'
observed = .!ismissing.(Y_b)
println("max error on observed entries: ", maximum(abs.(B_lss_mean[observed] .- B_true[observed])))
println("max error on missing entries: ", maximum(abs.(B_lss_mean[.!observed] .- B_true[.!observed])))
println("learned sensor precision: ", mean(lss_result.posteriors[:τ]), " (true τ_y = ", τ_y, ")")
end
# ╔═╡ e5c70a74-fd7c-56ed-8f83-d08cce125326
begin
let
p_lss = plot(layout = (2, 2), size = (1100, 760), legend = :outertop, legend_columns = 2, margin = 4Plots.mm)
for k in 1:K
haskey(dropouts, k) && vspan!(p_lss[k], [first(dropouts[k]), last(dropouts[k])], color = :gray, alpha = 0.2, label = "offline")
plot!(p_lss[k], 1:T, B_lss_mean[:, k], ribbon = 2 .* B_lss_std[:, k], color = k, fillalpha = 0.25, label = "posterior ± 2σ")
plot!(p_lss[k], 1:T, B_true[:, k], color = :black, linestyle = :dash, label = "true bₖ(t)")
plot!(p_lss[k], title = "source: $(source_names[k])", xlabel = k > 2 ? "time" : "", ylabel = "temperature")
end
p_lss
end
end
# ╔═╡ 2985f942-7c20-5bf1-9ef0-8055266e536e
md"""
While a sensor is offline the state-space prior carries the source using its observations before and after the gap and the learned coupling ``\Lambda`` to the other three sources. Uncertainty typically increases across these finite gaps. A random walk does not guarantee bounded uncertainty during an indefinitely long outage without informative observations. Now the same submodel goes into the grid.
"""
# ╔═╡ 69c6df8e-2715-55f3-9b5e-2a0b2d8b8654
begin
# Eliminate x analytically: J_b = J_LSS - β (I_T ⊗ M'A⁻¹M).
# This is a validation calculation, not part of the message-passing solver.
function gaussian_reference(A, source_cells, Y, Λ, τ, β)
steps, channels = size(Y)
n = size(A, 1)
M = zeros(n, channels)
for k in 1:channels
M[source_cells[k], k] .= 1.0
end
A_factor = cholesky(Symmetric(A))
H = A_factor \ M
S = Symmetric(M' * H)
J = zeros(steps * channels, steps * channels)
η = zeros(steps * channels)
block(t) = ((t - 1) * channels + 1):(t * channels)
J[block(1), block(1)] .+= 0.01 .* Matrix{Float64}(I, channels, channels)
for t in 1:steps
current = block(t)
J[current, current] .-= β .* S
if t > 1
previous = block(t - 1)
J[current, current] .+= Λ
J[previous, previous] .+= Λ
J[current, previous] .-= Λ
J[previous, current] .-= Λ
end
for k in 1:channels
if !ismissing(Y[t, k])
i = first(current) + k - 1
J[i, i] += τ
η[i] += τ * Y[t, k]
end
end
end
precision = Symmetric(J)
factor = cholesky(precision) # fails if the uncancelled model is improper
source_mean = permutedims(reshape(factor \ η, channels, steps))
source_covariance = factor \ Matrix{Float64}(I, steps * channels, steps * channels)
source_variance = permutedims(reshape(diag(source_covariance), channels, steps))
field_mean = source_mean * H'
conditional_variance = diag(A_factor \ Matrix{Float64}(I, n, n)) ./ β
field_variance = zeros(steps, n)
for t in 1:steps
C = source_covariance[block(t), block(t)]
field_variance[t, :] = conditional_variance + vec(sum((H * C) .* H; dims = 2))
end
return (; M, H, S, precision, η, source_mean, source_covariance,
source_variance, field_mean, field_variance, min_precision = eigmin(precision))
end
nothing
end
# ╔═╡ 51db8ee1-8fb7-54f1-ac41-d1040766f865
begin
β = 0.003 # original inverse temperature; the latent-source normaliser is retained
Λ_hat = mean(lss_result.posteriors[:Λ])
τ_hat = mean(lss_result.posteriors[:τ])
exact_grid = gaussian_reference(A_grid, source_cells, Y_b, Λ_hat, τ_hat, β)
longest_gap_source = first(sort(collect(keys(dropouts)); by = k -> length(dropouts[k]), rev = true))
longest_gap = length(dropouts[longest_gap_source])
gap_precision = 2 * Λ_hat[longest_gap_source, longest_gap_source] * (1 - cos(π / (longest_gap + 1)))
nothing
end
# ╔═╡ fcddcee4-39af-5d42-afa8-9a49512e1784
md"""
### The model
Conditional on its sources, each time slice is the linear system from the first half, scaled by an inverse temperature ``\beta``. Multiplying both ``A`` and ``b`` by ``\beta`` leaves the solution ``A^{-1}b`` untouched, but changes the covariance to ``(\beta A)^{-1}``: ``\beta`` says how tightly the field is expected to follow the steady-state equation. Interior cells keep the Gaussian self-potential of `linear_system_model` with ``b_i = 0``, written with a static right-hand side `b0` of zeros. The source cells are driven by the state instead: all cells of segment ``k`` share the temperature ``\langle b_t, e_k \rangle``, and the coupling ``\exp(\beta\, b_{t,k} x_i)`` is one `GaussianCoupling` node between that temperature and each cell of the segment. The prior on ``b_t`` is the submodel from the previous section, plugged in with one line. Its hyperparameters are not learned a second time: the first loopy sweeps through the grid send wild messages to the sources, and a variational update of ``\tau`` or ``\Lambda`` that sees them can become unstable. This example avoids those early updates by learning the hyperparameters in the standalone model first. The posterior means from the standalone run go in as constants instead, through the same two interfaces. The field itself has no dynamics of its own: the conditional mean of every slice is the steady state for its sources, and all memory lives in ``b_t``.
One subtlety of the potential remains. ``\exp(b^\mathsf{T}x - \frac{1}{2}x^\mathsf{T}Ax)`` integrates over ``x`` to ``\exp(\frac{1}{2}b^\mathsf{T}A^{-1}b)`` up to a constant. When ``b`` is data this is irrelevant; when ``b`` is latent it is a *negative* precision ``\beta\,M^\mathsf{T}A^{-1}M`` on the segment inputs (``M`` is the cell-to-segment indicator). On this grid the diagonal of ``M^\mathsf{T}A^{-1}M`` is **$(fmt(exact_grid.S[1, 1]))** per segment; it also has nonzero off-diagonal entries. Cancelling it exactly would need ``A^{-1}M``, four extra solves. We keep the original example's uncancelled potential and small ``\beta``, and check its precision explicitly below.
While a sensor is online its **learned** precision of **$(fmt(τ_hat))** outweighs the pull; **$(fmt(τ_y))** is the precision used to generate the observations. During a dropout the source relies on the temporal prior and the other sources. For a single missing channel with its neighbours and endpoints held fixed, the slowest random-walk mode across a gap of ``L`` entries has precision ``2\Lambda_{kk}(1-\cos(\pi/(L+1))) \approx \Lambda_{kk}\pi^2/(L+1)^2``. This scalar approximation is **$(fmt(gap_precision))** for the **$(longest_gap)**-entry outage. It is a diagnostic, not a guarantee for the coupled model: the full reduced source precision must be positive definite.
``\beta =`` **$(β)** gives a diagonal negative-precision contribution of **$(fmt(β * exact_grid.S[1, 1]))**, with largest eigenvalue **$(fmt(β * eigmax(exact_grid.S)))**. The reduced posterior precision is positive definite in this run (smallest eigenvalue **$(fmt(exact_grid.min_precision))**). Scaling by ``\beta`` leaves the conditional field mean ``A^{-1}Mb_t`` unchanged, but the uncancelled normaliser changes the posterior over latent sources. The effect during an outage is a bias, not new information; finite message-passing sweeps also introduce convergence error.
"""
# ╔═╡ de3604ce-2a5e-5cce-913c-86cd675327ce
@model function heat_grid_model(y_b, b0, A, β, Λ, τ, source_cells)
T, n, K = size(y_b, 1), size(A, 1), length(source_cells)
local x, bs # declared once so all loops below refer to the same variables
# The state-space prior on the sources, with the sensor likelihood inside it.
# Λ and τ are plain matrices/numbers here: the interfaces accept constants as well as variables.
b ~ lss_prior(y_b = y_b, Λ = Λ, τ = τ, dim_latent = K)
for t in 1:T
# The linear system at time t, exactly as before, scaled by β.
for i in 1:n
x[t, i] ~ Normal(mean = b0[i] / A[i, i], precision = β * A[i, i])
end
# Every cell of segment k is driven by component k of the state.
for k in 1:K
bs[t, k] := dot(b[t], unit(k, K))
for i in source_cells[k]
x[t, i] ~ GaussianCoupling(bs[t, k], β)
end
end
for i in 1:(n - 1), j in (i + 1):n
if !iszero(A[i, j])
x[t, j] ~ GaussianCoupling(x[t, i], -β * A[i, j])
end
end
end
end
# ╔═╡ 3eaeffed-1a08-5565-8955-b80eb9ded72a
md"""
The graph has loops within each time slice and a chain across time for the sources, so we again start from broad messages, on the field `x` and on the segment temperatures `bs`. With ``\Lambda`` and ``\tau`` fixed the whole model is Gaussian and no factorisation constraints are needed.
"""
# ╔═╡ 7ee9a9b6-cb79-58d0-9680-49bfb27ba2cd
begin
grid_init = @initialization begin
μ(x) = NormalMeanVariance(0.0, 1e6)
μ(bs) = NormalMeanVariance(0.0, 1e6)
end
grid_result = infer(
model = heat_grid_model(A = A_grid, β = β, Λ = Λ_hat, τ = τ_hat, source_cells = source_cells),
data = (y_b = Y_b, b0 = zeros(n_grid)),
initialization = grid_init,
returnvars = (x = KeepLast(), b = KeepLast()),
options = (limit_stack_depth = 100,), # the graph is deep: 120 slices × 400 cells
iterations = 10,
)
x_post = grid_result.posteriors[:x] # T × n matrix of beliefs
b_post = grid_result.posteriors[:b] # T beliefs over the 4-vector of source temperatures
X_mean, X_std = mean.(x_post), std.(x_post)
B_mean = reduce(hcat, mean.(b_post))' # T × 4
B_std = reduce(hcat, (sqrt.(var(q)) for q in b_post))' # T × 4
println("largest deviation of the posterior mean from the true field: ", maximum(abs.(X_mean - X_true)))
println("all variances proper: ", all(>(0), var.(x_post)) && all(q -> all(>(0), var(q)), b_post))
end
# ╔═╡ 6b4f64ed-954d-5d75-bc5a-9e7c8e29d364
md"""
### Sources recovered through the grid
"""
# ╔═╡ 265f8aca-588f-5823-8252-145569856248
begin
let
p_grid_sources = plot(layout = (2, 2), size = (1100, 760), legend = :outertop, legend_columns = 2, margin = 4Plots.mm)
for k in 1:K
haskey(dropouts, k) && vspan!(p_grid_sources[k], [first(dropouts[k]), last(dropouts[k])], color = :gray, alpha = 0.2, label = "offline")
plot!(p_grid_sources[k], 1:T, B_mean[:, k], ribbon = 2 .* B_std[:, k], color = k, fillalpha = 0.25, label = "through the grid ± 2σ")
plot!(p_grid_sources[k], 1:T, B_lss_mean[:, k], color = :gray30, linewidth = 1, label = "state-space model alone")
plot!(p_grid_sources[k], 1:T, B_true[:, k], color = :black, linestyle = :dash, label = "true bₖ(t)")
plot!(p_grid_sources[k], title = "source: $(source_names[k])", xlabel = k > 2 ? "time" : "", ylabel = "temperature")
end
p_grid_sources
end
end
# ╔═╡ 0920f1bc-7bcf-5bbf-a9e6-996adcb7dcf5
begin
reference_errors = (
sources = maximum(abs.(B_mean - exact_grid.source_mean)),
field = maximum(abs.(X_mean - exact_grid.field_mean)),
truth = maximum(abs.(X_mean - X_true)),
)
source_variance_ratio = B_std .^ 2 ./ exact_grid.source_variance
field_variance_ratio = X_std .^ 2 ./ exact_grid.field_variance
(; reference_errors, source_variance_ratio = extrema(source_variance_ratio),
field_variance_ratio = extrema(field_variance_ratio),
minimum_precision_eigenvalue = exact_grid.min_precision)
end
# ╔═╡ 0c5ec67b-19c0-59c1-bebd-817438aced7f
md"""
The estimates follow the state-space model alone closely. With no thermometers inside the plate the field has nothing of its own to say about the sources; the only route back from the grid is the un-cancelled normaliser, which has a small effect while the sensors are online and can push a blind source further away from zero during an outage. Any improvement in error against the simulated truth is incidental: this is a bias, not information. What the grid adds is the field: every slice is solved together with the sources, dropouts included.
### Are the numbers right?
The whole model is one big sparse Gaussian with ``120\times(400+4) =`` **$(T * (n_grid + K))** coordinates, counting the field and the four source states at each time. An exact reference does **not** require a dense matrix of that size: sparse factorisation is possible, and eliminating the field gives a **$(T * K) by $(T * K)** system for the sources. The check here forms that reduced precision, including the negative term above, and then recovers field means and marginal variances using ``A^{-1}M``. It does not form the inverse of the full joint precision.
On the earlier small examples, converged Gaussian message passing returned exact means and optimistic variances on the loopy graph. With the original **10** sweeps in the composed model, the maximum mean discrepancy from the exact Gaussian reference is **$(fmt(reference_errors.sources))** for the sources and **$(fmt(reference_errors.field))** for the field. The message-passing/exact marginal-variance ratio ranges from **$(fmt_range(source_variance_ratio))** for the sources and **$(fmt_range(field_variance_ratio))** for the field. These are numerical approximation checks, not errors against the simulated truth. In particular, ten sweeps do not establish convergence of this large graph; check the discrepancy from the exact reference before treating its mean as a solved linear system.
The maximum field-mean error against that truth is **$(fmt(reference_errors.truth))**. It also includes imperfect source recovery from noisy and missing data. Source uncertainty and errors typically grow during outages, but neither is restricted to times when a sensor is offline. The ribbons show approximate posterior uncertainty, not a guaranteed error bound.
"""
# ╔═╡ 5a0f1e28-b6d7-57ae-b800-635815aa800c
begin
let
p_field_error = plot(xlabel = "time", ylabel = "|posterior mean − truth|", title = "Source errors and sensor outages")
for k in 1:4
plot!(p_field_error, 1:T, abs.(B_mean[:, k] - B_true[:, k]), color = k, label = "source: $(source_names[k])")
end
for (k, window) in dropouts
vspan!(p_field_error, [first(window), last(window)], color = :gray, alpha = 0.2, label = k == 2 ? "sensor offline" : false)
end
p_field_error
end
end
# ╔═╡ da904129-d8e6-5c7b-9658-a34d3f96f535
md"""
### The field, live
Finally, the whole picture in motion: the true temperature field, the posterior mean, where the two differ, and how far the posterior standard deviation of every cell rises above its minimum over the displayed time window, with the source temperatures below. The absolute standard deviation is dominated by the field noise ``(\beta A)^{-1}`` and barely changes in time; the increase caused by a blind sensor is small but spreads from that wall into the plate. Thick bars mark the source segments, red while their sensor is offline. As the day goes on, the hot and cold spots rotate around the grid.
That baseline is computed separately for each cell; it is not a separate all-sensors-online inference run. The temperature scale clips at the 98th percentile of absolute non-source temperatures, and the error scale at the 99.5th percentile, so saturated pixels can exceed the colour-bar limit.
"""
# ╔═╡ 66403346-f86c-5a9e-ab33-3d5cc30f88ed
begin
to_grid(v) = reshape(v, m, m)' # row r, column c ↔ cell(r, c)
rc(i) = (div(i - 1, m) + 1, mod(i - 1, m) + 1)
interior = setdiff(1:n_grid, reduce(vcat, source_cells))
temp_max = sort(abs.(vec(X_true[:, interior])))[round(Int, 0.98 * end)] # let the hottest cells saturate
temp_lims = (-temp_max, temp_max)
X_err = abs.(X_mean - X_true)
err_lims = (0.0, sort(vec(X_err))[round(Int, 0.995 * end)])
X_extra_std = X_std .- minimum(X_std; dims = 1) # increase over each cell's minimum across this time window
std_lims = (0.0, maximum(X_extra_std))
function field_panel(values, title; clims, color)
p = heatmap(values; title, clims, color, aspect_ratio = :equal, yflip = true, framestyle = :none, colorbar = true, titlefontsize = 11)
return p
end
function mark_cells!(p, t)
for (k, cells) in enumerate(source_cells)
offline = haskey(dropouts, k) && t in dropouts[k]
plot!(p, last.(rc.(cells)), first.(rc.(cells)), color = offline ? :red : :black, linewidth = offline ? 5 : 3, label = false)
end
return p
end
function field_frame(t)
offline = count(k -> t in dropouts[k], keys(dropouts))
p_true = mark_cells!(field_panel(to_grid(X_true[t, :]), "true temperature"; clims = temp_lims, color = cgrad(:RdBu, rev = true)), t)
p_mean = mark_cells!(field_panel(to_grid(X_mean[t, :]), "posterior mean"; clims = temp_lims, color = cgrad(:RdBu, rev = true)), t)
p_err = mark_cells!(field_panel(to_grid(X_err[t, :]), "|posterior mean − truth|"; clims = err_lims, color = :viridis), t)
p_std = mark_cells!(field_panel(to_grid(X_extra_std[t, :]), "std above temporal minimum"; clims = std_lims, color = :viridis), t)
p_b = plot(legend = :outertop, legend_columns = 4, xlabel = "time", ylabel = "source temperature", xlims = (1, T), ylims = (-12, 12), left_margin = 6Plots.mm)
for (k, window) in dropouts
vspan!(p_b, [first(window), last(window)], color = :gray, alpha = 0.15, label = false)
end
for k in 1:K
plot!(p_b, 1:t, B_mean[1:t, k], ribbon = 2 .* B_std[1:t, k], color = k, fillalpha = 0.2, label = source_names[k])
plot!(p_b, 1:t, B_true[1:t, k], color = k, linestyle = :dot, linewidth = 1, label = false)
end
vline!(p_b, [t], color = :black, linewidth = 1, label = false)
plot(p_true, p_mean, p_err, p_std, p_b; layout = @layout([a b c d; e{0.42h}]), size = (1300, 580), margin = 2Plots.mm, bottom_margin = 5Plots.mm,
plot_title = "t = $t · $offline sensor(s) offline", plot_titlefontsize = 12)
end
nothing
end
# ╔═╡ e1a3304c-28ff-54df-9bc6-adbece63a05b
begin
anim = @animate for t in [1:T; fill(T, 12)]
field_frame(t)
end
animation_path = joinpath(mktempdir(), "heat-grid-smoother.gif")
gif(anim, animation_path, fps = 6, show_msg = false)
LocalResource(animation_path)
end
# ╔═╡ 2a317bbc-4c78-5760-9cbf-2fcaa7422da0
md"""
A single frame from the middle of the first outage:
"""
# ╔═╡ 8b367441-e479-585a-949f-d6e2dfc499a5
begin
field_frame(55)
end
# ╔═╡ dcbad74e-c9a4-51ca-a15d-4d5fe978cdd3
md"""
## Takeaways
- `GaussianCoupling(x[i], -A[i, j])` turns a matrix off-diagonal into a pairwise Gaussian coupling.
- Sparse matrices become sparse factor graphs, and local message updates recover the global solution.
- Loopy graphs need initial messages; `μ(x) = NormalMeanVariance(0, 1e6)` is a broad computational starting point, not an additional model prior.
- Under suitable convergence conditions, the converged means solve ``Ax=b`` exactly.
- Converged variances are exact on trees, but generally approximate on graphs with cycles.
- The same graph solves a whole sequence of systems: give ``b`` a temporal prior, feed it observations (missing ones included) and message passing turns into a smoother over time.
- A prior written as a submodel is reused with one `~` line; its hyperparameters go on the Markov blanket, so the caller decides whether they are learned (with constraints scoped by `for q in submodel`) or plugged in as constants.
- When ``b`` becomes latent, the normaliser of ``\exp(b^\mathsf{T}x - \frac{1}{2}x^\mathsf{T}Ax)`` turns into a negative precision on ``b`` that scales with ``\beta``. Its effect is most noticeable where the data go quiet: keep ``\beta`` small and verify that the posterior precision stays positive definite through the longest dropout. Message passing itself never needs ``A^{-1}``; the reference check uses four grid solves to form ``A^{-1}M``.
For production linear solves, Julia's specialized factorizations remain the natural default. The message-passing formulation becomes interesting when the graph is distributed, when local updates matter, or when the linear system is one component of a larger probabilistic model.
"""
# ╔═╡ 8e840433-3350-5e68-adcb-cad78d2a2e65
md"""
## Environment
This notebook was validated with Julia 1.11.9. Its embedded Project and Manifest specify the complete environment. The versions running now are:
| Package | Version |
|:--|:--|
| Julia | $(VERSION) |
| RxInfer | $(pkgversion(RxInfer)) |
| ReactiveMP | $(pkgversion(ReactiveMP)) |
| Plots | $(pkgversion(Plots)) |
| PlutoUI | $(pkgversion(PlutoUI)) |
| StableRNGs | $(pkgversion(StableRNGs)) |
| OneHotArrays | $(pkgversion(OneHotArrays)) |
"""
# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
OneHotArrays = "0b1bfda6-eb8a-41d2-88d8-f5af5cad476f"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
ReactiveMP = "a194aa59-28ba-4574-a09c-4a745416d6e3"
RxInfer = "86711068-29c9-4ff7-b620-ae75d7495b3d"
StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
[compat]
OneHotArrays = "~0.2.11"
Plots = "~1.41.7"
PlutoUI = "~0.7.83"
ReactiveMP = "~6.5.0"
RxInfer = "~5.5.2"
StableRNGs = "~1.0.4"
julia = "1.11"
"""
# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised
julia_version = "1.11.9"
manifest_format = "2.0"
project_hash = "0af848a8805b767c98c91c9281ea9d1f99d9e169"
[[deps.ADTypes]]
deps = ["PrecompileTools"]
git-tree-sha1 = "629de23e1c16911b439dabd2303c08af9575b226"
uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
version = "1.24.0"
[deps.ADTypes.extensions]
ADTypesChainRulesCoreExt = "ChainRulesCore"
ADTypesConstructionBaseExt = "ConstructionBase"
ADTypesEnzymeCoreExt = "EnzymeCore"
[deps.ADTypes.weakdeps]
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869"
[[deps.AbstractPlutoDingetjes]]
git-tree-sha1 = "6c3913f4e9bdf6ba3c08041a446fb1332716cbc2"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.4.0"
[[deps.Accessors]]
deps = ["CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "MacroTools"]
git-tree-sha1 = "7063ad1083578215c7c4bf410368150abe8d5524"
uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
version = "0.1.45"
[deps.Accessors.extensions]
AxisKeysExt = "AxisKeys"
IntervalSetsExt = "IntervalSets"
LinearAlgebraExt = "LinearAlgebra"
StaticArraysExt = "StaticArrays"
StructArraysExt = "StructArrays"
TestExt = "Test"
UnitfulExt = "Unitful"
[deps.Accessors.weakdeps]
AxisKeys = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5"
IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "daa72978cd7a624246e894a4f4f067706d4e17e2"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "4.7.0"
weakdeps = ["SparseArrays", "StaticArrays"]
[deps.Adapt.extensions]
AdaptSparseArraysExt = "SparseArrays"
AdaptStaticArraysExt = "StaticArrays"
[[deps.AliasTables]]
deps = ["PtrArrays", "Random"]
git-tree-sha1 = "9876e1e164b144ca45e9e3198d0b689cadfed9ff"
uuid = "66dad0bd-aa9a-41b7-9441-69ab47430ed8"
version = "1.1.3"
[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
version = "1.1.2"
[[deps.ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "d57bd3762d308bded22c3b82d033bff85f6195c6"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.4.0"
[[deps.ArrayInterface]]
deps = ["Adapt", "LinearAlgebra"]
git-tree-sha1 = "daf5b2aab5b1c1fdcb65b05883cdb4b18abac1b9"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "7.30.1"
[deps.ArrayInterface.extensions]
ArrayInterfaceAMDGPUExt = "AMDGPU"
ArrayInterfaceBandedMatricesExt = "BandedMatrices"
ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices"
ArrayInterfaceCUDAExt = "CUDA"
ArrayInterfaceCUDSSExt = ["CUDSS", "CUDA"]
ArrayInterfaceChainRulesCoreExt = "ChainRulesCore"
ArrayInterfaceChainRulesExt = "ChainRules"
ArrayInterfaceFillArraysExt = "FillArrays"
ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
ArrayInterfaceGPUArraysCoreTrackerExt = ["GPUArraysCore", "Tracker"]
ArrayInterfaceMetalExt = "Metal"
ArrayInterfaceReverseDiffExt = "ReverseDiff"
ArrayInterfaceSparseArraysExt = "SparseArrays"
ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore"
ArrayInterfaceTrackerExt = "Tracker"
[deps.ArrayInterface.weakdeps]
AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e"
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e"
ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
[[deps.ArrayLayouts]]
deps = ["FillArrays", "LinearAlgebra", "StaticArrays"]
git-tree-sha1 = "e0b47732a192dd59b9d079a06d04235e2f833963"
uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
version = "1.12.2"
weakdeps = ["SparseArrays"]
[deps.ArrayLayouts.extensions]
ArrayLayoutsSparseArraysExt = "SparseArrays"
[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
version = "1.11.0"
[[deps.Atomix]]
deps = ["UnsafeAtomics"]
git-tree-sha1 = "b8651b2eb5796a386b0398a20b519a6a6150f75c"
uuid = "a9b6321e-bd34-4604-b9c9-b65b8de01458"
version = "1.1.3"
[deps.Atomix.extensions]
AtomixCUDAExt = "CUDA"
AtomixMetalExt = "Metal"
AtomixOpenCLExt = "OpenCL"
AtomixoneAPIExt = "oneAPI"
[deps.Atomix.weakdeps]
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
OpenCL = "08131aa3-fb12-5dee-8b74-c09406e224a2"
oneAPI = "8f75cd03-7ff8-4ecb-9b8f-daf728133b1b"
[[deps.BFloat16s]]
deps = ["LinearAlgebra", "Printf", "Random"]
git-tree-sha1 = "e386db8b4753b42caac75ac81d0a4fe161a68a97"
uuid = "ab4f0b2a-ad5b-11e8-123f-65d77653426b"
version = "0.6.1"
[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
version = "1.11.0"
[[deps.BayesBase]]
deps = ["Distributions", "DomainSets", "LinearAlgebra", "Random", "SpecialFunctions", "StaticArrays", "Statistics", "StatsAPI", "StatsBase", "StatsFuns", "TinyHugeNumbers"]
git-tree-sha1 = "e644fb61dcc7b0df1be931230d2b9bc3a96d912f"
uuid = "b4ee3484-f114-42fe-b91c-797d54a0c67e"
version = "1.5.9"
weakdeps = ["FastCholesky"]
[deps.BayesBase.extensions]
FastCholeskyExt = "FastCholesky"
[[deps.BitFlags]]
git-tree-sha1 = "bbe1079eecf9c9fbb52765193ad2bae27ae09bc8"
uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35"
version = "0.1.10"
[[deps.BitSetTuples]]
deps = ["TupleTools"]
git-tree-sha1 = "aa19428fb6ad21db22f8568f068de4f443d3bacc"
uuid = "0f2f92aa-23a3-4d05-b791-88071d064721"
version = "1.1.5"
[[deps.BlockArrays]]
deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra"]
git-tree-sha1 = "75c9c4d41f387b58ac7ecac17a02062f4cf8e92a"
uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
version = "1.10.0"
[deps.BlockArrays.extensions]
BlockArraysAdaptExt = "Adapt"
BlockArraysBandedMatricesExt = "BandedMatrices"