Merge pull request #852 from AlgebraicJulia/jpf/sheaves

Add Sheaf implementation

docs/literate/sketches/sheaves.jl

@@ -0,0 +1,111 @@
+# # Basic Sheaf Constructions
+# This example shows how to use a sheaf Catlab.
+# We use the DiagramTopology on FinSet and the Free Vector Space functor
+# to illustrate this process.
+
+using Test
+using Catlab.CategoricalAlgebra
+using Catlab.CategoricalAlgebra.Categories
+using Catlab.Sheaves
+import Catlab.CategoricalAlgebra.Categories: do_ob_map, do_hom_map
+using Catlab.CategoricalAlgebra.Matrices: MatrixDom
+import Catlab.Sheaves: pullback_matrix, FinSetCat
+
+# The Topology that we are using is the `DiagramTopology`, which says that a cover of an object `c` is any diagram with `c` as its colimit.
+# The cocone of the diagram gives you a basis for the sieve that covers `c`.
+# The legs of the cocone can be interpreted as open subsets of `c` whose union is `c`.
+# Because colimits are functorial and the base category is adhesive, the diagram topology is well defined.
+#
+# We have two different implementations of equivalent sheaves.
+# The `VectSheaf` implementation uses julia functions and the `VectSheafMat` implementation uses sparse matrices.
+
+VectSheaf = Sheaf(DiagramTopology(), FVectPullback())
+VectSheafMat = Sheaf(DiagramTopology(), FMatPullback())
+
+# We want to introduce a cover of the set ${1...4}$ with two open sets ${1,2,3}$ and ${1,2,4}$.
+# To do that, we will construct a pushout.
+f = FinFunction([1,2], 3)
+g = FinFunction([1,2], 3)
+S = ColimCover(pushout(f,g))
+
+# The main API of a sheaf is that you can:
+# 1. restrict sections along a morphism,
+# 1. check if a collection of local sections match on the overlaps,
+# 1. and extend a collection of local sections if they match.
+
+# We can restrict a global section along the first leg
+v = [1,2,3,4]
+v₁ = restrict(VectSheaf, v, legs(S)[1])
+
+# Or the second leg.
+v₂ = restrict(VectSheaf, v, legs(S)[2])
+
+# Notice that we got two different values, but that the first two entries of those restricted vectors agree.
+# That is because our two legs have an overlap in the first two dimensions.
+# If we restrict again by the arrows in our diagram (f,g) we should get the same answer.
+
+restrict(VectSheaf, v₁, f)  == restrict(VectSheaf, v₂, g)
+
+# A sheaf requires:
+# - a notion of dividing up a space into pieces (the open coverage),
+# - a way of measuring data over the pieces (the set of sections given by the functor on objects), along with
+# - ways of relating measurements on a piece to measurements on a subpiece (the restriction maps given by the functor).
+# Functoriality of the the sheaf tells you that restricting the same data
+# along commuting paths will always give you the same data on the common piece.
+
+# The sheaf condition requires that you can compute an extension of local data to global data.
+# There is no way to derive this operation from the contravariant functor of the sheaf
+# so providing it is part of the API in defining a new sheaf.
+extend(VectSheaf, S, [[1.0, 2,3], [1,2.0,6]])
+
+# If the local sections don't match, then extending will fail.
+@test_throws MatchingError extend(VectSheaf, S, [[1.0, 2,3], [1,3.0,6]])
+# extend(VectSheaf, S, [[1.0, 2,3], [1,3.0,6]])
+
+# Pushouts are the covers with only two subobjects, but the sheaf works on diagrams of any size.
+
+D = FreeDiagram(FinSet.([3,2,3,3]), # list of objects
+ [ # list of (hom, src, tgt) tuples
+  (FinFunction([1,2], 3), 2,1),
+  (FinFunction([1,2], 3), 2,3),
+  (FinFunction([1,2], 3), 2,4)
+  ]
+)
+
+K = ColimCover(D)
+
+section_data = [Float64[1,2,3],
+   Float64[1,2],
+   Float64[1,2,5],
+   Float64[1,2,6]]
+
+v = extend(VectSheaf, K, section_data; debug=true)
+
+# Our two sheaves should agree, because they are just two different implementations of the same sheaf.
+
+global_section = extend(VectSheafMat, K, section_data)
+v == global_section
+
+# If you put in bad data, you get MatchingErrors
+section_data_bad = [Float64[1,2,3],
+   Float64[1,2],
+   Float64[1,3,5],
+   Float64[1,3,6]]
+
+@test_throws MatchingError extend(VectSheaf, K, section_data_bad)
+@test_throws MatchingError extend(VectSheafMat, K, section_data_bad)
+
+# if we disable the checks, VectSheafMat will solve a least squares problem instead of last write wins.
+extend(VectSheafMat, K, section_data_bad, check=false)
+
+# Last write wins definition is different.
+extend(VectSheaf, K, section_data_bad, check=false) 
+
+# You can diagnose a matching error from the exception that extend throws.
+try
+  extend(VectSheaf, K, section_data_bad)
+catch e
+  e
+end
+
+# This concludes our discussion of the sheaf API.

docs/make.jl

