Implemented HUGIN and idempotent architectures. Support for semiring …

…matrices.

docs/Project.toml

@@ -1,6 +1,5 @@
 [deps]
 BayesNets = "ba4760a4-c768-5bed-964b-cf806dc591cb"
-BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
@@ -9,3 +8,4 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
+UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"

docs/literate/kalman.jl

@@ -1,12 +1,10 @@
 # # Kalman Filter
 using AlgebraicInference
-using BenchmarkTools
 using Catlab.Graphics, Catlab.Programs, Catlab.WiringDiagrams
 using Distributions
 using FillArrays
 using LinearAlgebra
 using Random
-using StatsPlots
 # A Kalman filter with ``n`` steps is a probability distribution over states
 # ``(x_1, \dots, x_n)`` and measurements ``(y_1, \dots, y_n)`` determined by the equations
 # ```math
@@ -104,8 +102,8 @@ end
 to_graphviz(make_diagram(5), box_labels=:name; junction_labels=:variable)
 # We generate ``100`` points of data and solve the prediction problem.
 n = 100
-
 diagram = make_diagram(n)
+data = generate_data(n)
 
 hom_map = Dict{String, DenseGaussianSystem{Float64}}(
     "state" => normal(Zeros(2), 100I(2)),
@@ -116,40 +114,14 @@ ob_map = Dict(
     "X" => 2,
     "Y" => 2)
 
-data = generate_data(n)
-
 for i in 1:n
     hom_map["y$i"] = normal(data[:, i], Zeros(2, 2))
 end
 
 problem = InferenceProblem(diagram, hom_map, ob_map)
 
-solver = init(problem)
+Σ = solve(problem)
 
-Σ = solve!(solver)
-
-m = mean(Σ)
-
-round.(m; digits=4)
+round.(mean(Σ); digits=4)
 #
-V = cov(Σ)
-
-round.(V; digits=4)
-# The smoothing problem involves finding the posterior means and covariances of the states
-# ``(x_1, \dots, x_{n - 1})`` given observations of ``(y_1, \dots, y_n)``.
-#
-# Calling `mean(solver)` computes a dictionary with the posterior mean of every variable in
-# the model.
-ms = mean(solver)
-
-x = Matrix{Float64}(undef, 2, n)
-y = Matrix{Float64}(undef, 2, n)
-
-for i in 1:n
-    x[:, i] .= ms[1, i]
-    y[:, i] .= ms[2, i]
-end
-
-plot()
-plot!(x[1, :], label="x₁")
-plot!(x[2, :], label="x₂")
+round.(cov(Σ); digits=4)

docs/literate/pixel.jl

@@ -0,0 +1,86 @@
+# # Pixel Arrays
+using AlgebraicInference
+using Catlab
+using UnicodePlots
+
+using Catlab.CategoricalAlgebra.FinRelations: BoolRig
+# Spivak, Dobson, Kumari, and Wu, *Pixel Arrays*: A fast and elementary method for solving nonlinear systems.
+#
+# A pixel array is an array with entries in the boolean semiring.
+const PixelArray{N} = Array{BoolRig, N}
+
+function PixelArray(f::Function, xdim::NamedTuple, ydim::NamedTuple, tol::Real)
+    rows = xdim.resolution
+    cols = ydim.resolution
+
+    result = PixelArray{2}(undef, rows, cols)
+
+    xstep = (xdim.upper - xdim.lower) / xdim.resolution
+    ystep = (ydim.upper - ydim.lower) / ydim.resolution
+
+    for i in 1:rows, j in 1:cols
+        xval = xdim.lower + (i - 0.5) * xstep
+        yval = ydim.lower + (j - 0.5) * ystep
+
+        try
+            result[i, j] = -tol < f(xval, yval) < tol
+        catch
+            result[i, j] = false
+        end
+    end
+
+    result
+end
+# For plotting.
+function Base.isless(x::Int, y::BoolRig)
+    y == true ? x < 1 : x < 0
+end;
+# Example 2.4.1
+function f₁(x::Real, z::Real)
+    0 - (cos(log(z^2 + x / 10^3)) − x + 1 / z / 10^5)
+end
+
+function f₂(w::Real, y::Real)
+    2 - (cosh(w + y / 10^3) + y + w / 10^4)
+end
+
+function f₃(x::Real, y::Real)
+    1 - (tan(x + y) / (x - 2) / (x + 3) / y^2)
+end
+
+w = x = y = z = (lower = -3, upper = 3, resolution = 125)
+
+tol = .2
+
+A₁ = PixelArray(f₁, x, z, tol)
+A₂ = PixelArray(f₂, w, y, tol)
+A₃ = PixelArray(f₃, x, y, tol);
+# f₁(x, z)
+spy(A₁)
+# f₂(w, y)
+spy(A₂)
+# f₃(x, y)
+spy(A₃)
+#
+diagram = @relation (w, z) where (w::n, x::n, y::n, z::n) begin
+    f₁(x, z)
+    f₂(w, y)
+    f₃(x, y)
+end
+
+to_graphviz(diagram; box_labels=:name, junction_labels=:variable)
+#
+hom_map = Dict{Symbol, PixelArray}(
+    :f₁ => A₁,
+    :f₂ => A₂,
+    :f₃ => A₃)
+
+ob_map = Dict(
+    :n => 125)
+
+problem = InferenceProblem(diagram, hom_map, ob_map)
+
+A = solve(problem)
+
+spy(A)
+#

docs/src/api.md

@@ -36,8 +36,9 @@ init(::InferenceProblem, ::EliminationAlgorithm, ::SupernodeType, ::Architecture
 InferenceSolver
 
 solve!(::InferenceSolver)
-mean(::InferenceSolver)
-rand(::AbstractRNG, ::InferenceSolver)
+mean(::InferenceSolver{AncestralSampler()})
+rand(::InferenceSolver{AncestralSampler()})
+rand(::AbstractRNG, ::InferenceSolver{AncestralSampler()})
 ```
 
 ## Elimination
@@ -61,4 +62,7 @@ MaximalSupernode
 ArchitectureType
 ShenoyShafer
 LauritzenSpiegelhalter
+HUGIN
+Idempotent
+AncestralSampler
 ```

docs/src/generated/kalman.md

@@ -6,13 +6,11 @@ EditURL = "../../literate/kalman.jl"
 
 ````@example kalman
 using AlgebraicInference
-using BenchmarkTools
 using Catlab.Graphics, Catlab.Programs, Catlab.WiringDiagrams
 using Distributions
 using FillArrays
 using LinearAlgebra
 using Random
-using StatsPlots
 ````
 
 A Kalman filter with ``n`` steps is a probability distribution over states
@@ -123,8 +121,8 @@ We generate ``100`` points of data and solve the prediction problem.
 
 ````@example kalman
 n = 100
-
 diagram = make_diagram(n)
+data = generate_data(n)
 
 hom_map = Dict{String, DenseGaussianSystem{Float64}}(
     "state" => normal(Zeros(2), 100I(2)),
@@ -135,48 +133,18 @@ ob_map = Dict(
     "X" => 2,
     "Y" => 2)
 
-data = generate_data(n)
-
 for i in 1:n
     hom_map["y$i"] = normal(data[:, i], Zeros(2, 2))
 end
 
 problem = InferenceProblem(diagram, hom_map, ob_map)
 
-solver = init(problem)
-
-Σ = solve!(solver)
+Σ = solve(problem)
 
-m = mean(Σ)
-
-round.(m; digits=4)
+round.(mean(Σ); digits=4)
 ````
 
 ````@example kalman
-V = cov(Σ)
-
-round.(V; digits=4)
-````
-
-The smoothing problem involves finding the posterior means and covariances of the states
-``(x_1, \dots, x_{n - 1})`` given observations of ``(y_1, \dots, y_n)``.
-
-Calling `mean(solver)` computes a dictionary with the posterior mean of every variable in
-the model.
-
-````@example kalman
-ms = mean(solver)
-
-x = Matrix{Float64}(undef, 2, n)
-y = Matrix{Float64}(undef, 2, n)
-
-for i in 1:n
-    x[:, i] .= ms[1, i]
-    y[:, i] .= ms[2, i]
-end
-
-plot()
-plot!(x[1, :], label="x₁")
-plot!(x[2, :], label="x₂")
+round.(cov(Σ); digits=4)
 ````
 

docs/src/generated/pixel.md

@@ -0,0 +1,123 @@
+```@meta
+EditURL = "../../literate/pixel.jl"
+```
+
+# Pixel Arrays
+
+````@example pixel
+using AlgebraicInference
+using Catlab
+using UnicodePlots
+
+using Catlab.CategoricalAlgebra.FinRelations: BoolRig
+````
+
+Spivak, Dobson, Kumari, and Wu, *Pixel Arrays*: A fast and elementary method for solving nonlinear systems.
+
+A pixel array is an array with entries in the boolean semiring.
+
+````@example pixel
+const PixelArray{N} = Array{BoolRig, N}
+
+function PixelArray(f::Function, xdim::NamedTuple, ydim::NamedTuple, tol::Real)
+    rows = xdim.resolution
+    cols = ydim.resolution
+
+    result = PixelArray{2}(undef, rows, cols)
+
+    xstep = (xdim.upper - xdim.lower) / xdim.resolution
+    ystep = (ydim.upper - ydim.lower) / ydim.resolution
+
+    for i in 1:rows, j in 1:cols
+        xval = xdim.lower + (i - 0.5) * xstep
+        yval = ydim.lower + (j - 0.5) * ystep
+
+        try
+            result[i, j] = -tol < f(xval, yval) < tol
+        catch
+            result[i, j] = false
+        end
+    end
+
+    result
+end
+````
+
+For plotting.
+
+````@example pixel
+function Base.isless(x::Int, y::BoolRig)
+    y == true ? x < 1 : x < 0
+end;
+nothing #hide
+````
+
+Example 2.4.1
+
+````@example pixel
+function f₁(x::Real, z::Real)
+    0 - (cos(log(z^2 + x / 10^3)) − x + 1 / z / 10^5)
+end
+
+function f₂(w::Real, y::Real)
+    2 - (cosh(w + y / 10^3) + y + w / 10^4)
+end
+
+function f₃(x::Real, y::Real)
+    1 - (tan(x + y) / (x - 2) / (x + 3) / y^2)
+end
+
+w = x = y = z = (lower = -3, upper = 3, resolution = 125)
+
+tol = .2
+
+A₁ = PixelArray(f₁, x, z, tol)
+A₂ = PixelArray(f₂, w, y, tol)
+A₃ = PixelArray(f₃, x, y, tol);
+nothing #hide
+````
+
+f₁(x, z)
+
+````@example pixel
+spy(A₁)
+````
+
+f₂(w, y)
+
+````@example pixel
+spy(A₂)
+````
+
+f₃(x, y)
+
+````@example pixel
+spy(A₃)
+````
+
+````@example pixel
+diagram = @relation (w, z) where (w::n, x::n, y::n, z::n) begin
+    f₁(x, z)
+    f₂(w, y)
+    f₃(x, y)
+end
+
+to_graphviz(diagram; box_labels=:name, junction_labels=:variable)
+````
+
+````@example pixel
+hom_map = Dict{Symbol, PixelArray}(
+    :f₁ => A₁,
+    :f₂ => A₂,
+    :f₃ => A₃)
+
+ob_map = Dict(
+    :n => 125)
+
+problem = InferenceProblem(diagram, hom_map, ob_map)
+
+A = solve(problem)
+
+spy(A)
+````
+