Sheaf Morphisms

This example demonstrates constructing sheaf morphisms between Euclidean sheaves built over a cycle graph, including checking the naturality condition and composing morphisms.

The types SheafMorphism and ComplexMorphism, together with is_morphism and compose, are part of the CellularSheaves package — they do not need to be defined here.

using CellularSheaves
using LinearAlgebra

Build a cycle sheaf by adding stalks and edges

Instead of first creating a Graph and calling sheaf_from_graph, we instantiate an empty Euclidean sheaf and populate it using add_vertex_stalk! and add_sheaf_edge!.

n = 6
stalk_dim = 3
3

Start with an empty Euclidean sheaf (no vertex stalks)

F = EuclideanSheaf{Float64}(Int[])
A network sheaf with 0 vertex stalks and 0 edge stalks.

Add n vertex stalks of dimension stalk_dim.

for i in 1:n
    add_vertex_stalk!(F, stalk_dim)
end

Add edges to form a cycle. Restriction maps are stalk_dim x stalk_dim identity matrices, so each edge stalk also has dimension stalk_dim.

for i in 1:n
    v1 = i
    v2 = (i % n) + 1
    rm = Matrix{Float64}(I, stalk_dim, stalk_dim) # Identity map as restriction map
    add_sheaf_edge!(F, v1, v2, rm, rm)
end

println(F)
A network sheaf with 6 vertex stalks and 6 edge stalks.

Inspect the sheaf Laplacian and nearest global section

L = sheaf_laplacian_matrix(F)
println("Sheaf Laplacian is a ", size(L, 1), " x ", size(L, 2), " matrix.")
Sheaf Laplacian is a 18 x 18 matrix.

Start from a random 0-cochain and project to the nearest global section

x0 = rand(sum(vertex_stalks(F)))
gs = try
    nearest_global_section(F, x0)
catch e
    println("Could not compute nearest global section during docs build: ", e)
    zeros(length(x0))
end
println("Norm of coboundary on nearest global section: ", norm(coboundary_map(F) * gs))
Norm of coboundary on nearest global section: 1.5992607561452414e-15

The printed zero (or tiny numerical value) above indicates gs is close to a true global section.

Start with an empty Euclidean sheaf (no vertex stalks)

G = EuclideanSheaf{Float64}(Int[])
A network sheaf with 0 vertex stalks and 0 edge stalks.

Add n vertex stalks of dimension stalk_dim.

for i in 1:n
    add_vertex_stalk!(G, stalk_dim)
end

Add edges to form a cycle. Restriction maps are (stalk_dim - 1) x stalk_dim, so each edge stalk has dimension stalk_dim - 1.

for i in 1:n
    v1 = i
    v2 = (i % n) + 1
    rm = Matrix{Float64}(I, stalk_dim - 1, stalk_dim)
    add_sheaf_edge!(G, v1, v2, rm, rm)
end

println(G)

Vmap = collect(1:n)
Emap = collect(1:n)
Vmaps = [Matrix{Float64}(I, stalk_dim, stalk_dim) for i in 1:n]
Emaps = [Matrix{Float64}(I, stalk_dim - 1, stalk_dim) for i in 1:n]

spec = SheafMorphism(Vmap, Emap, Vmaps, Emaps)
cm = ComplexMorphism(F, G, spec)

is_morphism(cm, F, G)
true

The morphism condition holds, so the coboundary maps commute with the morphism matrices.

dF = coboundary_map(F)
dG = coboundary_map(G)
12×18 BlockSparseMatrixCSC{Float64, Int64}:
 1.0  0.0  0.0  -1.0  -0.0  -0.0   0.0  …   0.0   0.0   0.0   0.0   0.0   0.0
 0.0  1.0  0.0  -0.0  -1.0  -0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
 1.0  0.0  0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0  -1.0  -0.0  -0.0
 0.0  1.0  0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0  -0.0  -1.0  -0.0
 0.0  0.0  0.0   1.0   0.0   0.0  -1.0      0.0   0.0   0.0   0.0   0.0   0.0
 0.0  0.0  0.0   0.0   1.0   0.0  -0.0  …   0.0   0.0   0.0   0.0   0.0   0.0
 0.0  0.0  0.0   0.0   0.0   0.0   1.0      0.0   0.0   0.0   0.0   0.0   0.0
 0.0  0.0  0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
 0.0  0.0  0.0   0.0   0.0   0.0   0.0     -1.0  -0.0  -0.0   0.0   0.0   0.0
 0.0  0.0  0.0   0.0   0.0   0.0   0.0     -0.0  -1.0  -0.0   0.0   0.0   0.0
 0.0  0.0  0.0   0.0   0.0   0.0   0.0  …   1.0   0.0   0.0  -1.0  -0.0  -0.0
 0.0  0.0  0.0   0.0   0.0   0.0   0.0      0.0   1.0   0.0  -0.0  -1.0  -0.0

Build a third sheaf H on the same graph topology but with 1-dimensional edge stalks. Vertex stalks match G (stalk_dim each). Edge restriction maps are 1 x stalk_dim projection matrices that pick the first coordinate.

H = EuclideanSheaf{Float64}(Int[])
for i in 1:n
  add_vertex_stalk!(H, stalk_dim)
end
for i in 1:n
  v1 = i
  v2 = (i % n) + 1
  rmH = zeros(1, stalk_dim)  # restriction: project vertex R^d -> R^1 by taking first coordinate
  rmH[1,1] = 1.0
  add_sheaf_edge!(H, v1, v2, rmH, rmH)