@@ -53,6 +53,7 @@ makedocs(
         "generated/sketches/meets.md",
         "generated/sketches/smc.md",
         "generated/sketches/cat_elements.md",
+        "generated/sketches/sheaves.md",
         ],
       "Graphs" => Any[
         "generated/graphs/graphs.md",
@@ -80,6 +81,7 @@ makedocs(
       "apis/wiring_diagrams.md",
       "apis/graphics.md",
       "apis/programs.md",
+      "apis/sheaves.md",
     ],
     "Developer Docs" => Any[
       "devdocs/style.md",

docs/src/apis/sheaves.md

@@ -0,0 +1,11 @@
+# [Sheaves](@id sheaves)
+
+Sheaves are useful for modeling spatial data. The classic example is continuous functions over a topological space. If you have a function defined at every point in a topological space, then you can think about covers of your space and how that function restricts to subspaces of their domain. When you have locally defined functions that agree on overlapping domains, you can extend them to one function defined on the union of their domains.
+
+Sheaves are a generalization of this idea. The motivating example of sheaves of vectors is implemented and described in the vignettes.
+
+```@autodocs
+Modules = [
+  Sheaves
+  ]
+```

src/Catlab.jl

@@ -9,12 +9,14 @@ include("categorical_algebra/CategoricalAlgebra.jl")
 include("wiring_diagrams/WiringDiagrams.jl")
 include("graphics/Graphics.jl")
 include("programs/Programs.jl")
+include("sheaves/Sheaves.jl")
 
 @reexport using .Theories
 @reexport using .Graphs
 @reexport using .CategoricalAlgebra
 @reexport using .WiringDiagrams
 @reexport using .Graphics
 @reexport using .Programs
+@reexport using .Sheaves
 
 end

src/sheaves/FVect.jl

@@ -0,0 +1,90 @@
+struct FinVect <: Category{FinSet, Function, SmallCatSize} end
+
+"""    FVectPullback{T} where T <: Number
+
+The contravariant free vector space functor. 
+The returned function f✶ restricts via precomposition.
+"""
+struct FVectPullback <: Functor{FinSetOpT, FinVect} end
+dom(::FVectPullback) = FinSetOpT
+codom(::FVectPullback) = FinVect
+do_ob_map(::FVectPullback, n::FinSet) = n
+do_hom_map(::FVectPullback, f::FinFunction) = (v->v[f.(dom(f))])
+# as a callable functor this would be: FVectPullback = Functor(identity, f->(v->v[f.(dom(f))]), OppositeCat(FinSetCat()), FinVect())
+
+"""    FreeVect{T} where T <: Number
+
+The covariant free vector space functor. 
+The returned  function f✶ sums over preimages.
+"""
+FVectPushforward = Functor(identity, # identity on objects
+  # specify the pushforward construction
+  f->(v->map(dom(f)) do i
+    sum(v[j] for j in preimage(f,i))
+  end),
+  # covariant functor from FinSetCat to FinVect
+  FinSetCat(), FinVect()
+)
+
+const FinMat = TypeCat(MatrixDom, Matrix)
+const FinMatT = typeof(FinMat)
+
+# call a matrix on a vector multiplies by it.
+function (M::SparseArrays.SparseMatrixCSC{Int64, Int64})(x::AbstractVector)
+  M*x
+end
+
+function pullback_matrix(f::FinFunction)
+    n = length(dom(f))
+    sparse(1:n, f.(dom(f)), ones(Int,n), dom(f).n, codom(f).n)
+end
+
+function pushforward_matrix(f::FinFunction)
+  pullback_matrix(f)'
+end
+
+struct FMatPullback <: Functor{FinSetOpT, FinMatT} end
+dom(::FMatPullback) = FinSetOpT
+codom(::FMatPullback) = FinMat
+do_ob_map(::FMatPullback, n::FinSet) = MatrixDom{AbstractMatrix}(n.n)
+do_hom_map(::FMatPullback, f::FinFunction) = pullback_matrix(f)
+
+FMatPushforward = Functor(n->MatrixDom(n.n),
+  pushforward_matrix,
+  FinSetCat(), FinVect()
+)
+
+"""    extend(X::Sheaf{T, FVectPullback}, cover::ColimCover, sections::Vector{Vector{R}}; check=true, debug=false) where {T<:DiagramTopology, R}
+
+This method implements the extension operation for the diagram topology on FinSet for the Free Vector Space functor.
+The implementation does copies the value of the ith section into the jth spot as indicated by the legs of the cocone.
+"""
+function extend(X::Sheaf{T, FVectPullback}, cover::ColimCover, sections::Vector{Vector{R}}; check=true, debug=false) where {T<:DiagramTopology, R}
+  length(sections) == length(legs(cover)) || throw(ArgumentError("There are $(length(sections)) but only $(length(legs(cover))) legs in the cover"))
+  v = zeros(R, apex(cover))
+  if check
+    match_errors = diagnose_match(X, cover, sections; debug=debug)
+    length(match_errors) == 0 || throw(MatchingError(match_errors))
+  end
+  for (i, l) in enumerate(legs(cover))
+    for j in dom(l)
+      # copy the value of the ith section into the jth spot
+      # last write wins is fine,
+      # because the sections have to agree on the overlapping indices
+      v[l(j)] = sections[i][j]
+    end
+  end
+  return v
+end
+
+function extend(X::Sheaf{T, FMatPullback}, cover::ColimCover, sections::Vector{Vector{R}};check=true, debug=false) where {T<:DiagramTopology, R}
+  length(sections) == length(legs(cover)) || throw(ArgumentError("There are $(length(sections)) but only $(length(legs(cover))) legs in the cover"))
+  v = zeros(R, apex(cover))
+  if check
+    match_errors = diagnose_match(X, cover, sections; debug=debug)
+    length(match_errors) == 0 || throw(MatchingError(match_errors))
+  end
+  f = copair(legs(cover))
+  M = do_hom_map(functor(X), f)
+  return Float64.(M) \ direct_sum(sections)
+end