end

println(H)
A network sheaf with 6 vertex stalks and 6 edge stalks.

Build a SheafMorphism spec from G -> H. Vertex maps are identity (3x3), edge maps project G's edge stalk (size = stalk_dim - 1) to H's edge stalk (1)

Vmap_GH = collect(1:n)
Emap_GH = collect(1:n)
Vmaps_GH = [Matrix{Float64}(I, stalk_dim, stalk_dim) for i in 1:n]
edimsG = edge_stalk_dimensions(G)
Emaps_GH = [zeros(1, edimsG[i]) for i in 1:length(edimsG)]
for i in 1:length(edimsG)
  Emaps_GH[i][1,1] = 1.0  # project to first coordinate of G's edge stalk
end

specGH = make_sheaf_morphism_spec(Vmap_GH, Emap_GH, Vmaps_GH, Emaps_GH)
cmGH = ComplexMorphism(G, H, specGH)

is_morphism(cmGH, G, H)

dH = coboundary_map(H)
6×18 BlockSparseMatrixCSC{Float64, Int64}:
 1.0  0.0  0.0  -1.0  -0.0  -0.0   0.0  …   0.0   0.0   0.0   0.0   0.0   0.0
 1.0  0.0  0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0  -1.0  -0.0  -0.0
 0.0  0.0  0.0   1.0   0.0   0.0  -1.0      0.0   0.0   0.0   0.0   0.0   0.0
 0.0  0.0  0.0   0.0   0.0   0.0   1.0      0.0   0.0   0.0   0.0   0.0   0.0
 0.0  0.0  0.0   0.0   0.0   0.0   0.0     -1.0  -0.0  -0.0   0.0   0.0   0.0
 0.0  0.0  0.0   0.0   0.0   0.0   0.0  …   1.0   0.0   0.0  -1.0  -0.0  -0.0

Compose morphisms using the library's compose function

Compose the structured specs F -> G and G -> H to get a spec for F -> H.

composite_morphism = ComplexMorphism(F, H, compose(spec, specGH))
ComplexMorphism([1.0 0.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 1.0], [1.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])

The V component is

composite_morphism.V
6×6-blocked 18×18 BlockMatrix{Float64}:
 1.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0  0.0  …  0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0
 0.0  1.0  0.0  │  0.0  0.0  0.0  │  0.0  0.0     0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0
 0.0  0.0  1.0  │  0.0  0.0  0.0  │  0.0  0.0
 ───────────────┼─────────────────┼──────────     0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0
 ───────────────┼─────────────────┼─────
 0.0  0.0  0.0  │  1.0  0.0  0.0  │  0.0  0.0     0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0
 0.0  0.0  0.0  │  0.0  1.0  0.0  │  0.0  0.0     0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0
 0.0  0.0  0.0  │  0.0  0.0  1.0  │  0.0  0.0
 ───────────────┼─────────────────┼──────────  …  0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0
 ───────────────┼─────────────────┼─────
 0.0  0.0  0.0  │  0.0  0.0  0.0  │  1.0  0.0     0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0
 0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0  1.0     0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0
 0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0  0.0
 ───────────────┼─────────────────┼──────────     0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0
 ───────────────┼─────────────────┼─────
 0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0  0.0     0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0
 0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0  0.0  …  0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0
 0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0  0.0
 ───────────────┼─────────────────┼──────────     1.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0
 ───────────────┼─────────────────┼─────
 0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0  0.0     0.0  1.0  0.0  │  0.0  0.0  0.0  │  0.0
 0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0  0.0     0.0  0.0  1.0  │  0.0  0.0  0.0  │  0.0
 0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0  0.0
 ───────────────┼─────────────────┼──────────     0.0  0.0  0.0  │  1.0  0.0  0.0  │  0.0
 ───────────────┼─────────────────┼─────
 0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0  0.0  …  0.0  0.0  0.0  │  0.0  1.0  0.0  │  0.0
 0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0  0.0     0.0  0.0  0.0  │  0.0  0.0  1.0  │  0.0
 0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0  0.0     0.0  0.0  0.0  │  0.0  0.0  0.0  │  1.0

The E component is

composite_morphism.E
6×6-blocked 6×18 BlockMatrix{Float64}:
 1.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0  0.0
 ───────────────┼─────────────────┼──────────  …  0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0
 ───────────────┼─────────────────┼─────
 0.0  0.0  0.0  │  1.0  0.0  0.0  │  0.0  0.0
 ───────────────┼─────────────────┼──────────     0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0
 ───────────────┼─────────────────┼─────
 0.0  0.0  0.0  │  0.0  0.0  0.0  │  1.0  0.0
 ───────────────┼─────────────────┼──────────     0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0
 ───────────────┼─────────────────┼─────
 0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0  0.0
 ───────────────┼─────────────────┼──────────     0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0
 ───────────────┼─────────────────┼─────
 0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0  0.0
 ───────────────┼─────────────────┼──────────     0.0  1.0  0.0  │  0.0  0.0  0.0  │  0.0
 ───────────────┼─────────────────┼─────
 0.0  0.0  0.0  │  0.0  0.0  0.0  │  0.0  0.0  …  0.0  0.0  0.0  │  0.0  1.0  0.0  │  0.0

The composite is supposed to be the same as the result of directly multiplying the ComplexMorphism matrices:

ϕ = compose(cm, cmGH)
norm(ϕ.V - composite_morphism.V)
0.0

and the E component is

norm(ϕ.E - composite_morphism.E)
0.0