From 2703d703f787becc6f8853ccdaea246d8ddd8443 Mon Sep 17 00:00:00 2001 From: Luke Morris <70283489+lukem12345@users.noreply.github.com> Date: Tue, 24 Sep 2024 15:00:44 -0400 Subject: [PATCH 01/16] Add GeometryBasics to docs/Project.toml From df29cef08aca3caf82cc48a10960363c5660d8cb Mon Sep 17 00:00:00 2001 From: Luke Morris <70283489+lukem12345@users.noreply.github.com> Date: Tue, 24 Sep 2024 19:46:05 -0400 Subject: [PATCH 02/16] Add ACSets to docs/Project.toml From 2af28349208738ed02204b8a2d95ff154f3c5087 Mon Sep 17 00:00:00 2001 From: Owen Lynch Date: Thu, 30 May 2024 00:34:08 -0400 Subject: [PATCH 03/16] started work on simplicial complexes --- Project.toml | 1 + src/GeometricMaps.jl | 27 ++++++++++++++++++++++ src/SimplicialComplexes.jl | 47 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 75 insertions(+) create mode 100644 src/GeometricMaps.jl create mode 100644 src/SimplicialComplexes.jl diff --git a/Project.toml b/Project.toml index 7bf99b06..b4b1b417 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Catlab = "134e5e36-593f-5add-ad60-77f754baafbe" DataMigrations = "0c4ad18d-0c49-4bc2-90d5-5bca8f00d6ae" +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" diff --git a/src/GeometricMaps.jl b/src/GeometricMaps.jl new file mode 100644 index 00000000..2b928dfd --- /dev/null +++ b/src/GeometricMaps.jl @@ -0,0 +1,27 @@ +""" +A particularly simple definition for a geometric map from an n-truncated +simplicial set S to an n-truncated simplicial set T is just a function + +S_0 -> T_n x Δ^n + +This sends each vertex in S to a pair of an n-simplex in T and a point inside +the generic geometric n-simplex. + +However, this naive definition can run into problems. Specifically, given a +k-simplex in S, its vertices must all be sent into the same n-simplex. But +this would imply that connected components in S must all be sent into the +same n-simplex. In order to fix this, we should change the definition to + +S_0 -> \sum_{k=0}^n T_k x int(\Delta^k) + +Then the condition for a k-simplex in S +""" + +struct Point + simplex::Int + coordinates::Vector{Float64} +end + +struct GeometricMap + values::Vector{Point} +end diff --git a/src/SimplicialComplexes.jl b/src/SimplicialComplexes.jl new file mode 100644 index 00000000..5ef1e228 --- /dev/null +++ b/src/SimplicialComplexes.jl @@ -0,0 +1,47 @@ +module SimplicialComplexes + +function add_0_cells(d::HasDeltaSet, t::Trie{Int, Int}) + for v in vertices(d_set) + t[[v]] = v + end +end + +function add_1_cells(d::HasDeltaSet, t::Trie{Int, Int}) + for e in edges(d) + v1, v2 = src(d, e), tgt(d, e) + v1 < v2 || error("Degenerate or unsorted edge: $e") + haskey(t, [v1, v2]) && error("Duplicate edge: $e") + t[[v1, v2]] = e + end +end + +function add_2_cells(d::HasDeltaSet, t::Trie{Int, Int}) + for tr in triangles(d) + v1, v2, v3 = triangle_vertices(d, t) + v1 < v2 < v3 || error("Degenerate or unsorted trangle: $tr") + haskey(t, [v1, v2, v3]) && error("Duplicate triangle: $tr") + t[[v1, v2, v3]] = tr + end +end + +struct SimplicialComplex{D} + delta_set::D + complexes::Trie{Int, Int} + + function SimplicialComplex(d::D) where {D<:AbstractDeltaSet1D} + t = Trie{Int, Int}() + add_0_cells(d, t) + add_1_cells(d, t) + new{D}(d, t) + end + + function SimplicialComplex(d::D) where {D<:AbstractDeltaSet2D} + t = Trie{Int, Int}() + add_0_cells(d, t) + add_1_cells(d, t) + add_2_cells(d, t) + new{D}(d, t) + end +end + +end From e7c2f2eba941671eb1e96b0df062aeeb0f6d9b0d Mon Sep 17 00:00:00 2001 From: Owen Lynch Date: Wed, 5 Jun 2024 11:45:39 -0700 Subject: [PATCH 04/16] Added basic simplicial complexes, vendored trie implementation Co-authored-by: Kevin Carlson --- Project.toml | 1 - src/CombinatorialSpaces.jl | 4 + src/GeometricMaps.jl | 38 +++++++- src/SimplicialComplexes.jl | 63 +++++++++++-- src/SimplicialSets.jl | 17 +++- src/Tries.jl | 179 ++++++++++++++++++++++++++++++++++++ test/SimplicialComplexes.jl | 20 ++++ test/Tries.jl | 87 ++++++++++++++++++ test/runtests.jl | 8 ++ 9 files changed, 402 insertions(+), 15 deletions(-) create mode 100644 src/Tries.jl create mode 100644 test/SimplicialComplexes.jl create mode 100644 test/Tries.jl diff --git a/Project.toml b/Project.toml index b4b1b417..7bf99b06 100644 --- a/Project.toml +++ b/Project.toml @@ -10,7 +10,6 @@ Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Catlab = "134e5e36-593f-5add-ad60-77f754baafbe" DataMigrations = "0c4ad18d-0c49-4bc2-90d5-5bca8f00d6ae" -DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" diff --git a/src/CombinatorialSpaces.jl b/src/CombinatorialSpaces.jl index 96a59770..7908df2c 100644 --- a/src/CombinatorialSpaces.jl +++ b/src/CombinatorialSpaces.jl @@ -2,16 +2,20 @@ module CombinatorialSpaces using Reexport +include("Tries.jl") include("ArrayUtils.jl") include("CombinatorialMaps.jl") include("ExteriorCalculus.jl") include("SimplicialSets.jl") include("DiscreteExteriorCalculus.jl") include("MeshInterop.jl") +include("SimplicialComplexes.jl") include("FastDEC.jl") include("Meshes.jl") +@reexport using .Tries @reexport using .SimplicialSets +@reexport using .SimplicialComplexes @reexport using .DiscreteExteriorCalculus @reexport using .FastDEC @reexport using .MeshInterop diff --git a/src/GeometricMaps.jl b/src/GeometricMaps.jl index 2b928dfd..db65e494 100644 --- a/src/GeometricMaps.jl +++ b/src/GeometricMaps.jl @@ -22,6 +22,42 @@ struct Point coordinates::Vector{Float64} end -struct GeometricMap +function vertices(d::DeltaSet, dim::Int, i::Int) + if dim == 0 + [i] + elseif dim == 1 + [src(d, i), tgt(d, i)] + elseif dim == 2 + triangle_vertices(d, i) + else + error("we can't do above 2-dimensions") + end +end + +function vertices(d::DeltaSet, p::Point) + vertices(d, length(p.coordinates) - 1, p.simplex) +end + +function check_simplex_exists( + codom::SimplicialComplex{D'} + points::Vector{Point} +) where {D, D'} + verts = vcat([vertices(codom.delta_set, p) for p in points]) + sort!(verts) + unique!(verts) + haskey(codom.complexes, verts) || + error("edge $e cannot map into a valid complex") +end + +struct GeometricMap{D, D'} + dom::SimplicialComplex{D} + codom::SimplicialComplex{D'} values::Vector{Point} + function GeometricMap( + dom::SimplicialComplex{D}, + codom::SimplicialComplex{D'}, + values::Vector{Point} + ) where {D <: AbstractDeltaSet1D, D'} + + end end diff --git a/src/SimplicialComplexes.jl b/src/SimplicialComplexes.jl index 5ef1e228..80a89d27 100644 --- a/src/SimplicialComplexes.jl +++ b/src/SimplicialComplexes.jl @@ -1,32 +1,35 @@ module SimplicialComplexes +export SimplicialComplex, VertexList, has_simplex +using ..Tries +using ..SimplicialSets function add_0_cells(d::HasDeltaSet, t::Trie{Int, Int}) - for v in vertices(d_set) + for v in vertices(d) t[[v]] = v end end function add_1_cells(d::HasDeltaSet, t::Trie{Int, Int}) for e in edges(d) - v1, v2 = src(d, e), tgt(d, e) - v1 < v2 || error("Degenerate or unsorted edge: $e") - haskey(t, [v1, v2]) && error("Duplicate edge: $e") - t[[v1, v2]] = e + vs = sort([src(d, e), tgt(d, e)]) + allunique(vs) || error("Degenerate edge: $e") + haskey(t, vs) && error("Duplicate edge: $e") + t[vs] = e end end function add_2_cells(d::HasDeltaSet, t::Trie{Int, Int}) for tr in triangles(d) - v1, v2, v3 = triangle_vertices(d, t) - v1 < v2 < v3 || error("Degenerate or unsorted trangle: $tr") - haskey(t, [v1, v2, v3]) && error("Duplicate triangle: $tr") - t[[v1, v2, v3]] = tr + vs = sort(triangle_vertices(d, tr)) + allunique(vs) || error("Degenerate triangle: $tr") + haskey(t, vs) && error("Duplicate triangle: $tr") + t[vs] = tr end end struct SimplicialComplex{D} delta_set::D - complexes::Trie{Int, Int} + cache::Trie{Int, Int} function SimplicialComplex(d::D) where {D<:AbstractDeltaSet1D} t = Trie{Int, Int}() @@ -44,4 +47,44 @@ struct SimplicialComplex{D} end end +struct VertexList #XX parameterize by n? remember to replace sort! w sort + vs::Vector{Int} # must be sorted + function VertexList(vs::Vector{Int}; sorted=false) + new(sorted ? vs : sort(vs)) + end + function VertexList(d::HasDeltaSet, s::Simplex{n,0}) where n + new(sort(simplex_vertices(d,s))) + end +end + +Base.length(s::VertexList) = length(s.vs) +has_simplex(sc::SimplicialComplex,s::VertexList) = haskey(sc.cache, s.vs) + +Base.getindex(v::VertexList, i::Int) = v.vs[i] + +function Base.getindex(sc::SimplicialComplex, s::VertexList)::Simplex + has_simplex(sc,s) || error("Simplex not found: $s") + Simplex{length(s)}(sc.cache[s.vs]) end + +function Base.union(vs1::VertexList, vs2::VertexList) + out = Int[] + i, j = 1 + while (i <= length(vs1)) && (j <= length(vs2)) + v1, v2 = vs1[i], vs2[j] + if (v1 == v2) + push!(out, v1) + i += 1 + j += 1 + elseif (v1 <= v2) + push!(out, v1) + i += 1 + else + push!(out, v2) + j += 1 + end + end + VertexList(out, sorted=true) +end + +end \ No newline at end of file diff --git a/src/SimplicialSets.jl b/src/SimplicialSets.jl index da5b9b07..c707e035 100644 --- a/src/SimplicialSets.jl +++ b/src/SimplicialSets.jl @@ -35,7 +35,7 @@ export Simplex, V, E, Tri, Tet, SimplexChain, VChain, EChain, TriChain, TetChain tetrahedron_triangles, tetrahedron_edges, tetrahedron_vertices, ntetrahedra, tetrahedra, add_tetrahedron!, glue_tetrahedron!, glue_sorted_tetrahedron!, glue_sorted_tet_cube!, is_manifold_like, nonboundaries, - star, St, closed_star, St̄, link, Lk + star, St, closed_star, St̄, link, Lk, simplex_vertices using LinearAlgebra: det using SparseArrays @@ -602,8 +602,19 @@ const Tri = Simplex{2} """ const Tet = Simplex{3} -""" Wrapper for chain of oriented simplices of dimension `n`. -""" +# could generalize to Simplex{n, N} +function simplex_vertices(s::HasDeltaSet, x::Simplex{n,0}) where n + simplex_vertices(Val{n}, s, x) +end + +function simplex_vertices(::Type{Val{n}},s::HasDeltaSet,x::Simplex{n,0}) where n + n == 0 && return [x.data] + n == 1 && return edge_vertices(s, x.data) + n == 2 && return triangle_vertices(s, x.data) + n == 3 && return tetrahedron_vertices(s, x.data) +end + +""" Wrapper for simplex chain of dimension `n`.""" @vector_struct SimplexChain{n} const VChain = SimplexChain{0} diff --git a/src/Tries.jl b/src/Tries.jl new file mode 100644 index 00000000..d6278653 --- /dev/null +++ b/src/Tries.jl @@ -0,0 +1,179 @@ +module Tries +export Trie, keys_with_prefix, partial_path, +find_prefixes, subtrie + +# vendored in from https://github.com/JuliaCollections/DataStructures.jl/blob/master/src/trie.jl + +mutable struct Trie{K,V} + value::Union{Some{V}, Nothing} + children::Dict{K,Trie{K,V}} + function Trie{K,V}() where {K,V} + self = new{K,V}() + self.value = nothing + self.children = Dict{K,Trie{K,V}}() + return self + end +end + +function Trie{K,V}(ks, vs) where {K,V} + return Trie{K,V}(zip(ks, vs)) +end + +function Trie{K,V}(kv) where {K,V} + t = Trie{K,V}() + for (k,v) in kv + t[k] = v + end + return t +end + +Trie() = Trie{Any,Any}() +Trie(ks::AbstractVector{K}, vs::AbstractVector{V}) where {K,V} = Trie{eltype(K),V}(ks, vs) +Trie(kv::AbstractVector{Tuple{K,V}}) where {K,V} = Trie{eltype(K),V}(kv) +Trie(kv::AbstractDict{K,V}) where {K,V} = Trie{eltype(K),V}(kv) +Trie(ks::AbstractVector{K}) where {K} = Trie{eltype(K),Nothing}(ks, similar(ks, Nothing)) + +hasvalue(t::Trie) = !isnothing(t.value) + + +function Base.setindex!(t::Trie{K,V}, val, key) where {K,V} + value = convert(V, val) # we don't want to iterate before finding out it fails + node = t + for char in key + if !haskey(node.children, char) + node.children[char] = Trie{K,V}() + end + node = node.children[char] + end + node.value = Some(value) +end + +function Base.getindex(t::Trie, key) + node = subtrie(t, key) + if !isnothing(node) && hasvalue(node) + return something(node.value) + end + throw(KeyError("key not found: $key")) +end + +function subtrie(t::Trie, prefix) + node = t + for char in prefix + if !haskey(node.children, char) + return nothing + else + node = node.children[char] + end + end + return node +end + +function Base.haskey(t::Trie, key) + node = subtrie(t, key) + !isnothing(node) && hasvalue(node) +end + +function Base.get(t::Trie, key, notfound) + node = subtrie(t, key) + if !isnothing(node) && hasvalue(node) + return something(node.value) + end + return notfound +end + +_concat(prefix::String, char::Char) = string(prefix, char) +_concat(prefix::Vector{T}, char::T) where {T} = vcat(prefix, char) + +_empty_prefix(::Trie{Char,V}) where {V} = "" +_empty_prefix(::Trie{K,V}) where {K,V} = K[] + +function Base.keys(t::Trie{K,V}, + prefix=_empty_prefix(t), + found=Vector{typeof(prefix)}()) where {K,V} + if hasvalue(t) + push!(found, prefix) + end + for (char,child) in t.children + keys(child, _concat(prefix, char), found) + end + return found +end + +function keys_with_prefix(t::Trie, prefix) + st = subtrie(t, prefix) + st != nothing ? keys(st,prefix) : [] +end + +# The state of a TrieIterator is a pair (t::Trie, i::Int), +# where t is the Trie which was the output of the previous iteration +# and i is the index of the current character of the string. +# The indexing is potentially confusing; +# see the comments and implementation below for details. +struct TrieIterator + t::Trie + str +end + +# At the start, there is no previous iteration, +# so the first element of the state is undefined. +# We use a "dummy value" of it.t to keep the type of the state stable. +# The second element is 0 +# since the root of the trie corresponds to a length 0 prefix of str. +function Base.iterate(it::TrieIterator, (t, i) = (it.t, 0)) + if i == 0 + return it.t, (it.t, firstindex(it.str)) + elseif i > lastindex(it.str) || !(it.str[i] in keys(t.children)) + return nothing + else + t = t.children[it.str[i]] + return (t, (t, nextind(it.str, i))) + end +end + +partial_path(t::Trie, str) = TrieIterator(t, str) +Base.IteratorSize(::Type{TrieIterator}) = Base.SizeUnknown() + +""" + find_prefixes(t::Trie, str) + +Find all keys from the `Trie` that are prefix of the given string + +# Examples +```julia-repl +julia> t = Trie(["A", "ABC", "ABCD", "BCE"]) + +julia> find_prefixes(t, "ABCDE") +3-element Vector{AbstractString}: + "A" + "ABC" + "ABCD" + +julia> t′ = Trie([1:1, 1:3, 1:4, 2:4]); + +julia> find_prefixes(t′, 1:5) +3-element Vector{UnitRange{Int64}}: + 1:1 + 1:3 + 1:4 + +julia> find_prefixes(t′, [1,2,3,4,5]) +3-element Vector{Vector{Int64}}: + [1] + [1, 2, 3] + [1, 2, 3, 4] +``` +""" +function find_prefixes(t::Trie, str::T) where {T} + prefixes = T[] + it = partial_path(t, str) + idx = 0 + for t in it + if hasvalue(t) + push!(prefixes, str[firstindex(str):idx]) + end + idx = nextind(str, idx) + end + return prefixes +end + +end \ No newline at end of file diff --git a/test/SimplicialComplexes.jl b/test/SimplicialComplexes.jl new file mode 100644 index 00000000..eee6ad13 --- /dev/null +++ b/test/SimplicialComplexes.jl @@ -0,0 +1,20 @@ +module TestSimplicialComplexes +using Test +using CombinatorialSpaces.SimplicialSets +using CombinatorialSpaces.SimplicialComplexes + +# Triangulated commutative square. +s = DeltaSet2D() +add_vertices!(s, 4) +glue_triangle!(s, 1, 2, 3) +glue_triangle!(s, 1, 4, 3) + +sc = SimplicialComplex(s) +vl = VertexList(s,Simplex{2}(1)) + +@test vl.vs == [1,2,3] +@test has_simplex(sc,vl) +@test !has_simplex(sc,VertexList([1,2,4])) +@test sc[vl] == Simplex{2}(1) + +end \ No newline at end of file diff --git a/test/Tries.jl b/test/Tries.jl new file mode 100644 index 00000000..1d284c4f --- /dev/null +++ b/test/Tries.jl @@ -0,0 +1,87 @@ +module TestTries +# Originally taken from: https://github.com/JuliaCollections/DataStructures.jl/blob/master/test/test_trie.jl + +using Test +using CombinatorialSpaces.Tries + +@testset "Core Functionality" begin + t = Trie{Char,Int}() + t["amy"] = 56 + t["ann"] = 15 + t["emma"] = 30 + t["rob"] = 27 + t["roger"] = 52 + t["kevin"] = Int8(11) + + @test haskey(t, "roger") + @test get(t, "rob", nothing) == 27 + @test sort(keys(t)) == ["amy", "ann", "emma", "kevin", "rob", "roger"] + @test t["rob"] == 27 + @test sort(keys_with_prefix(t, "ro")) == ["rob", "roger"] +end + +@testset "Constructors" begin + ks = ["amy", "ann", "emma", "rob", "roger"] + vs = [56, 15, 30, 27, 52] + kvs = collect(zip(ks, vs)) + @test isa(Trie(ks, vs), Trie{Char,Int}) + @test isa(Trie(kvs), Trie{Char,Int}) + @test isa(Trie(Dict(kvs)), Trie{Char,Int}) + @test isa(Trie(ks), Trie{Char,Nothing}) +end + +@testset "partial_path iterator" begin + t = Trie{Char,Int}() + t["rob"] = 27 + t["roger"] = 52 + t["kevin"] = Int8(11) + t0 = t + t1 = t0.children['r'] + t2 = t1.children['o'] + t3 = t2.children['b'] + @test collect(partial_path(t, "b")) == [t0] + @test collect(partial_path(t, "rob")) == [t0, t1, t2, t3] + @test collect(partial_path(t, "robb")) == [t0, t1, t2, t3] + @test collect(partial_path(t, "ro")) == [t0, t1, t2] + @test collect(partial_path(t, "roa")) == [t0, t1, t2] +end + +@testset "partial_path iterator non-ascii" begin + t = Trie(["東京都"]) + t0 = t + t1 = t0.children['東'] + t2 = t1.children['京'] + t3 = t2.children['都'] + @test collect(partial_path(t, "西")) == [t0] + @test collect(partial_path(t, "東京都")) == [t0, t1, t2, t3] + @test collect(partial_path(t, "東京都渋谷区")) == [t0, t1, t2, t3] + @test collect(partial_path(t, "東京")) == [t0, t1, t2] + @test collect(partial_path(t, "東京スカイツリー")) == [t0, t1, t2] +end + +@testset "find_prefixes" begin + t = Trie(["A", "ABC", "ABD", "BCD"]) + prefixes = find_prefixes(t, "ABCDE") + @test prefixes == ["A", "ABC"] +end + +@testset "find_prefixes non-ascii" begin + t = Trie(["東京都", "東京都渋谷区", "東京都新宿区"]) + prefixes = find_prefixes(t, "東京都渋谷区東") + @test prefixes == ["東京都", "東京都渋谷区"] +end + +@testset "non-string indexing" begin + t = Trie{Int,Int}() + t[[1, 2, 3, 4]] = 1 + t[[1, 2]] = 2 + @test haskey(t, [1, 2]) + @test get(t, [1, 2], nothing) == 2 + st = subtrie(t, [1, 2, 3]) + @test keys(st) == [[4]] + @test st[[4]] == 1 + @test find_prefixes(t, [1, 2, 3, 5]) == [[1, 2]] + @test find_prefixes(t, 1:3) == [1:2] +end + +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 5d2f9798..84c8d0bd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,10 @@ module RunTests using Test +@testset "Tries" begin + include("Tries.jl") +end + @testset "CombinatorialMaps" begin include("CombinatorialMaps.jl") end @@ -9,6 +13,10 @@ end include("SimplicialSets.jl") end +@testset "SimplicialComplexes" begin + include("SimplicialComplexes.jl") +end + @testset "ExteriorCalculus" begin include("ExteriorCalculus.jl") include("DiscreteExteriorCalculus.jl") From abb0982d3f8f761fb8b5520fb866b763d9a7bfb8 Mon Sep 17 00:00:00 2001 From: Owen Lynch Date: Wed, 5 Jun 2024 11:51:33 -0700 Subject: [PATCH 05/16] fixed union and added tests --- src/SimplicialComplexes.jl | 12 ++++++++++-- test/SimplicialComplexes.jl | 4 ++++ 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/src/SimplicialComplexes.jl b/src/SimplicialComplexes.jl index 80a89d27..9764f42b 100644 --- a/src/SimplicialComplexes.jl +++ b/src/SimplicialComplexes.jl @@ -58,9 +58,10 @@ struct VertexList #XX parameterize by n? remember to replace sort! w sort end Base.length(s::VertexList) = length(s.vs) +Base.lastindex(s::VertexList) = lastindex(s.vs) has_simplex(sc::SimplicialComplex,s::VertexList) = haskey(sc.cache, s.vs) -Base.getindex(v::VertexList, i::Int) = v.vs[i] +Base.getindex(v::VertexList, i) = v.vs[i] function Base.getindex(sc::SimplicialComplex, s::VertexList)::Simplex has_simplex(sc,s) || error("Simplex not found: $s") @@ -69,7 +70,7 @@ end function Base.union(vs1::VertexList, vs2::VertexList) out = Int[] - i, j = 1 + i, j = 1, 1 while (i <= length(vs1)) && (j <= length(vs2)) v1, v2 = vs1[i], vs2[j] if (v1 == v2) @@ -84,7 +85,14 @@ function Base.union(vs1::VertexList, vs2::VertexList) j += 1 end end + if (i <= length(vs1)) + append!(out, vs1[i:end]) + end + if (j <= length(vs2)) + append!(out, vs2[j:end]) + end VertexList(out, sorted=true) end +#TODO: get a point by barycentric coordinates, maps end \ No newline at end of file diff --git a/test/SimplicialComplexes.jl b/test/SimplicialComplexes.jl index eee6ad13..74526e51 100644 --- a/test/SimplicialComplexes.jl +++ b/test/SimplicialComplexes.jl @@ -17,4 +17,8 @@ vl = VertexList(s,Simplex{2}(1)) @test !has_simplex(sc,VertexList([1,2,4])) @test sc[vl] == Simplex{2}(1) +vl′ = VertexList([1,2])∪VertexList([2,3]) +@test has_simplex(sc,vl′) +@test !has_simplex(sc,VertexList([1,2])∪VertexList([2,4])) + end \ No newline at end of file From 330c5c86d92f1d3ee4cc683e3b89cce6920be0a7 Mon Sep 17 00:00:00 2001 From: AlgebraicJulia Bot <129184742+algebraicjuliabot@users.noreply.github.com> Date: Tue, 14 May 2024 18:00:55 -0400 Subject: [PATCH 06/16] Set version to 0.6.4 From 91502426b126f92108ed6339ed8befac7b96aeae Mon Sep 17 00:00:00 2001 From: Kevin Carlson Date: Fri, 18 Oct 2024 15:33:46 -0400 Subject: [PATCH 07/16] You got you a category of simplicial complexes here --- Project.toml | 1 + src/SimplicialComplexes.jl | 81 +++++++++++++++++++++++++++++++++++-- src/SimplicialSets.jl | 19 +++++---- src/Tries.jl | 4 +- test/SimplicialComplexes.jl | 23 +++++++++++ 5 files changed, 116 insertions(+), 12 deletions(-) diff --git a/Project.toml b/Project.toml index 7bf99b06..58715d76 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.6.7" [deps] ACSets = "227ef7b5-1206-438b-ac65-934d6da304b8" +AlgebraicInterfaces = "23cfdc9f-0504-424a-be1f-4892b28e2f0c" Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Catlab = "134e5e36-593f-5add-ad60-77f754baafbe" diff --git a/src/SimplicialComplexes.jl b/src/SimplicialComplexes.jl index 9764f42b..f47ff1f4 100644 --- a/src/SimplicialComplexes.jl +++ b/src/SimplicialComplexes.jl @@ -1,7 +1,14 @@ +""" +The category of simplicial complexes and Kleisli maps for the convex space monad. +""" module SimplicialComplexes -export SimplicialComplex, VertexList, has_simplex +export SimplicialComplex, VertexList, has_simplex, GeometricPoint, has_point, has_span, GeometricMap, nv, +as_matrix, compose, id using ..Tries using ..SimplicialSets +import AlgebraicInterfaces: dom,codom,compose,id +import LinearAlgebra: I +#import ..SimplicialSets: nv,ne function add_0_cells(d::HasDeltaSet, t::Trie{Int, Int}) for v in vertices(d) @@ -21,7 +28,7 @@ end function add_2_cells(d::HasDeltaSet, t::Trie{Int, Int}) for tr in triangles(d) vs = sort(triangle_vertices(d, tr)) - allunique(vs) || error("Degenerate triangle: $tr") + allunique(vs) || error("Degenerate triangle: $tr") haskey(t, vs) && error("Duplicate triangle: $tr") t[vs] = tr end @@ -31,6 +38,12 @@ struct SimplicialComplex{D} delta_set::D cache::Trie{Int, Int} + function SimplicialComplex(d::DeltaSet0D) + t = Trie{Int, Int}() + add_0_cells(d, t) + new{DeltaSet0D}(d, t) + end + function SimplicialComplex(d::D) where {D<:AbstractDeltaSet1D} t = Trie{Int, Int}() add_0_cells(d, t) @@ -47,7 +60,14 @@ struct SimplicialComplex{D} end end -struct VertexList #XX parameterize by n? remember to replace sort! w sort +#nv(sc::SimplicialComplex) = nv(sc.delta_set) + +for f in [:nv,:ne] + @eval SimplicialSets.$f(sc::SimplicialComplex) = $f(sc.delta_set) +end + + +struct VertexList #XX parameterize by n? vs::Vector{Int} # must be sorted function VertexList(vs::Vector{Int}; sorted=false) new(sorted ? vs : sort(vs)) @@ -94,5 +114,58 @@ function Base.union(vs1::VertexList, vs2::VertexList) VertexList(out, sorted=true) end -#TODO: get a point by barycentric coordinates, maps +#A point in an unspecified simplicial complex, given by its barycentric coordinates. +#Constructed via a dense vector of coordinates. +#XX: This type is maybe more trouble than it's worth? +struct GeometricPoint + bcs::Vector{Float64} #XX: Need a sparse form? + function GeometricPoint(bcs, checked=true) + if checked + sum(bcs) ≈ 1 || error("Barycentric coordinates must sum to 1") + all(x -> 1 ≥ x ≥ 0, bcs) || error("Barycentric coordinates must be between 0 and 1") + end + new(bcs) + end +end +Base.show(p::GeometricPoint) = "GeometricPoint($(p.bcs))" +coords(p::GeometricPoint)=p.bcs + +""" +A simplicial complex contains a geometric point if and only if it contains the combinatorial simplex spanned +by the vertices wrt which the point has a nonzero coordinate. +""" +has_point(sc::SimplicialComplex, p::GeometricPoint) = has_simplex(sc, VertexList(findall(x->x>0,coords(p)))) +""" +A simplicial complex contains the geometric simplex spanned by a list of geometric points if and only if it +contains the combinatorial simplex spanned by all the vertices wrt which some geometric point has a nonzero coordinate. +""" +has_span(sc::SimplicialComplex,ps::Vector{GeometricPoint}) = has_simplex(sc,reduce(union,VertexList.(map(cs->findall(x->x>0,cs),(coords.(ps)))))) + +#geoemtric map between simplicial complexes, given as a list of geometric points in the codomain +#indexed by the 0-simplices of the domain. +struct GeometricMap{D,D′} + dom::SimplicialComplex{D} + cod::SimplicialComplex{D′} + points::Vector{GeometricPoint} + function GeometricMap(sc::SimplicialComplex{D}, sc′::SimplicialComplex{D′}, points::Vector{GeometricPoint},checked=true) where {D,D′} + length(points) == nv(sc) || error("Number of points must match number of vertices in domain") + all(map(x->has_span(sc′,points[x]),keys(sc.cache))) || error("Span of points in simplices of domain must lie in codomain") #lol wrong + new{D,D′}(sc, sc′, points) + end +end +GeometricMap(sc,sc′,points::AbstractArray) = GeometricMap(sc,sc′,GeometricPoint.(eachcol(points))) +dom(f::GeometricMap) = f.dom +codom(f::GeometricMap) = f.cod +#want f(n) to give values[n]? +""" +Returns the data-centric view of f as a matrix whose i-th column +is the coordinates of the image of the i-th vertex under f. +""" +as_matrix(f::GeometricMap) = hcat(coords.(f.points)...) +compose(f::GeometricMap, g::GeometricMap) = GeometricMap(f.dom, g.cod, as_matrix(g)*as_matrix(f)) +id(sc::SimplicialComplex) = GeometricMap(sc,sc,GeometricPoint.(eachcol(I(nv(sc))))) + + + +#TODO: composition of maps! end \ No newline at end of file diff --git a/src/SimplicialSets.jl b/src/SimplicialSets.jl index c707e035..cb3e9830 100644 --- a/src/SimplicialSets.jl +++ b/src/SimplicialSets.jl @@ -15,7 +15,7 @@ exterior derivative) operators. For additional operators, see the module SimplicialSets export Simplex, V, E, Tri, Tet, SimplexChain, VChain, EChain, TriChain, TetChain, SimplexForm, VForm, EForm, TriForm, TetForm, HasDeltaSet, - HasDeltaSet1D, AbstractDeltaSet1D, DeltaSet1D, SchDeltaSet1D, + HasDeltaSet1D, DeltaSet0D, AbstractDeltaSet1D, DeltaSet1D, SchDeltaSet1D, OrientedDeltaSet1D, SchOrientedDeltaSet1D, EmbeddedDeltaSet1D, SchEmbeddedDeltaSet1D, HasDeltaSet2D, AbstractDeltaSet2D, DeltaSet2D, SchDeltaSet2D, @@ -47,12 +47,6 @@ import Catlab.Graphs: src, tgt, nv, ne, vertices, edges, has_vertex, has_edge, add_vertex!, add_vertices!, add_edge!, add_edges! using ..ArrayUtils -# 0-D simplicial sets -##################### - -@present SchDeltaSet0D(FreeSchema) begin - V::Ob -end """ Abstract type for C-sets that contain a delta set of some dimension. @@ -69,6 +63,17 @@ has_vertex(s::HasDeltaSet, v) = has_part(s, :V, v) add_vertex!(s::HasDeltaSet; kw...) = add_part!(s, :V; kw...) add_vertices!(s::HasDeltaSet, n::Int; kw...) = add_parts!(s, :V, n; kw...) +# 0-D simplicial sets +##################### + +@present SchDeltaSet0D(FreeSchema) begin + V::Ob +end +""" A 0-dimensional delta set, aka a set of vertices. +""" +@acset_type DeltaSet0D(SchDeltaSet0D) <: HasDeltaSet +ne(::DeltaSet0D) = error("0-dimensional Δ-set lacks edges.") + # 1D simplicial sets #################### diff --git a/src/Tries.jl b/src/Tries.jl index d6278653..2654c0ef 100644 --- a/src/Tries.jl +++ b/src/Tries.jl @@ -73,7 +73,7 @@ function Base.haskey(t::Trie, key) !isnothing(node) && hasvalue(node) end -function Base.get(t::Trie, key, notfound) +function Base.get(t::Trie, key, notfound) node = subtrie(t, key) if !isnothing(node) && hasvalue(node) return something(node.value) @@ -176,4 +176,6 @@ function find_prefixes(t::Trie, str::T) where {T} return prefixes end +Base.show(io::IO,t::Trie) = print(io,"Trie($(map(k->(k,t[k]),sort(keys(t),by=x->(length(x),sort(x))))))") + end \ No newline at end of file diff --git a/test/SimplicialComplexes.jl b/test/SimplicialComplexes.jl index 74526e51..d9e9754f 100644 --- a/test/SimplicialComplexes.jl +++ b/test/SimplicialComplexes.jl @@ -2,6 +2,7 @@ module TestSimplicialComplexes using Test using CombinatorialSpaces.SimplicialSets using CombinatorialSpaces.SimplicialComplexes +using Catlab:@acset # Triangulated commutative square. s = DeltaSet2D() @@ -10,6 +11,7 @@ glue_triangle!(s, 1, 2, 3) glue_triangle!(s, 1, 4, 3) sc = SimplicialComplex(s) +@test nv(sc) == 4 && ne(sc) == 5 vl = VertexList(s,Simplex{2}(1)) @test vl.vs == [1,2,3] @@ -21,4 +23,25 @@ vl′ = VertexList([1,2])∪VertexList([2,3]) @test has_simplex(sc,vl′) @test !has_simplex(sc,VertexList([1,2])∪VertexList([2,4])) +p = GeometricPoint([0.2,0,0.5,0.3]) +q = GeometricPoint([0,0.2,0.5,0.3]) +r = GeometricPoint([0.5,0,0,0.5]) +s = GeometricPoint([0.5,0,0.5,0]) +t = GeometricPoint([0,0.5,0,0.5]) +@test has_point(sc,p) && !has_point(sc,q) +@test has_span(sc,[r,s]) && !has_span(sc,[r,t]) + +Δ⁰ = SimplicialComplex(@acset DeltaSet0D begin V=1 end) +Δ¹ = SimplicialComplex(@acset DeltaSet1D begin V=2; E=1; ∂v0 = [2]; ∂v1 = [1] end) +f = GeometricMap(Δ⁰,Δ¹,[1/3,2/3]) +A = [0.2 0.4 + 0 0 + 0.5 0 + 0.3 0.6] +g = GeometricMap(Δ¹,sc,A) +@test A == as_matrix(g) +h = compose(f,g) +@test as_matrix(h) ≈ [1/3, 0, 1/6, 1/2] +isc = id(sc) +@test as_matrix(h) == as_matrix(compose(h,isc)) end \ No newline at end of file From c6f203589a42d8b886174ae9a8cbd8aae4cc273d Mon Sep 17 00:00:00 2001 From: Kevin Arlin Date: Tue, 25 Jun 2024 13:53:50 -0700 Subject: [PATCH 08/16] 2D subdivision, mainly --- src/DiscreteExteriorCalculus.jl | 5 +- src/SimplicialComplexes.jl | 129 +++++++++++++++++++++++++------- src/SimplicialSets.jl | 77 ++++++++++++++++++- src/Tries.jl | 21 +++++- test/SimplicialComplexes.jl | 22 +++--- test/SimplicialSets.jl | 3 +- 6 files changed, 212 insertions(+), 45 deletions(-) diff --git a/src/DiscreteExteriorCalculus.jl b/src/DiscreteExteriorCalculus.jl index 87a38435..d2bb9354 100644 --- a/src/DiscreteExteriorCalculus.jl +++ b/src/DiscreteExteriorCalculus.jl @@ -49,7 +49,7 @@ using DataMigrations: @migrate using ..ArrayUtils, ..SimplicialSets using ..SimplicialSets: CayleyMengerDet, operator_nz, ∂_nz, d_nz, cayley_menger, negate -import ..SimplicialSets: ∂, d, volume +import ..SimplicialSets: ∂, d, volume, subdivide abstract type DiscreteFlat end struct DPPFlat <: DiscreteFlat end @@ -199,7 +199,7 @@ end """ Subdivide a 1D delta set. """ -function subdivide(s::HasDeltaSet1D) +function subdivide(s::HasDeltaSet1D) #Should perhaps be in simplicialsets @migrate typeof(s) s begin V => @cases begin v::V @@ -220,6 +220,7 @@ function subdivide(s::HasDeltaSet1D) end end + # 1D oriented dual complex #------------------------- diff --git a/src/SimplicialComplexes.jl b/src/SimplicialComplexes.jl index f47ff1f4..96132624 100644 --- a/src/SimplicialComplexes.jl +++ b/src/SimplicialComplexes.jl @@ -7,6 +7,7 @@ as_matrix, compose, id using ..Tries using ..SimplicialSets import AlgebraicInterfaces: dom,codom,compose,id +import StaticArrays: MVector import LinearAlgebra: I #import ..SimplicialSets: nv,ne @@ -35,38 +36,75 @@ function add_2_cells(d::HasDeltaSet, t::Trie{Int, Int}) end struct SimplicialComplex{D} - delta_set::D - cache::Trie{Int, Int} + delta_set::D + cache::Trie{Int,Int} - function SimplicialComplex(d::DeltaSet0D) - t = Trie{Int, Int}() - add_0_cells(d, t) - new{DeltaSet0D}(d, t) - end + function SimplicialComplex(d::DeltaSet0D) + t = Trie{Int,Int}() + add_0_cells(d, t) + new{DeltaSet0D}(d, t) + end - function SimplicialComplex(d::D) where {D<:AbstractDeltaSet1D} - t = Trie{Int, Int}() - add_0_cells(d, t) - add_1_cells(d, t) - new{D}(d, t) - end + function SimplicialComplex(d::D) where {D<:AbstractDeltaSet1D} + t = Trie{Int,Int}() + add_0_cells(d, t) + add_1_cells(d, t) + new{D}(d, t) + end + + function SimplicialComplex(d::D) where {D<:AbstractDeltaSet2D} + t = Trie{Int,Int}() + add_0_cells(d, t) + add_1_cells(d, t) + add_2_cells(d, t) + new{D}(d, t) + end - function SimplicialComplex(d::D) where {D<:AbstractDeltaSet2D} - t = Trie{Int, Int}() - add_0_cells(d, t) - add_1_cells(d, t) - add_2_cells(d, t) - new{D}(d, t) + #XX Make this work for oriented types, maybe error for embedded types + """ + Build a simplicial complex without a pre-existing delta-set. + + In this case any initial values in the trie are meaningless and will be overwritten. + If you apply this to the cache of a simplicial complex, you may get a non-isomorphic + Δ-set, but it will be isomorphic as a simplicial complex (i.e. after symmetrizing.) + `simplices[1]`` is sorted just for predictability of the output--this guarantees that + the result will have the same indexing for the vertices in its `cache` as in its + `delta_set`. + """ + function SimplicialComplex(dt::Type{D}, t::Trie{Int,Int}) where {D<:HasDeltaSet} + n = dimension(D) + simplices = MVector{n + 1,Vector{Vector{Int}}}(fill([], n + 1)) + for k in keys(t) + push!(simplices[length(k)], k) end + d = D() + for v in sort(simplices[1]) + t[v] = add_vertex!(d) + end + n > 0 && for e in simplices[2] + t[e] = add_edge!(d, t[e[1]], t[e[2]]) + end + n > 1 && for tri in simplices[3] + t[tri] = glue_triangle!(d, t[tri[1]], t[tri[2]], t[tri[3]]) + end + n > 2 && for tet in simplices[4] + t[tet] = glue_tetrahedron!(d, t[tet[1]], t[tet[2]], t[tet[3]], t[tet[4]]) + end + new{D}(d, t) + end end -#nv(sc::SimplicialComplex) = nv(sc.delta_set) +#XX: Should this output something oriented? +""" +Build a simplicial complex from a trie, constructing a delta-set of the minimal +dimension consistent with the trie. +""" +SimplicialComplex(t::Trie{Int,Int}) = SimplicialComplex(DeltaSet(max(height(t)-1,0)),t) -for f in [:nv,:ne] +for f in [:nv,:ne,:ntriangles,:dimension] @eval SimplicialSets.$f(sc::SimplicialComplex) = $f(sc.delta_set) end - struct VertexList #XX parameterize by n? vs::Vector{Int} # must be sorted function VertexList(vs::Vector{Int}; sorted=false) @@ -77,6 +115,43 @@ struct VertexList #XX parameterize by n? end end +""" +Iterator over proper subsimplices of a simplex in reversed binary order. + +Example +```julia-repl +julia> vl = VertexList([1,4,9]) +VertexList([1, 4, 9]) +julia> iter = SubsimplexIterator(vl) +SubsimplexIterator(VertexList([1, 4, 9]), 7) +julia> collect(iter) +7-element Vector{VertexList}: + VertexList(Int64[]) + VertexList([1]) + VertexList([4]) + VertexList([1, 4]) + VertexList([9]) + VertexList([1, 9]) + VertexList([4, 9]) +``` +""" +struct SubsimplexIterator + vl::VertexList + length::Int + #Note that an n-simplex has 2^(n+1)-1 subsimplices, with n+1 vertices. + SubsimplexIterator(vl::VertexList) = new(vl, 2^length(vl.vs)-1) +end +Base.length(iter::SubsimplexIterator) = iter.length +Base.eltype(iter::SubsimplexIterator) = VertexList +function Base.iterate(iter::SubsimplexIterator,i=0) + if i >= iter.length + return nothing + end + ds = digits(i,base=2) + mask = Bool[ds;fill(0,length(iter.vl.vs)-length(ds))] + (VertexList(iter.vl.vs[mask]),i+1) +end + Base.length(s::VertexList) = length(s.vs) Base.lastindex(s::VertexList) = lastindex(s.vs) has_simplex(sc::SimplicialComplex,s::VertexList) = haskey(sc.cache, s.vs) @@ -114,9 +189,12 @@ function Base.union(vs1::VertexList, vs2::VertexList) VertexList(out, sorted=true) end -#A point in an unspecified simplicial complex, given by its barycentric coordinates. -#Constructed via a dense vector of coordinates. + #XX: This type is maybe more trouble than it's worth? +""" +A point in an unspecified simplicial complex, given by its barycentric coordinates. +Constructed via a dense vector of coordinates. +""" struct GeometricPoint bcs::Vector{Float64} #XX: Need a sparse form? function GeometricPoint(bcs, checked=true) @@ -149,7 +227,7 @@ struct GeometricMap{D,D′} points::Vector{GeometricPoint} function GeometricMap(sc::SimplicialComplex{D}, sc′::SimplicialComplex{D′}, points::Vector{GeometricPoint},checked=true) where {D,D′} length(points) == nv(sc) || error("Number of points must match number of vertices in domain") - all(map(x->has_span(sc′,points[x]),keys(sc.cache))) || error("Span of points in simplices of domain must lie in codomain") #lol wrong + all(map(x->has_span(sc′,points[x]),keys(sc.cache))) || error("Span of points in simplices of domain must lie in codomain") new{D,D′}(sc, sc′, points) end end @@ -167,5 +245,4 @@ id(sc::SimplicialComplex) = GeometricMap(sc,sc,GeometricPoint.(eachcol(I(nv(sc)) -#TODO: composition of maps! end \ No newline at end of file diff --git a/src/SimplicialSets.jl b/src/SimplicialSets.jl index cb3e9830..be253f6c 100644 --- a/src/SimplicialSets.jl +++ b/src/SimplicialSets.jl @@ -15,7 +15,7 @@ exterior derivative) operators. For additional operators, see the module SimplicialSets export Simplex, V, E, Tri, Tet, SimplexChain, VChain, EChain, TriChain, TetChain, SimplexForm, VForm, EForm, TriForm, TetForm, HasDeltaSet, - HasDeltaSet1D, DeltaSet0D, AbstractDeltaSet1D, DeltaSet1D, SchDeltaSet1D, + HasDeltaSet1D, DeltaSet, DeltaSet0D, AbstractDeltaSet1D, DeltaSet1D, SchDeltaSet1D, OrientedDeltaSet1D, SchOrientedDeltaSet1D, EmbeddedDeltaSet1D, SchEmbeddedDeltaSet1D, HasDeltaSet2D, AbstractDeltaSet2D, DeltaSet2D, SchDeltaSet2D, @@ -31,11 +31,11 @@ export Simplex, V, E, Tri, Tet, SimplexChain, VChain, EChain, TriChain, TetChain add_vertex!, add_vertices!, add_edge!, add_edges!, add_sorted_edge!, add_sorted_edges!, triangle_edges, triangle_vertices, ntriangles, triangles, - add_triangle!, glue_triangle!, glue_sorted_triangle!, + add_triangle!, glue_triangle!, glue_triangles!, glue_sorted_triangle!, tetrahedron_triangles, tetrahedron_edges, tetrahedron_vertices, ntetrahedra, tetrahedra, add_tetrahedron!, glue_tetrahedron!, glue_sorted_tetrahedron!, glue_sorted_tet_cube!, is_manifold_like, nonboundaries, - star, St, closed_star, St̄, link, Lk, simplex_vertices + star, St, closed_star, St̄, link, Lk, simplex_vertices, dimension, subdivide using LinearAlgebra: det using SparseArrays @@ -63,6 +63,18 @@ has_vertex(s::HasDeltaSet, v) = has_part(s, :V, v) add_vertex!(s::HasDeltaSet; kw...) = add_part!(s, :V; kw...) add_vertices!(s::HasDeltaSet, n::Int; kw...) = add_parts!(s, :V, n; kw...) +""" +Calculate the dimension of a delta set from its acset schema. +Assumes that vertices, edges, triangles, and tetrahedra are +named :V, :E, :Tri, and :Tet respectively. +""" +function dimension(d::HasDeltaSet) + S = acset_schema(d) + :E in ob(S) ? :Tri in ob(S) ? :Tet in ob(S) ? 3 : 2 : 1 : 0 +end +dimension(dt::Type{D}) where {D<:HasDeltaSet} = dimension(D()) + + # 0-D simplicial sets ##################### @@ -330,6 +342,11 @@ function get_edge!(s::HasDeltaSet1D, src::Int, tgt::Int) es = edges(s, src, tgt) isempty(es) ? add_edge!(s, src, tgt) : first(es) end +function glue_triangles!(s,v₀s,v₁s,v₂s; kw...) + for (v₀,v₁,v₂) in zip(v₀s,v₁s,v₂s) + glue_triangle!(s, v₀, v₁, v₂; kw...) + end +end """ Glue a triangle onto a simplicial set, respecting the order of the vertices. """ @@ -338,6 +355,58 @@ function glue_sorted_triangle!(s::HasDeltaSet2D, v₀::Int, v₁::Int, v₂::Int glue_triangle!(s, v₀, v₁, v₂; kw...) end +""" +Designed for simplicial complexes. +""" +function subdivide(s::HasDeltaSet2D) + d = typeof(s)() + #vertices for all simplices + add_vertices!(d, nv(s)) + add_vertices!(d, ne(s)) + add_vertices!(d, ntriangles(s)) + + e_to_v(e) = e + nv(s) + ts_as_vs = (1:ntriangles(s)) .+ (nv(s) + ne(s)) + #edges from source of edge to edge and target of edge to edge + add_edges!(d, subpart(s, :∂v1), e_to_v.(1:ne(s))) + add_edges!(d, subpart(s, :∂v0), e_to_v.(1:ne(s))) + #edges from vertex of triangle to triangle + add_edges!(d, subpart(s, [:∂e2, :∂v1]), ts_as_vs) + add_edges!(d, subpart(s, [:∂e2, :∂v0]), ts_as_vs) + add_edges!(d, subpart(s, [:∂e1, :∂v0]), ts_as_vs) + #edges from edge of triangle to triangle + add_edges!(d, e_to_v.(subpart(s, :∂e2)), ts_as_vs) + add_edges!(d, e_to_v.(subpart(s, :∂e1)), ts_as_vs) + add_edges!(d, e_to_v.(subpart(s, :∂e0)), ts_as_vs) + #triangles from vertex of edge of triangle to triangle + glue_triangles!(d, + subpart(s, [:∂e2, :∂v1]), + e_to_v.(subpart(s, :∂e2)), + ts_as_vs) + glue_triangles!(d, + subpart(s, [:∂e2, :∂v0]), + e_to_v.(subpart(s, :∂e2)), + ts_as_vs) + glue_triangles!(d, + subpart(s, [:∂e1, :∂v1]), + e_to_v.(subpart(s, :∂e1)), + ts_as_vs) + glue_triangles!(d, + subpart(s, [:∂e1, :∂v0]), + e_to_v.(subpart(s, :∂e1)), + ts_as_vs) + glue_triangles!(d, + subpart(s, [:∂e0, :∂v1]), + e_to_v.(subpart(s, :∂e0)), + ts_as_vs) + glue_triangles!(d, + subpart(s, [:∂e0, :∂v0]), + e_to_v.(subpart(s, :∂e0)), + ts_as_vs) + ##XX: orientations? + d +end + # 2D oriented simplicial sets #---------------------------- @@ -585,6 +654,8 @@ volume(::Type{Val{3}}, s::HasDeltaSet3D, t::Int, ::CayleyMengerDet) = # General operators ################### +DeltaSet(n) = eval(Symbol("DeltaSet$(n)D")) + """ Wrapper for simplex or simplices of dimension `n`. See also: [`V`](@ref), [`E`](@ref), [`Tri`](@ref). diff --git a/src/Tries.jl b/src/Tries.jl index 2654c0ef..ec37cb9e 100644 --- a/src/Tries.jl +++ b/src/Tries.jl @@ -1,6 +1,6 @@ module Tries export Trie, keys_with_prefix, partial_path, -find_prefixes, subtrie +find_prefixes, subtrie, height,subdivide_trie # vendored in from https://github.com/JuliaCollections/DataStructures.jl/blob/master/src/trie.jl @@ -35,7 +35,6 @@ Trie(ks::AbstractVector{K}) where {K} = Trie{eltype(K),Nothing}(ks, similar(ks, hasvalue(t::Trie) = !isnothing(t.value) - function Base.setindex!(t::Trie{K,V}, val, key) where {K,V} value = convert(V, val) # we don't want to iterate before finding out it fails node = t @@ -82,7 +81,8 @@ function Base.get(t::Trie, key, notfound) end _concat(prefix::String, char::Char) = string(prefix, char) -_concat(prefix::Vector{T}, char::T) where {T} = vcat(prefix, char) +#was wrong if T is itself a vector type! +_concat(prefix::Vector{T}, char::T) where {T} = vcat(prefix, [char]) _empty_prefix(::Trie{Char,V}) where {V} = "" _empty_prefix(::Trie{K,V}) where {K,V} = K[] @@ -101,8 +101,21 @@ end function keys_with_prefix(t::Trie, prefix) st = subtrie(t, prefix) - st != nothing ? keys(st,prefix) : [] + !isnothing(st) ? keys(st,prefix) : [] +end + +height(t::Trie) = isempty(t.children) ? 0 : 1 + maximum(height.(values(t.children))) + +#= +function subdivide_trie(t::Trie{K,V},help=K[]) where {K,V} + sdt = Trie{Vector{K},V}() + sdt.value = t.value + for k in [k for k in keys(t) if length(k) > 0] + sdt.children[vcat(help,k)] = subdivide(subtrie(t,k),vcat(help,k)) + end + sdt end +=# # The state of a TrieIterator is a pair (t::Trie, i::Int), # where t is the Trie which was the output of the previous iteration diff --git a/test/SimplicialComplexes.jl b/test/SimplicialComplexes.jl index d9e9754f..39344113 100644 --- a/test/SimplicialComplexes.jl +++ b/test/SimplicialComplexes.jl @@ -1,19 +1,23 @@ module TestSimplicialComplexes using Test -using CombinatorialSpaces.SimplicialSets -using CombinatorialSpaces.SimplicialComplexes +using CombinatorialSpaces using Catlab:@acset # Triangulated commutative square. -s = DeltaSet2D() -add_vertices!(s, 4) -glue_triangle!(s, 1, 2, 3) -glue_triangle!(s, 1, 4, 3) +ss = DeltaSet2D() +add_vertices!(ss, 4) +glue_triangle!(ss, 1, 2, 3) +glue_triangle!(ss, 1, 4, 3) -sc = SimplicialComplex(s) -@test nv(sc) == 4 && ne(sc) == 5 -vl = VertexList(s,Simplex{2}(1)) +sc = SimplicialComplex(ss) +@test nv(sc) == 4 && ne(sc) == 5 && ntriangles(sc) == 2 +sc′ = SimplicialComplex(DeltaSet2D,sc.cache).delta_set +@test nv(sc′) == 4 && ne(sc′) == 5 && ntriangles(sc′) == 2 #identifies this up to iso +#awkward height=0 edge case, technically can think of the empty sset as -1-dimensional. +sc′′=SimplicialComplex(Trie{Int,Int}()) +@test dimension(sc′′) == 0 && nv(sc′′) == 0 +vl = VertexList(ss,Simplex{2}(1)) @test vl.vs == [1,2,3] @test has_simplex(sc,vl) @test !has_simplex(sc,VertexList([1,2,4])) diff --git a/test/SimplicialSets.jl b/test/SimplicialSets.jl index 797bd825..280dac56 100644 --- a/test/SimplicialSets.jl +++ b/test/SimplicialSets.jl @@ -121,7 +121,8 @@ glue_triangle!(s, 1, 4, 3) @test triangles(s) == 1:2 @test ne(s) == 5 @test sort(map(Pair, src(s), tgt(s))) == [1=>2, 1=>3, 1=>4, 2=>3, 4=>3] - +sd = subdivide(s) +@test ntriangles(sd) == 12 && ne(sd) == 22 && nv(sd) == 11 # 2D oriented simplicial sets #---------------------------- From e67ded6930afc3f96b8b81b2faf5b44b7b52d4b4 Mon Sep 17 00:00:00 2001 From: Kevin Arlin Date: Tue, 25 Jun 2024 14:52:05 -0700 Subject: [PATCH 09/16] GeometricMap constructor for barycentric subdivision --- src/SimplicialComplexes.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/SimplicialComplexes.jl b/src/SimplicialComplexes.jl index 96132624..51eac506 100644 --- a/src/SimplicialComplexes.jl +++ b/src/SimplicialComplexes.jl @@ -9,6 +9,7 @@ using ..SimplicialSets import AlgebraicInterfaces: dom,codom,compose,id import StaticArrays: MVector import LinearAlgebra: I +import ..DiscreteExteriorCalculus: Barycenter #import ..SimplicialSets: nv,ne function add_0_cells(d::HasDeltaSet, t::Trie{Int, Int}) @@ -243,6 +244,15 @@ as_matrix(f::GeometricMap) = hcat(coords.(f.points)...) compose(f::GeometricMap, g::GeometricMap) = GeometricMap(f.dom, g.cod, as_matrix(g)*as_matrix(f)) id(sc::SimplicialComplex) = GeometricMap(sc,sc,GeometricPoint.(eachcol(I(nv(sc))))) +function GeometricMap(sc::SimplicialComplex,::Barycenter) + dom = SimplicialComplex(subdivide(sc.delta_set)) + #Vertices of dom correspond to vertices, edges, triangles of sc. + mat = zeros(Float64,nv(sc),nv(dom)) + for i in 1:nv(sc) mat[i,i] = 1 end + for i in 1:ne(sc) for n in edge_vertices(sc.delta_set,i) mat[n,nv(sc)+i] = 1/2 end end + for i in 1:ntriangles(sc) for n in triangle_vertices(sc.delta_set,i) mat[n,nv(sc)+ne(sc)+i] = 1/3 end end + GeometricMap(dom,sc,mat) +end end \ No newline at end of file From 6fe8968637daf006997f7ccd7ee7c73975464ef6 Mon Sep 17 00:00:00 2001 From: Kevin Arlin Date: Wed, 26 Jun 2024 13:09:07 -0700 Subject: [PATCH 10/16] primal vector field migrations with dumb subdivision --- src/ArrayUtils.jl | 3 +- src/DiscreteExteriorCalculus.jl | 34 +++++++------------- src/SimplicialComplexes.jl | 8 +++++ src/SimplicialSets.jl | 57 ++++++++++++++++++++++++++++----- 4 files changed, 70 insertions(+), 32 deletions(-) diff --git a/src/ArrayUtils.jl b/src/ArrayUtils.jl index f8f45e82..c6887240 100644 --- a/src/ArrayUtils.jl +++ b/src/ArrayUtils.jl @@ -130,7 +130,8 @@ lazy(::typeof(vcat), args...) = ApplyArray(vcat, args...) # Wrapped arrays ################ -""" Generate struct for a named vector struct. +""" Generate struct for a named vector struct. A special case of `@parts_array_struct` +when we need only vectors of the type. """ macro vector_struct(struct_sig) name, params = parse_struct_signature(struct_sig) diff --git a/src/DiscreteExteriorCalculus.jl b/src/DiscreteExteriorCalculus.jl index d2bb9354..b6a0e0a0 100644 --- a/src/DiscreteExteriorCalculus.jl +++ b/src/DiscreteExteriorCalculus.jl @@ -197,28 +197,7 @@ function dual_triangle_vertices(s::HasDeltaSet2D, t...) s[s[t..., :D_∂e0], :D_∂v0]) end -""" Subdivide a 1D delta set. -""" -function subdivide(s::HasDeltaSet1D) #Should perhaps be in simplicialsets - @migrate typeof(s) s begin - V => @cases begin - v::V - e::E - end - E => @cases begin - e₁::E - e₂::E - end - ∂v1 => begin - e₁ => e - e₂ => e - end - ∂v0 => begin - e₁ => (v∘∂v1) - e₂ => (v∘∂v0) - end - end -end + # 1D oriented dual complex @@ -282,7 +261,7 @@ end make_dual_simplices_1d!(s::AbstractDeltaDualComplex1D) = make_dual_simplices_1d!(s, E) -""" Make dual vertice and edges for dual complex of dimension ≧ 1. +""" Make dual vertices and edges for dual complex of dimension ≧ 1. Although zero-dimensional duality is geometrically trivial (subdividing a vertex gives back the same vertex), we treat the dual vertices as disjoint from the @@ -484,6 +463,15 @@ end tri_center::Hom(Tri, DualV) end + +#type_by_name(header::String,n::Int) = eval(Symbol(header * string(n)*"D")) + +#= +function extract_subdivided_primal_migration(n) + S = +end +=# + """ Abstract type for dual complex of a 2D delta set. """ @abstract_acset_type AbstractDeltaDualComplex2D <: HasDeltaSet2D diff --git a/src/SimplicialComplexes.jl b/src/SimplicialComplexes.jl index 51eac506..618f72ec 100644 --- a/src/SimplicialComplexes.jl +++ b/src/SimplicialComplexes.jl @@ -7,9 +7,11 @@ as_matrix, compose, id using ..Tries using ..SimplicialSets import AlgebraicInterfaces: dom,codom,compose,id +import Base:* import StaticArrays: MVector import LinearAlgebra: I import ..DiscreteExteriorCalculus: Barycenter +import ..DiscreteExteriorCalculus: PrimalVectorField #import ..SimplicialSets: nv,ne function add_0_cells(d::HasDeltaSet, t::Trie{Int, Int}) @@ -254,5 +256,11 @@ function GeometricMap(sc::SimplicialComplex,::Barycenter) GeometricMap(dom,sc,mat) end +function pullback_primal(f::GeometricMap, v::PrimalVectorField{T}) where T + nv(f.cod) == length(v) || error("Vector field must have same number of vertices as codomain") + PrimalVectorField(T.(eachcol(hcat(v.data...)*as_matrix(f)))) +end +*(f::GeometricMap,v::PrimalVectorField) = pullback_primal(f,v) + end \ No newline at end of file diff --git a/src/SimplicialSets.jl b/src/SimplicialSets.jl index be253f6c..53b30bd4 100644 --- a/src/SimplicialSets.jl +++ b/src/SimplicialSets.jl @@ -12,7 +12,7 @@ all geometric applications, namely the boundary and coboundary (discrete exterior derivative) operators. For additional operators, see the `DiscreteExteriorCalculus` module. """ -module SimplicialSets +module SimplicialSets export Simplex, V, E, Tri, Tet, SimplexChain, VChain, EChain, TriChain, TetChain, SimplexForm, VForm, EForm, TriForm, TetForm, HasDeltaSet, HasDeltaSet1D, DeltaSet, DeltaSet0D, AbstractDeltaSet1D, DeltaSet1D, SchDeltaSet1D, @@ -35,7 +35,8 @@ export Simplex, V, E, Tri, Tet, SimplexChain, VChain, EChain, TriChain, TetChain tetrahedron_triangles, tetrahedron_edges, tetrahedron_vertices, ntetrahedra, tetrahedra, add_tetrahedron!, glue_tetrahedron!, glue_sorted_tetrahedron!, glue_sorted_tet_cube!, is_manifold_like, nonboundaries, - star, St, closed_star, St̄, link, Lk, simplex_vertices, dimension, subdivide + star, St, closed_star, St̄, link, Lk, simplex_vertices, dimension, subdivide, + DeltaSet, OrientedDeltaSet, EmbeddedDeltaSet using LinearAlgebra: det using SparseArrays @@ -53,7 +54,8 @@ using ..ArrayUtils This dimension could be zero, in which case the delta set consists only of vertices (0-simplices). """ -@abstract_acset_type HasDeltaSet +@abstract_acset_type HasDeltaSet +const HasDeltaSet0D = HasDeltaSet vertices(s::HasDeltaSet) = parts(s, :V) nv(s::HasDeltaSet) = nparts(s, :V) @@ -84,7 +86,6 @@ end """ A 0-dimensional delta set, aka a set of vertices. """ @acset_type DeltaSet0D(SchDeltaSet0D) <: HasDeltaSet -ne(::DeltaSet0D) = error("0-dimensional Δ-set lacks edges.") # 1D simplicial sets #################### @@ -111,10 +112,12 @@ More generally, this type implements the graphs interface in `Catlab.Graphs`. """ @acset_type DeltaSet1D(SchDeltaSet1D, index=[:∂v0,:∂v1]) <: AbstractDeltaSet1D +edges(::HasDeltaSet) = error("0D simplicial sets have no edges") edges(s::HasDeltaSet1D) = parts(s, :E) edges(s::HasDeltaSet1D, src::Int, tgt::Int) = (e for e in coface(1,1,s,src) if ∂(1,0,s,e) == tgt) +ne(::HasDeltaSet) = 0 ne(s::HasDeltaSet1D) = nparts(s, :E) nsimplices(::Type{Val{1}}, s::HasDeltaSet1D) = ne(s) @@ -185,7 +188,7 @@ function d_nz(::Type{Val{0}}, s::HasDeltaSet1D, v::Int) (lazy(vcat, e₀, e₁), lazy(vcat, sign(1,s,e₀), -sign(1,s,e₁))) end -# 1D embedded simplicial sets +# 1D embedded, oriented simplicial sets #---------------------------- @present SchEmbeddedDeltaSet1D <: SchOrientedDeltaSet1D begin @@ -275,6 +278,7 @@ function triangles(s::HasDeltaSet2D, v₀::Int, v₁::Int, v₂::Int) end ntriangles(s::HasDeltaSet2D) = nparts(s, :Tri) +ntriangles(s::HasDeltaSet) = 0 nsimplices(::Type{Val{2}}, s::HasDeltaSet2D) = ntriangles(s) face(::Type{Val{(2,0)}}, s::HasDeltaSet2D, args...) = subpart(s, args..., :∂e0) @@ -355,8 +359,32 @@ function glue_sorted_triangle!(s::HasDeltaSet2D, v₀::Int, v₁::Int, v₂::Int glue_triangle!(s, v₀, v₁, v₂; kw...) end +""" Subdivide a 1D delta set. Note that this is written as if it'll work +for any type of 1D delta-set, but it can't handle orientations or embeddings. """ -Designed for simplicial complexes. +function subdivide(s::HasDeltaSet1D) + @migrate typeof(s) s begin + V => @cases begin + v::V + e::E + end + E => @cases begin + e₁::E + e₂::E + end + ∂v1 => begin + e₁ => e + e₂ => e + end + ∂v0 => begin + e₁ => (v∘∂v1) + e₂ => (v∘∂v0) + end + end +end +""" +Subdivision of a 2D simplicial set, relies on glue_triangle! so not good +for arbitrary simplicial sets. """ function subdivide(s::HasDeltaSet2D) d = typeof(s)() @@ -654,7 +682,16 @@ volume(::Type{Val{3}}, s::HasDeltaSet3D, t::Int, ::CayleyMengerDet) = # General operators ################### -DeltaSet(n) = eval(Symbol("DeltaSet$(n)D")) +DeltaSetTypes = Dict{Tuple{Symbol,Int},Type}() +add_type!(s,n) = DeltaSetTypes[(s,n)] = eval(Symbol(string(s)*string(n)*"D")) +for symb in [:DeltaSet,:EmbeddedDeltaSet,:OrientedDeltaSet] + for n in 1:3 + add_type!(symb,n) + end + #defines eg DeltaSet(2) = DeltaSet2D + eval(Expr(:(=),Expr(:call,symb,:n),Expr(:ref,:DeltaSetTypes,Expr(:tuple,QuoteNode(symb),:n)))) +end + """ Wrapper for simplex or simplices of dimension `n`. @@ -690,7 +727,11 @@ function simplex_vertices(::Type{Val{n}},s::HasDeltaSet,x::Simplex{n,0}) where n n == 3 && return tetrahedron_vertices(s, x.data) end -""" Wrapper for simplex chain of dimension `n`.""" +""" Wrapper for simplex chain of dimension `n`. + +Example: EChain([2,-1,1]) represents the chain 2a-b+c in the +simplicial set with edges a,b,c. +""" @vector_struct SimplexChain{n} const VChain = SimplexChain{0} From 40c2234ecf5ff9cde67ac439f4b74512882578ac Mon Sep 17 00:00:00 2001 From: Kevin Carlson Date: Fri, 18 Oct 2024 15:34:58 -0400 Subject: [PATCH 11/16] Implementing subdivision functions for all DeltaSet types. Changes "D_" to "dual_" in DualComplex schemas. There's some handedness bug still in the comparison between this and the `subdivide` function. --- src/DiscreteExteriorCalculus.jl | 203 +++++++++++++++++++++---------- src/FastDEC.jl | 48 +++++--- src/SimplicialComplexes.jl | 3 + src/SimplicialSets.jl | 4 +- test/DiscreteExteriorCalculus.jl | 57 +++++---- test/SimplicialComplexes.jl | 2 + 6 files changed, 202 insertions(+), 115 deletions(-) diff --git a/src/DiscreteExteriorCalculus.jl b/src/DiscreteExteriorCalculus.jl index b6a0e0a0..5fd6316d 100644 --- a/src/DiscreteExteriorCalculus.jl +++ b/src/DiscreteExteriorCalculus.jl @@ -18,6 +18,7 @@ export DualSimplex, DualV, DualE, DualTri, DualTet, DualChain, DualForm, AbstractDeltaDualComplex3D, DeltaDualComplex3D, SchDeltaDualComplex3D, OrientedDeltaDualComplex3D, SchOrientedDeltaDualComplex3D, EmbeddedDeltaDualComplex3D, SchEmbeddedDeltaDualComplex3D, + DeltaDualComplex, EmbeddedDeltaDualComplex, OrientedDeltaDualComplex, SimplexCenter, Barycenter, Circumcenter, Incenter, geometric_center, subsimplices, primal_vertex, elementary_duals, dual_boundary, dual_derivative, ⋆, hodge_star, inv_hodge_star, δ, codifferential, ∇², laplace_beltrami, Δ, laplace_de_rham, @@ -26,7 +27,7 @@ export DualSimplex, DualV, DualE, DualTri, DualTet, DualChain, DualForm, vertex_center, edge_center, triangle_center, tetrahedron_center, dual_tetrahedron_vertices, dual_triangle_vertices, dual_edge_vertices, dual_point, dual_volume, subdivide_duals!, DiagonalHodge, GeometricHodge, subdivide, PPSharp, AltPPSharp, DesbrunSharp, LLSDDSharp, de_sign, - ♭♯, ♭♯_mat, flat_sharp, flat_sharp_mat + ♭♯, ♭♯_mat, flat_sharp, flat_sharp_mat, dualize, dual_extractor, extract_dual import Base: ndims import Base: * @@ -42,6 +43,7 @@ const Point3D = SVector{3,Float64} using ACSets.DenseACSets: attrtype_type using Catlab, Catlab.CategoricalAlgebra.CSets using Catlab.CategoricalAlgebra.FinSets: deleteat +using Catlab.CategoricalAlgebra.FunctorialDataMigrations: DeltaMigration, migrate import Catlab.CategoricalAlgebra.CSets: ∧ import Catlab.Theories: Δ using DataMigrations: @migrate @@ -123,7 +125,7 @@ end @present SchDeltaDualComplex1D <: SchDeltaSet1D begin # Dual vertices and edges. (DualV, DualE)::Ob - (D_∂v0, D_∂v1)::Hom(DualE, DualV) + (dual_∂v0, dual_∂v1)::Hom(DualE, DualV) # Centers of primal simplices are dual vertices. vertex_center::Hom(V, DualV) @@ -133,10 +135,10 @@ end # # (∂v0_dual, ∂v1_dual)::Hom(E,DualE) # - # ∂v0_dual ⋅ D_∂v1 == ∂v0 ⋅ vertex_center - # ∂v1_dual ⋅ D_∂v1 == ∂v1 ⋅ vertex_center - # ∂v0_dual ⋅ D_∂v0 == edge_center - # ∂v1_dual ⋅ D_∂v0 == edge_center + # ∂v0_dual ⋅ dual_∂v1 == ∂v0 ⋅ vertex_center + # ∂v1_dual ⋅ dual_∂v1 == ∂v1 ⋅ vertex_center + # ∂v0_dual ⋅ dual_∂v0 == edge_center + # ∂v1_dual ⋅ dual_∂v0 == edge_center # # We could, and arguably should, track these through dedicated morphisms, as # in the commented code above. We don't because it scales badly in dimension: @@ -157,7 +159,7 @@ The data structure includes both the primal complex and the dual complex, as well as the mapping between them. """ @acset_type DeltaDualComplex1D(SchDeltaDualComplex1D, - index=[:∂v0,:∂v1,:D_∂v0,:D_∂v1]) <: AbstractDeltaDualComplex1D + index=[:∂v0,:∂v1,:dual_∂v0,:dual_∂v1]) <: AbstractDeltaDualComplex1D """ Dual vertex corresponding to center of primal vertex. """ @@ -168,12 +170,12 @@ vertex_center(s::HasDeltaSet, args...) = s[args..., :vertex_center] edge_center(s::HasDeltaSet1D, args...) = s[args..., :edge_center] subsimplices(::Type{Val{1}}, s::HasDeltaSet1D, e::Int) = - SVector{2}(incident(s, edge_center(s, e), :D_∂v0)) + SVector{2}(incident(s, edge_center(s, e), :dual_∂v0)) -primal_vertex(::Type{Val{1}}, s::HasDeltaSet1D, e...) = s[e..., :D_∂v1] +primal_vertex(::Type{Val{1}}, s::HasDeltaSet1D, e...) = s[e..., :dual_∂v1] elementary_duals(::Type{Val{0}}, s::AbstractDeltaDualComplex1D, v::Int) = - incident(s, vertex_center(s,v), :D_∂v1) + incident(s, vertex_center(s,v), :dual_∂v1) elementary_duals(::Type{Val{1}}, s::AbstractDeltaDualComplex1D, e::Int) = SVector(edge_center(s,e)) @@ -182,8 +184,8 @@ elementary_duals(::Type{Val{1}}, s::AbstractDeltaDualComplex1D, e::Int) = This accessor assumes that the simplicial identities for the dual hold. """ function dual_edge_vertices(s::HasDeltaSet1D, t...) - SVector(s[t..., :D_∂v0], - s[t..., :D_∂v1]) + SVector(s[t..., :dual_∂v0], + s[t..., :dual_∂v1]) end @@ -191,10 +193,10 @@ end This accessor assumes that the simplicial identities for the dual hold. """ -function dual_triangle_vertices(s::HasDeltaSet2D, t...) - SVector(s[s[t..., :D_∂e1], :D_∂v1], - s[s[t..., :D_∂e0], :D_∂v1], - s[s[t..., :D_∂e0], :D_∂v0]) +function dual_triangle_vertices(s::HasDeltaSet1D, t...) + SVector(s[s[t..., :dual_∂e1], :dual_∂v1], + s[s[t..., :dual_∂e0], :dual_∂v1], + s[s[t..., :dual_∂e0], :dual_∂v0]) end @@ -206,13 +208,13 @@ end @present SchOrientedDeltaDualComplex1D <: SchDeltaDualComplex1D begin Orientation::AttrType edge_orientation::Attr(E, Orientation) - D_edge_orientation::Attr(DualE, Orientation) + dual_edge_orientation::Attr(DualE, Orientation) end """ Oriented dual complex of an oriented 1D delta set. """ @acset_type OrientedDeltaDualComplex1D(SchOrientedDeltaDualComplex1D, - index=[:∂v0,:∂v1,:D_∂v0,:D_∂v1]) <: AbstractDeltaDualComplex1D + index=[:∂v0,:∂v1,:dual_∂v0,:dual_∂v1]) <: AbstractDeltaDualComplex1D dual_boundary_nz(::Type{Val{1}}, s::AbstractDeltaDualComplex1D, x::Int) = # Boundary vertices of dual 1-cell ↔ @@ -275,9 +277,9 @@ function make_dual_simplices_1d!(s::HasDeltaSet1D, ::Type{Simplex{n}}) where n # Make dual vertices and edges. s[:vertex_center] = vcenters = add_parts!(s, :DualV, nv(s)) s[:edge_center] = ecenters = add_parts!(s, :DualV, ne(s)) - D_edges = map((0,1)) do i + dual_edges = map((0,1)) do i add_parts!(s, :DualE, ne(s); - D_∂v0 = ecenters, D_∂v1 = view(vcenters, ∂(1,i,s))) + dual_∂v0 = ecenters, dual_∂v1 = view(vcenters, ∂(1,i,s))) end # Orient elementary dual edges. @@ -293,11 +295,11 @@ function make_dual_simplices_1d!(s::HasDeltaSet1D, ::Type{Simplex{n}}) where n end end edge_orient = s[:edge_orientation] - s[D_edges[1], :D_edge_orientation] = negate.(edge_orient) - s[D_edges[2], :D_edge_orientation] = edge_orient + s[dual_edges[1], :dual_edge_orientation] = negate.(edge_orient) + s[dual_edges[2], :dual_edge_orientation] = edge_orient end - D_edges + dual_edges end # TODO: Instead of copying-and-pasting the DeltaSet1D version: @@ -352,7 +354,7 @@ Although they are redundant information, the lengths of the primal and dual edges are precomputed and stored. """ @acset_type EmbeddedDeltaDualComplex1D(SchEmbeddedDeltaDualComplex1D, - index=[:∂v0,:∂v1,:D_∂v0,:D_∂v1]) <: AbstractDeltaDualComplex1D + index=[:∂v0,:∂v1,:dual_∂v0,:dual_∂v1]) <: AbstractDeltaDualComplex1D """ Point associated with dual vertex of complex. """ @@ -370,7 +372,7 @@ dual_volume(::Type{Val{1}}, s::HasDeltaSet1D, e, ::PrecomputedVol) = s[e, :dual_length] dual_volume(::Type{Val{1}}, s::HasDeltaSet1D, e::Int, ::CayleyMengerDet) = - volume(dual_point(s, SVector(s[e,:D_∂v0], s[e,:D_∂v1]))) + volume(dual_point(s, SVector(s[e,:dual_∂v0], s[e,:dual_∂v1]))) hodge_diag(::Type{Val{0}}, s::AbstractDeltaDualComplex1D, v::Int) = sum(dual_volume(Val{1}, s, elementary_duals(Val{0},s,v))) @@ -449,13 +451,13 @@ end @present SchDeltaDualComplex2D <: SchDeltaSet2D begin # Dual vertices, edges, and triangles. (DualV, DualE, DualTri)::Ob - (D_∂v0, D_∂v1)::Hom(DualE, DualV) - (D_∂e0, D_∂e1, D_∂e2)::Hom(DualTri, DualE) + (dual_∂v0, dual_∂v1)::Hom(DualE, DualV) + (dual_∂e0, dual_∂e1, dual_∂e2)::Hom(DualTri, DualE) # Simplicial identities for dual simplices. - D_∂e1 ⋅ D_∂v1 == D_∂e2 ⋅ D_∂v1 - D_∂e0 ⋅ D_∂v1 == D_∂e2 ⋅ D_∂v0 - D_∂e0 ⋅ D_∂v0 == D_∂e1 ⋅ D_∂v0 + dual_∂e1 ⋅ dual_∂v1 == dual_∂e2 ⋅ dual_∂v1 + dual_∂e0 ⋅ dual_∂v1 == dual_∂e2 ⋅ dual_∂v0 + dual_∂e0 ⋅ dual_∂v0 == dual_∂e1 ⋅ dual_∂v0 # Centers of primal simplices are dual vertices. vertex_center::Hom(V, DualV) @@ -476,25 +478,26 @@ end """ @abstract_acset_type AbstractDeltaDualComplex2D <: HasDeltaSet2D +const AbstractDeltaDualComplex = Union{AbstractDeltaDualComplex1D, AbstractDeltaDualComplex2D} """ Dual complex of a two-dimensional delta set. """ @acset_type DeltaDualComplex2D(SchDeltaDualComplex2D, - index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:D_∂v0,:D_∂v1,:D_∂e0,:D_∂e1,:D_∂e2]) <: AbstractDeltaDualComplex2D + index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:dual_∂v0,:dual_∂v1,:dual_∂e0,:dual_∂e1,:dual_∂e2]) <: AbstractDeltaDualComplex2D """ Dual vertex corresponding to center of primal triangle. """ triangle_center(s::HasDeltaSet2D, args...) = s[args..., :tri_center] subsimplices(::Type{Val{2}}, s::HasDeltaSet2D, t::Int) = - SVector{6}(incident(s, triangle_center(s,t), @SVector [:D_∂e1, :D_∂v0])) + SVector{6}(incident(s, triangle_center(s,t), @SVector [:dual_∂e1, :dual_∂v0])) primal_vertex(::Type{Val{2}}, s::HasDeltaSet2D, t...) = - primal_vertex(Val{1}, s, s[t..., :D_∂e2]) + primal_vertex(Val{1}, s, s[t..., :dual_∂e2]) elementary_duals(::Type{Val{0}}, s::AbstractDeltaDualComplex2D, v::Int) = - incident(s, vertex_center(s,v), @SVector [:D_∂e1, :D_∂v1]) + incident(s, vertex_center(s,v), @SVector [:dual_∂e1, :dual_∂v1]) elementary_duals(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, e::Int) = - incident(s, edge_center(s,e), :D_∂v1) + incident(s, edge_center(s,e), :dual_∂v1) elementary_duals(::Type{Val{2}}, s::AbstractDeltaDualComplex2D, t::Int) = SVector(triangle_center(s,t)) @@ -505,14 +508,14 @@ elementary_duals(::Type{Val{2}}, s::AbstractDeltaDualComplex2D, t::Int) = Orientation::AttrType edge_orientation::Attr(E, Orientation) tri_orientation::Attr(Tri, Orientation) - D_edge_orientation::Attr(DualE, Orientation) - D_tri_orientation::Attr(DualTri, Orientation) + dual_edge_orientation::Attr(DualE, Orientation) + dual_tri_orientation::Attr(DualTri, Orientation) end """ Oriented dual complex of an oriented 2D delta set. """ @acset_type OrientedDeltaDualComplex2D(SchOrientedDeltaDualComplex2D, - index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:D_∂v0,:D_∂v1,:D_∂e0,:D_∂e1,:D_∂e2]) <: AbstractDeltaDualComplex2D + index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:dual_∂v0,:dual_∂v1,:dual_∂e0,:dual_∂e1,:dual_∂e2]) <: AbstractDeltaDualComplex2D dual_boundary_nz(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, x::Int) = # Boundary vertices of dual 1-cell ↔ @@ -570,24 +573,24 @@ relative to the primal triangles they subdivide. """ function make_dual_simplices_2d!(s::HasDeltaSet2D, ::Type{Simplex{n}}) where n # Make dual vertices and edges. - D_edges01 = make_dual_simplices_1d!(s) + dual_edges01 = make_dual_simplices_1d!(s) s[:tri_center] = tri_centers = add_parts!(s, :DualV, ntriangles(s)) - D_edges12 = map((0,1,2)) do e + dual_edges12 = map((0,1,2)) do e add_parts!(s, :DualE, ntriangles(s); - D_∂v0=tri_centers, D_∂v1=edge_center(s, ∂(2,e,s))) + dual_∂v0=tri_centers, dual_∂v1=edge_center(s, ∂(2,e,s))) end - D_edges02 = map(triangle_vertices(s)) do vs + dual_edges02 = map(triangle_vertices(s)) do vs add_parts!(s, :DualE, ntriangles(s); - D_∂v0=tri_centers, D_∂v1=vertex_center(s, vs)) + dual_∂v0=tri_centers, dual_∂v1=vertex_center(s, vs)) end # Make dual triangles. # Counterclockwise order in drawing with vertices 0, 1, 2 from left to right. - D_triangle_schemas = ((0,1,1),(0,2,1),(1,2,0),(1,0,1),(2,0,0),(2,1,0)) - D_triangles = map(D_triangle_schemas) do (v,e,ev) + dual_triangle_schemas = ((0,1,1),(0,2,1),(1,2,0),(1,0,1),(2,0,0),(2,1,0)) + dual_triangles = map(dual_triangle_schemas) do (v,e,ev) add_parts!(s, :DualTri, ntriangles(s); - D_∂e0=D_edges12[e+1], D_∂e1=D_edges02[v+1], - D_∂e2=view(D_edges01[ev+1], ∂(2,e,s))) + dual_∂e0=dual_edges12[e+1], dual_∂e1=dual_edges02[v+1], + dual_∂e2=view(dual_edges01[ev+1], ∂(2,e,s))) end if has_subpart(s, :tri_orientation) @@ -604,21 +607,21 @@ function make_dual_simplices_2d!(s::HasDeltaSet2D, ::Type{Simplex{n}}) where n # Orient elementary dual triangles. tri_orient = s[:tri_orientation] rev_tri_orient = negate.(tri_orient) - for (i, D_tris) in enumerate(D_triangles) - s[D_tris, :D_tri_orientation] = isodd(i) ? rev_tri_orient : tri_orient + for (i, dual_tris) in enumerate(dual_triangles) + s[dual_tris, :dual_tri_orientation] = isodd(i) ? rev_tri_orient : tri_orient end # Orient elementary dual edges. for e in (0,1,2) - s[D_edges12[e+1], :D_edge_orientation] = relative_sign.( + s[dual_edges12[e+1], :dual_edge_orientation] = relative_sign.( s[∂(2,e,s), :edge_orientation], isodd(e) ? rev_tri_orient : tri_orient) end # Remaining dual edges are oriented arbitrarily. - s[lazy(vcat, D_edges02...), :D_edge_orientation] = one(attrtype_type(s, :Orientation)) + s[lazy(vcat, dual_edges02...), :dual_edge_orientation] = one(attrtype_type(s, :Orientation)) end - D_triangles + dual_triangles end relative_sign(x, y) = sign(x*y) @@ -643,7 +646,7 @@ Although they are redundant information, the lengths and areas of the primal/dual edges and triangles are precomputed and stored. """ @acset_type EmbeddedDeltaDualComplex2D(SchEmbeddedDeltaDualComplex2D, - index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:D_∂v0,:D_∂v1,:D_∂e0,:D_∂e1,:D_∂e2]) <: AbstractDeltaDualComplex2D + index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:dual_∂v0,:dual_∂v1,:dual_∂e0,:dual_∂e1,:dual_∂e2]) <: AbstractDeltaDualComplex2D volume(::Type{Val{n}}, s::EmbeddedDeltaDualComplex2D, x) where n = volume(Val{n}, s, x, PrecomputedVol()) @@ -655,9 +658,9 @@ dual_volume(::Type{Val{2}}, s::HasDeltaSet2D, t, ::PrecomputedVol) = s[t, :dual_area] function dual_volume(::Type{Val{2}}, s::HasDeltaSet2D, t::Int, ::CayleyMengerDet) - dual_vs = SVector(s[s[t, :D_∂e1], :D_∂v1], - s[s[t, :D_∂e2], :D_∂v0], - s[s[t, :D_∂e0], :D_∂v0]) + dual_vs = SVector(s[s[t, :dual_∂e1], :dual_∂v1], + s[s[t, :dual_∂e2], :dual_∂v0], + s[s[t, :dual_∂e0], :dual_∂v0]) volume(dual_point(s, dual_vs)) end @@ -687,7 +690,7 @@ function ♭(s::AbstractDeltaDualComplex2D, X::AbstractVector, ::DPPFlat) # For each of these dual edges: mapreduce(+, dual_edges, dual_lengths) do dual_e, dual_length # Get the vector at the center of the triangle this edge is pointing at. - X_vec = X[tri_map[s[dual_e, :D_∂v0]]] + X_vec = X[tri_map[s[dual_e, :dual_∂v0]]] # Take their dot product and multiply by the length of this dual edge. dual_length * dot(X_vec, e_vec) # When done, sum these weights up and divide by the total length. @@ -714,7 +717,7 @@ function ♭_mat(s::AbstractDeltaDualComplex2D, p2s) # The dual edge pointing to each triangle. des = map(dvs) do dv # (This is the edges(s,src,tgt) function.) - only(de for de in incident(s, dv, :D_∂v0) if s[de, :D_∂v1] == center) + only(de for de in incident(s, dv, :dual_∂v0) if s[de, :dual_∂v1] == center) end # The lengths of those dual edges. dels = volume(s, DualE(des)) @@ -743,7 +746,7 @@ function ♯(s::AbstractDeltaDualComplex2D, α::AbstractVector, DS::DiscreteShar for e in deleteat(tri_edges, i) v, sgn = src(s,e) == v₀ ? (tgt(s,e), -1) : (src(s,e), +1) dual_area = sum(dual_volume(2,s,d) for d in elementary_duals(0,s,v) - if s[s[d, :D_∂e0], :D_∂v0] == tri_center) + if s[s[d, :dual_∂e0], :dual_∂v0] == tri_center) area = ♯_denominator(s, v, t, DS) α♯[v] += sgn * sign(1,s,e) * α[e] * (dual_area / area) * out_vec end @@ -796,7 +799,7 @@ function ♯_assign!(♯_mat::AbstractSparseMatrix, s::AbstractDeltaDualComplex2 v, sgn = src(s,e) == v₀ ? (tgt(s,e), -1) : (src(s,e), +1) # | ⋆vₓ ∩ σⁿ | dual_area = sum(dual_volume(2,s,d) for d in elementary_duals(0,s,v) - if s[s[d, :D_∂e0], :D_∂v0] == tri_center) + if s[s[d, :dual_∂e0], :dual_∂v0] == tri_center) area = ♯_denominator(s, v, t, DS) ♯_mat[v,e] += sgn * sign(1,s,e) * (dual_area / area) * out_vec end @@ -809,7 +812,7 @@ function ♯_assign!(♯_mat::AbstractSparseMatrix, s::AbstractDeltaDualComplex2 sgn = v == tgt(s,e₀) ? -1 : +1 # | ⋆vₓ ∩ σⁿ | dual_area = sum(dual_volume(2,s,d) for d in elementary_duals(0,s,v) - if s[s[d, :D_∂e0], :D_∂v0] == tri_center) + if s[s[d, :dual_∂e0], :dual_∂v0] == tri_center) area = ♯_denominator(s, v, t, DS) ♯_mat[v,e₀] += sgn * sign(1,s,e₀) * (dual_area / area) * out_vec end @@ -837,7 +840,7 @@ function ♯_mat(s::AbstractDeltaDualComplex2D, DS::DiscreteSharp) ♯_mat end -de_sign(s,de) = s[de, :D_edge_orientation] ? +1 : -1 +de_sign(s,de) = s[de, :dual_edge_orientation] ? +1 : -1 """ function ♯_mat(s::AbstractDeltaDualComplex2D, ::LLSDDSharp) @@ -857,12 +860,12 @@ function ♯_mat(s::AbstractDeltaDualComplex2D, ::LLSDDSharp) # | ⋆eₓ ∩ σⁿ | star_e_cap_t = map(tri_edges) do e only(filter(elementary_duals(1,s,e)) do de - s[de, :D_∂v0] == tri_center + s[de, :dual_∂v0] == tri_center end) end de_vecs = map(star_e_cap_t) do de de_sign(s,de) * - (dual_point(s,s[de, :D_∂v0]) - dual_point(s,s[de, :D_∂v1])) + (dual_point(s,s[de, :dual_∂v0]) - dual_point(s,s[de, :dual_∂v1])) end weights = s[star_e_cap_t, :dual_length] ./ map(tri_edges) do e @@ -886,10 +889,10 @@ function ∧(::Type{Tuple{1,1}}, s::HasDeltaSet2D, α, β, x::Int) # XXX: This calculation of the volume coefficients is awkward due to the # design decision described in `SchDeltaDualComplex1D`. dual_vs = vertex_center(s, triangle_vertices(s, x)) - dual_es = sort(SVector{6}(incident(s, triangle_center(s, x), :D_∂v0)), - by=e -> s[e,:D_∂v1] .== dual_vs, rev=true)[1:3] + dual_es = sort(SVector{6}(incident(s, triangle_center(s, x), :dual_∂v0)), + by=e -> s[e,:dual_∂v1] .== dual_vs, rev=true)[1:3] coeffs = map(dual_es) do e - sum(dual_volume(2, s, SVector{2}(incident(s, e, :D_∂e1)))) + sum(dual_volume(2, s, SVector{2}(incident(s, e, :dual_∂e1)))) end / volume(2, s, x) # Wedge product of two primal 1-forms, as in (Hirani 2003, Example 7.1.2). @@ -1466,6 +1469,72 @@ function dual_derivative(::Type{Val{n}}, s::HasDeltaSet, args...) where n end end +""" +Checks whethere a DeltaSet is embedded by (droolingly) searching for 'Embedded' +in the name of its type. This could also check for 'Point' in the schema, which +would feel better but be less trustworthy. +""" +is_embedded(d::HasDeltaSet) = is_embedded(typeof(t)) +is_embedded(t::Type{T}) where {T<:HasDeltaSet} = !isnothing(findfirst("Embedded",string(t.name.name))) +const REPLACEMENT_FOR_DUAL_TYPE = "Set" => "DualComplex" +rename_to_dual(s::Symbol) = Symbol(replace(string(s),REPLACEMENT_FOR_DUAL_TYPE)) +rename_from_dual(s::Symbol) = Symbol(replace(string(s),reverse(REPLACEMENT_FOR_DUAL_TYPE))) + +#adds the Real type for lengths in the Embedded case. Will need further customization +#if we add another type whose dual has more parameters than its primal. +function dual_param_list(t) + is_embedded(t) ? [t.parameters[1],eltype(t.parameters[2]),t.parameters[2]] : t.parameters +end +""" +Calls the constructor for d's dual type on d, including parameters. +Does not call `subdivide_duals!` on the result. +Should work out of the box on new DeltaSet types if (1) their dual type +has the same name as their primal type with "Set" substituted by "DualComplex" +and (2) their dual type has the same parameter set as their primal type. At the +time of writing (July 2024) only "Embedded" types fail criterion (2) and get special treatment. + +# Examples +s = EmbeddedDeltaSet2D{Bool,SVector{Float64,Float64}}() +dualize(s)::EmbeddedDeltaDualComplex2D{Bool,Float64,SVector{Float64,Float64}} +""" +function dualize(d::HasDeltaSet) #add embedded keyword? + t = typeof(d) + n = dualize_type_name(t) + ps = dual_param_list(t) + length(ps) > 0 ? n{ps...}(d) : n(d) +end +#Note this can only successfully dualize types whose duals are defined in the scope of this module! +dualize_type_name(t::Type) = eval(rename_to_dual(t.name.name)) + +dual_extractor(d::HasDeltaSet) = dual_extractor(typeof(d)) +function dual_extractor(t::Type{T}) where {T <: HasDeltaSet} + n = t.name.name + sch_name = Symbol("Sch",n) + sch = eval(sch_name) + dual_sch = eval(rename_to_dual(sch_name)) + C = FinCat(sch) + D = FinCat(dual_sch) + o = Dict(nameof(a) => Symbol("Dual",nameof(a)) for a in sch.generators[:Ob]) + for a in sch.generators[:AttrType] o[nameof(a)] = nameof(a) end + h = Dict(nameof(a) => Symbol("dual_",nameof(a)) for a in hom_generators(C)) + DeltaMigration(FinDomFunctor(o,h,C,D)) +end +dual_extractor(t::Type{T}) where {T <: AbstractDeltaDualComplex} = dual_extractor(eval(rename_from_dual(t.name.name))) + +""" +Get the subdivision of a `DeltaSet` `d` as a `DeltaSet` of the same type, +or extract just the subdivided part of a `DeltaDualComplex` as a `DeltaSet`. +""" +extract_dual(d::HasDeltaSet) = migrate(typeof(d),dualize(d),dual_extractor(d)) +function extract_dual(d::HasDeltaSet,alg) + is_embedded(d) || error("Cannot subdivide duals for a non-embedded DeltaSet.") + s = dualize(d) + subdivide_duals!(s,alg) + migrate(typeof(d),s,dual_extractor(d)) +end +extract_dual(d::AbstractDeltaDualComplex) = migrate(eval(rename_from_dual(typeof(d).name.name)),d,dual_extractor(d)) + + """ Hodge star operator from primal ``n``-forms to dual ``N-n``-forms. !!! note diff --git a/src/FastDEC.jl b/src/FastDEC.jl index 75eaa818..5f704444 100644 --- a/src/FastDEC.jl +++ b/src/FastDEC.jl @@ -186,7 +186,7 @@ function wedge_dd_01_mat(sd::HasDeltaSet) m = spzeros(ne(sd), ntriangles(sd)) for e in edges(sd) des = elementary_duals(1,sd,e) - dvs = sd[des, :D_∂v0] + dvs = sd[des, :dual_∂v0] tris = only.(incident(sd, dvs, :tri_center)) ws = sd[des, :dual_length] ./ sum(sd[des, :dual_length]) for (w,t) in zip(ws,tris) @@ -223,7 +223,7 @@ function wedge_pd_01_mat(sd::HasDeltaSet) for e in edges(sd) α, β = edge_vertices(sd,e) des = elementary_duals(1,sd,e) - dvs = sd[des, :D_∂v0] + dvs = sd[des, :dual_∂v0] tris = only.(incident(sd, dvs, :tri_center)) γδ = map(tris) do t only(filter(x -> x ∉ [α,β], triangle_vertices(sd,t))) @@ -398,12 +398,15 @@ end #-------------------- function dec_p_hodge_diag(::Type{Val{0}}, sd::EmbeddedDeltaDualComplex1D{Bool, float_type, _p} where _p) where float_type - nvsd = nv(sd) - h_0 = zeros(float_type, nvsd) - for de in parts(sd, :DualE) - v1 = sd[de, :D_∂v1] - if 1 <= v1 <= nvsd - h_0[v1] += sd[de, :dual_length] + num_v_sd = nv(sd) + + hodge_diag_0 = zeros(float_type, num_v_sd) + + for d_edge_idx in parts(sd, :DualE) + v1 = sd[d_edge_idx, :dual_∂v1] + if (1 <= v1 <= num_v_sd) + hodge_diag_0[v1] += sd[d_edge_idx, :dual_length] + end end end h_0 @@ -416,21 +419,26 @@ end function dec_p_hodge_diag(::Type{Val{0}}, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, _p} where _p) where float_type - h_0 = zeros(float_type, nv(sd)) - for dt in parts(sd, :DualTri) - v = sd[sd[dt, :D_∂e1], :D_∂v1] - h_0[v] += sd[dt, :dual_area] - end - h_0 + hodge_diag_0 = zeros(float_type, nv(sd)) + + for dual_tri in parts(sd, :DualTri) + v = sd[sd[dual_tri, :dual_∂e1], :dual_∂v1] + hodge_diag_0[v] += sd[dual_tri, :dual_area] + end + return hodge_diag_0 end function dec_p_hodge_diag(::Type{Val{1}}, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, _p} where _p) where float_type - nvsd, nesd = nv(sd), ne(sd) - h_1 = zeros(float_type, nesd) - for de in parts(sd, :DualE) - v1_shift = sd[de, :D_∂v1] - nvsd - if (1 <= v1_shift <= nesd) - h_1[v1_shift] += sd[de, :dual_length] / sd[v1_shift, :length] + num_v_sd = nv(sd) + num_e_sd = ne(sd) + + hodge_diag_1 = zeros(float_type, num_e_sd) + + for d_edge_idx in parts(sd, :DualE) + v1_shift = sd[d_edge_idx, :dual_∂v1] - num_v_sd + if (1 <= v1_shift <= num_e_sd) + hodge_diag_1[v1_shift] += sd[d_edge_idx, :dual_length] / sd[v1_shift, :length] + end end end h_1 diff --git a/src/SimplicialComplexes.jl b/src/SimplicialComplexes.jl index 618f72ec..e97affe0 100644 --- a/src/SimplicialComplexes.jl +++ b/src/SimplicialComplexes.jl @@ -194,6 +194,7 @@ end #XX: This type is maybe more trouble than it's worth? +#yeah kill """ A point in an unspecified simplicial complex, given by its barycentric coordinates. Constructed via a dense vector of coordinates. @@ -246,6 +247,7 @@ as_matrix(f::GeometricMap) = hcat(coords.(f.points)...) compose(f::GeometricMap, g::GeometricMap) = GeometricMap(f.dom, g.cod, as_matrix(g)*as_matrix(f)) id(sc::SimplicialComplex) = GeometricMap(sc,sc,GeometricPoint.(eachcol(I(nv(sc))))) +#make sparse matrices function GeometricMap(sc::SimplicialComplex,::Barycenter) dom = SimplicialComplex(subdivide(sc.delta_set)) #Vertices of dom correspond to vertices, edges, triangles of sc. @@ -255,6 +257,7 @@ function GeometricMap(sc::SimplicialComplex,::Barycenter) for i in 1:ntriangles(sc) for n in triangle_vertices(sc.delta_set,i) mat[n,nv(sc)+ne(sc)+i] = 1/3 end end GeometricMap(dom,sc,mat) end +#accessors for the nonzeros in a column of the matrix function pullback_primal(f::GeometricMap, v::PrimalVectorField{T}) where T nv(f.cod) == length(v) || error("Vector field must have same number of vertices as codomain") diff --git a/src/SimplicialSets.jl b/src/SimplicialSets.jl index 53b30bd4..055e07b8 100644 --- a/src/SimplicialSets.jl +++ b/src/SimplicialSets.jl @@ -46,6 +46,7 @@ using ACSets.DenseACSets: attrtype_type using Catlab, Catlab.CategoricalAlgebra, Catlab.Graphs import Catlab.Graphs: src, tgt, nv, ne, vertices, edges, has_vertex, has_edge, add_vertex!, add_vertices!, add_edge!, add_edges! +import DataMigrations: @migrate using ..ArrayUtils @@ -112,7 +113,7 @@ More generally, this type implements the graphs interface in `Catlab.Graphs`. """ @acset_type DeltaSet1D(SchDeltaSet1D, index=[:∂v0,:∂v1]) <: AbstractDeltaSet1D -edges(::HasDeltaSet) = error("0D simplicial sets have no edges") +edges(::HasDeltaSet) = error("0-dimensional simplicial sets have no edges") edges(s::HasDeltaSet1D) = parts(s, :E) edges(s::HasDeltaSet1D, src::Int, tgt::Int) = (e for e in coface(1,1,s,src) if ∂(1,0,s,e) == tgt) @@ -684,6 +685,7 @@ volume(::Type{Val{3}}, s::HasDeltaSet3D, t::Int, ::CayleyMengerDet) = DeltaSetTypes = Dict{Tuple{Symbol,Int},Type}() add_type!(s,n) = DeltaSetTypes[(s,n)] = eval(Symbol(string(s)*string(n)*"D")) +add_type!(:DeltaSet,0) for symb in [:DeltaSet,:EmbeddedDeltaSet,:OrientedDeltaSet] for n in 1:3 add_type!(symb,n) diff --git a/test/DiscreteExteriorCalculus.jl b/test/DiscreteExteriorCalculus.jl index 50b0d438..bda557f1 100644 --- a/test/DiscreteExteriorCalculus.jl +++ b/test/DiscreteExteriorCalculus.jl @@ -22,7 +22,8 @@ const Point3D = SVector{3,Float64} primal_s = DeltaSet1D() add_vertices!(primal_s, 5) add_edges!(primal_s, 1:4, repeat([5], 4)) -s = DeltaDualComplex1D(primal_s) +s = dualize(primal_s) +@test extract_dual(primal_s) == extract_dual(s) @test nparts(s, :DualV) == nv(primal_s) + ne(primal_s) @test nparts(s, :DualE) == 2 * ne(primal_s) @@ -32,10 +33,11 @@ dual_v = elementary_duals(1,s,4) dual_es = elementary_duals(0,s,5) @test length(dual_es) == 4 -@test s[dual_es, :D_∂v0] == edge_center(s, 1:4) +@test s[dual_es, :dual_∂v0] == edge_center(s, 1:4) @test elementary_duals(s, V(5)) == DualE(dual_es) primal_s′ = subdivide(primal_s) +@test is_isomorphic(primal_s′,extract_dual(s)) #why is this failing? @test nv(primal_s′) == nv(primal_s) + ne(primal_s) @test ne(primal_s′) == 2*ne(primal_s) @@ -45,9 +47,9 @@ primal_s′ = subdivide(primal_s) primal_s = OrientedDeltaSet1D{Bool}() add_vertices!(primal_s, 3) add_edges!(primal_s, [1,2], [2,3], edge_orientation=[true,false]) -s = OrientedDeltaDualComplex1D{Bool}(primal_s) -@test s[only(elementary_duals(0,s,1)), :D_edge_orientation] == true -@test s[only(elementary_duals(0,s,3)), :D_edge_orientation] == true +s = dualize(primal_s) +@test s[only(elementary_duals(0,s,1)), :dual_edge_orientation] == true +@test s[only(elementary_duals(0,s,3)), :dual_edge_orientation] == true @test ∂(s, DualChain{1}([1,0,1])) isa DualChain{0} @test d(s, DualForm{0}([1,1])) isa DualForm{1} @@ -73,7 +75,7 @@ add_vertices!(implicit_s, 3, point=[Point2D(1,0), Point2D(0,0), Point2D(0,2)]) add_edges!(implicit_s, [1,2], [2,3]) for primal_s in [explicit_s, implicit_s] - s = EmbeddedDeltaDualComplex1D{Bool,Float64,Point2D}(primal_s) + s = dualize(primal_s) subdivide_duals!(s, Barycenter()) @test dual_point(s, edge_center(s, [1,2])) ≈ [Point2D(0.5,0), Point2D(0,1)] @test volume(s, E(1:2)) ≈ [1.0, 2.0] @@ -97,7 +99,7 @@ end primal_s = EmbeddedDeltaSet1D{Bool,Point2D}() add_vertices!(primal_s, 5, point=[Point2D(i,0) for i in -2:2]) add_edges!(primal_s, 1:4, 2:5, edge_orientation=true) -s = EmbeddedDeltaDualComplex1D{Bool,Float64,Point2D}(primal_s) +s = dualize(primal_s) subdivide_duals!(s, Barycenter()) @test ∇²(s, VForm([0,0,1,0,0])) ≈ VForm([0,-1,2,-1,0]) @test ∇²(0,s) ≈ [ 2 -2 0 0 0; @@ -129,7 +131,7 @@ primal_s = DeltaSet2D() add_vertices!(primal_s, 4) glue_triangle!(primal_s, 1, 2, 3) glue_triangle!(primal_s, 1, 3, 4) -s = DeltaDualComplex2D(primal_s) +s = dualize(primal_s) @test nparts(s, :DualV) == nv(primal_s) + ne(primal_s) + ntriangles(primal_s) @test nparts(s, :DualE) == 2*ne(primal_s) + 6*ntriangles(primal_s) @test nparts(s, :DualTri) == 6*ntriangles(primal_s) @@ -138,8 +140,8 @@ s = DeltaDualComplex2D(primal_s) dual_vs = elementary_duals(2,s,2) @test dual_vs == [triangle_center(s,2)] @test elementary_duals(s, Tri(2)) == DualV(dual_vs) -@test s[elementary_duals(1,s,2), :D_∂v1] == [edge_center(s,2)] -@test s[elementary_duals(1,s,3), :D_∂v1] == repeat([edge_center(s,3)], 2) +@test s[elementary_duals(1,s,2), :dual_∂v1] == [edge_center(s,2)] +@test s[elementary_duals(1,s,3), :dual_∂v1] == repeat([edge_center(s,3)], 2) @test [length(elementary_duals(s, V(i))) for i in 1:4] == [4,2,4,2] @test dual_triangle_vertices(s, 1) == [1,7,10] @test dual_edge_vertices(s, 1) == [5,2] @@ -161,11 +163,11 @@ glue_triangle!(implicit_s, 1, 2, 3) glue_triangle!(implicit_s, 1, 3, 4) for primal_s in [explicit_s, implicit_s] - s = OrientedDeltaDualComplex2D{Bool}(primal_s) - @test sum(s[:D_tri_orientation]) == nparts(s, :DualTri) ÷ 2 - @test [sum(s[elementary_duals(0,s,i), :D_tri_orientation]) + s = dualize(primal_s) + @test sum(s[:dual_tri_orientation]) == nparts(s, :DualTri) ÷ 2 + @test [sum(s[elementary_duals(0,s,i), :dual_tri_orientation]) for i in 1:4] == [2,1,2,1] - @test sum(s[elementary_duals(1,s,3), :D_edge_orientation]) == 1 + @test sum(s[elementary_duals(1,s,3), :dual_edge_orientation]) == 1 for k in 0:1 @test dual_boundary(2-k,s) == (-1)^k * ∂(k+1,s)' @@ -199,7 +201,7 @@ primal_s = get_regular_polygon(6) # Rotate counter-clockwise by pi/6 to match the Hirani figure. θ = -pi/6 primal_s[:point] = [[[cos(θ), -sin(θ), 0];; [sin(θ), cos(θ), 0];; [0,0,1]] * p for p in primal_s[:point]] -s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(primal_s) +s = dualize(primal_s) subdivide_duals!(s, Circumcenter()) X = map([8,4,0,8,4,0]) do i SVector(unit_vector(i*(2pi/12))) @@ -238,7 +240,7 @@ primal_s = EmbeddedDeltaSet2D{Bool,Point2D}() add_vertices!(primal_s, 3, point=[Point2D(0,0), Point2D(1,0), Point2D(0,1)]) glue_triangle!(primal_s, 1, 2, 3, tri_orientation=true) primal_s[:edge_orientation] = true -s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(primal_s) +s = dualize(primal_s) subdivide_duals!(s, Barycenter()) @test dual_point(s, triangle_center(s, 1)) ≈ Point2D(1/3, 1/3) @@ -263,7 +265,7 @@ subdivide_duals!(s, Barycenter()) # geometric hodge star) flipped_ps = deepcopy(primal_s) orient_component!(flipped_ps, 1, false) -flipped_s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(flipped_ps) +flipped_s = dualize(flipped_ps) subdivide_duals!(flipped_s, Barycenter()) @test ⋆(1,s) ≈ ⋆(1,flipped_s) @@ -357,6 +359,7 @@ function test_♯(s, covector::SVector; atol=1e-8) X♯ = ♯(s, X, AltPPSharp()) @test all(isapprox.(X♯, [covector])) # Test that the matrix and non-matrix versions yield the same result. + # XXX: Why do all these primal forms come out constant?! @test all(isapprox.(♯_mat(s, PPSharp()) * X, ♯(s, X, PPSharp()))) @test all(isapprox.(♯_mat(s, AltPPSharp()) * X, ♯(s, X, AltPPSharp()))) end @@ -382,7 +385,7 @@ foreach(vf -> test_♯(s, vf), vfs) # Triangulated regular dodecagon. primal_s = get_regular_polygon(12) primal_s[:point] = [Point3D(1/4,1/5,0) + p for p in primal_s[:point]] -s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(primal_s) +s = dualize(primal_s) subdivide_duals!(s, Circumcenter()) foreach(vf -> test_♯(s, vf), vfs) # TODO: Compute results for Desbrun's ♯ by hand. @@ -394,7 +397,7 @@ add_vertices!(primal_s, 4, point=[Point2D(-1,+1), Point2D(+1,+1), glue_triangle!(primal_s, 1, 2, 3, tri_orientation=true) glue_triangle!(primal_s, 1, 3, 4, tri_orientation=true) primal_s[:edge_orientation] = true -s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(primal_s) +s = dualize(primal_s) subdivide_duals!(s, Barycenter()) ♭_m = ♭_mat(s) @@ -443,7 +446,7 @@ primal_s = EmbeddedDeltaSet2D{Bool,Point2D}() add_vertices!(primal_s, 3, point=[Point2D(0,0), Point2D(1,0), Point2D(0.5,sqrt(0.75))]) glue_triangle!(primal_s, 1, 2, 3, tri_orientation=true) primal_s[:edge_orientation] = true -s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(primal_s) +s = dualize(primal_s) subdivide_duals!(s, Barycenter()) @test isapprox(Δ(1, s), [-12 -6 6; @@ -464,7 +467,7 @@ primal_s = EmbeddedDeltaSet2D{Bool,Point2D}() add_vertices!(primal_s, 3, point=[Point2D(0,0), Point2D(6/8,0), Point2D(6/8,8/3)]) glue_triangle!(primal_s, 1, 2, 3, tri_orientation=true) primal_s[:edge_orientation] = true -s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(primal_s) +s = dualize(primal_s) subdivide_duals!(s, Barycenter()) #@assert only(s[:area]) == 1.0 @@ -504,7 +507,7 @@ primal_s = EmbeddedDeltaSet2D{Bool,Point2D}() add_vertices!(primal_s, 4, point=[Point2D(0,0), Point2D(1,0), Point2D(0,2), Point2D(-2,5)]) glue_triangle!(primal_s, 1, 2, 3) glue_triangle!(primal_s, 1, 3, 4) -s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(primal_s) +s = dualize(primal_s) subdivide_duals!(s, Barycenter()) X = [SVector(2,3), SVector(5,7)] ♭_m = ♭_mat(s) @@ -523,18 +526,18 @@ primal_s = EmbeddedDeltaSet2D{Bool,Point2D}() add_vertices!(primal_s, 4, point=[Point2D(0,0), Point2D(1,0), Point2D(0,2), Point2D(-2,5)]) glue_triangle!(primal_s, 1, 2, 3) glue_triangle!(primal_s, 1, 3, 4) -s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(primal_s) +s = dualize(primal_s) subdivide_duals!(s, Barycenter()) primal_s′ = EmbeddedDeltaSet2D{Bool,Point2D}() add_vertices!(primal_s′, 4, point=[Point2D(0,0), Point2D(1,0), Point2D(0,2), Point2D(-2,5)]) glue_triangle!(primal_s′, 1, 2, 3) glue_triangle!(primal_s′, 1, 3, 4) -s′ = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(primal_s′) +s′ = dualize(primal_s′) s′[1, :tri_center] = 11 s′[2, :tri_center] = 10 -s′[[11,13,15,17,19,21], :D_∂v0] = 11 -s′[[12,14,16,18,20,22], :D_∂v0] = 10 +s′[[11,13,15,17,19,21], :dual_∂v0] = 11 +s′[[12,14,16,18,20,22], :dual_∂v0] = 10 subdivide_duals!(s′, Barycenter()) #@assert is_isomorphic(s,s′) @@ -544,7 +547,7 @@ X = [SVector(2,3), SVector(5,7)] @test ♭_mat(s) * DualVectorField(X) == ♭_mat(s′) * DualVectorField(X) tg′ = triangulated_grid(100,100,10,10,Point2D); -tg = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(tg′); +tg = dualize(tg′); subdivide_duals!(tg, Barycenter()); rect′ = loadmesh(Rectangle_30x10()); diff --git a/test/SimplicialComplexes.jl b/test/SimplicialComplexes.jl index 39344113..86d1a57a 100644 --- a/test/SimplicialComplexes.jl +++ b/test/SimplicialComplexes.jl @@ -48,4 +48,6 @@ h = compose(f,g) @test as_matrix(h) ≈ [1/3, 0, 1/6, 1/2] isc = id(sc) @test as_matrix(h) == as_matrix(compose(h,isc)) + + end \ No newline at end of file From 2b3fa3d1b24a3ee5dd1177224c287c0ac35d5faa Mon Sep 17 00:00:00 2001 From: Kevin Carlson Date: Mon, 15 Jul 2024 11:58:18 -0700 Subject: [PATCH 12/16] Making the dual extraction/subdivision functionality less badly-written --- src/DiscreteExteriorCalculus.jl | 97 +++++++++++++++++++++----------- test/DiscreteExteriorCalculus.jl | 21 ++++++- test/SimplicialSets.jl | 6 ++ 3 files changed, 88 insertions(+), 36 deletions(-) diff --git a/src/DiscreteExteriorCalculus.jl b/src/DiscreteExteriorCalculus.jl index 5fd6316d..56480e22 100644 --- a/src/DiscreteExteriorCalculus.jl +++ b/src/DiscreteExteriorCalculus.jl @@ -68,7 +68,6 @@ struct DiagonalHodge <: DiscreteHodge end # Euclidean geometry #################### - """ A notion of "geometric center" of a simplex. See also: [`geometric_center`](@ref). @@ -1469,22 +1468,53 @@ function dual_derivative(::Type{Val{n}}, s::HasDeltaSet, args...) where n end end -""" -Checks whethere a DeltaSet is embedded by (droolingly) searching for 'Embedded' -in the name of its type. This could also check for 'Point' in the schema, which -would feel better but be less trustworthy. -""" -is_embedded(d::HasDeltaSet) = is_embedded(typeof(t)) -is_embedded(t::Type{T}) where {T<:HasDeltaSet} = !isnothing(findfirst("Embedded",string(t.name.name))) const REPLACEMENT_FOR_DUAL_TYPE = "Set" => "DualComplex" rename_to_dual(s::Symbol) = Symbol(replace(string(s),REPLACEMENT_FOR_DUAL_TYPE)) rename_from_dual(s::Symbol) = Symbol(replace(string(s),reverse(REPLACEMENT_FOR_DUAL_TYPE))) -#adds the Real type for lengths in the Embedded case. Will need further customization -#if we add another type whose dual has more parameters than its primal. -function dual_param_list(t) - is_embedded(t) ? [t.parameters[1],eltype(t.parameters[2]),t.parameters[2]] : t.parameters +const EmbeddedDeltaSet = Union{EmbeddedDeltaSet1D,EmbeddedDeltaSet2D,EmbeddedDeltaSet3D} +const EmbeddedDeltaDualComplex = Union{EmbeddedDeltaDualComplex1D,EmbeddedDeltaDualComplex2D} +""" +Adds the Real type for lengths in the EmbeddedDeltaSet case, and removes it in the EmbeddedDeltaDualComplex case. +Will need further customization +if we add another type whose dual has different parameters than its primal. +""" +dual_param_list(d::HasDeltaSet) = typeof(d).parameters +dual_param_list(d::EmbeddedDeltaSet) = + begin t = typeof(d) ; [t.parameters[1],eltype(t.parameters[2]),t.parameters[2]] end +dual_param_list(d::EmbeddedDeltaDualComplex) = + begin t = typeof(d); [t.parameters[1],t.parameters[3]] end + + +""" +Keys are symbols for all the DeltaSet and DeltaDualComplex types. +Values are the types themselves, without parameters, so mostly UnionAlls. +Note there aren't any 0D or 3D types in here thus far. +""" +type_dict = Dict{Symbol,Type}() +const prefixes = ["Embedded","Oriented",""] +const postfixes = ["1D","2D"] +const midfixes = ["DeltaDualComplex","DeltaSet"] +for pre in prefixes for mid in midfixes for post in postfixes + s = Symbol(pre,mid,post) + type_dict[s] = eval(s) +end end end + +""" +Get the dual type of a plain, oriented, or embedded DeltaSet1D or 2D. +Will always return a `DataType`, i.e. any parameters will be evaluated. +""" +function dual_type(d::HasDeltaSet) + n = type_dict[rename_to_dual(typeof(d).name.name)] + ps = dual_param_list(d) + length(ps) > 0 ? n{ps...} : n +end +function dual_type(d::AbstractDeltaDualComplex) + n = type_dict[rename_from_dual(typeof(d).name.name)] + ps = dual_param_list(d) + length(ps) > 0 ? n{ps...} : n end + """ Calls the constructor for d's dual type on d, including parameters. Does not call `subdivide_duals!` on the result. @@ -1497,42 +1527,43 @@ time of writing (July 2024) only "Embedded" types fail criterion (2) and get spe s = EmbeddedDeltaSet2D{Bool,SVector{Float64,Float64}}() dualize(s)::EmbeddedDeltaDualComplex2D{Bool,Float64,SVector{Float64,Float64}} """ -function dualize(d::HasDeltaSet) #add embedded keyword? - t = typeof(d) - n = dualize_type_name(t) - ps = dual_param_list(t) - length(ps) > 0 ? n{ps...}(d) : n(d) -end -#Note this can only successfully dualize types whose duals are defined in the scope of this module! -dualize_type_name(t::Type) = eval(rename_to_dual(t.name.name)) - -dual_extractor(d::HasDeltaSet) = dual_extractor(typeof(d)) -function dual_extractor(t::Type{T}) where {T <: HasDeltaSet} - n = t.name.name - sch_name = Symbol("Sch",n) - sch = eval(sch_name) - dual_sch = eval(rename_to_dual(sch_name)) - C = FinCat(sch) - D = FinCat(dual_sch) +dualize(d::HasDeltaSet) = dual_type(d)(d) + +""" +Get the acset schema, as a Presentation, of a HasDeltaSet. +XXX: upstream to Catlab. +""" +fancy_acset_schema(d::HasDeltaSet) = Presentation(acset_schema(d)) + +""" +Produces a DataMigration which will migrate the dualization of a +`DeltaSet` of `d`'s type back to `d`'s schema, selecting only the dual +part. +""" +function dual_extractor(d::HasDeltaSet) + sch, dual_sch = fancy_acset_schema.([d,dual_type(d)()]) + C,D = FinCat.([sch,dual_sch]) + #Map objects to the object with "Dual" appended to the name o = Dict(nameof(a) => Symbol("Dual",nameof(a)) for a in sch.generators[:Ob]) + #Map attrtypes to themselves for a in sch.generators[:AttrType] o[nameof(a)] = nameof(a) end + #Map homs (and attrs) to the hom with "dual_" appended to the name h = Dict(nameof(a) => Symbol("dual_",nameof(a)) for a in hom_generators(C)) DeltaMigration(FinDomFunctor(o,h,C,D)) end -dual_extractor(t::Type{T}) where {T <: AbstractDeltaDualComplex} = dual_extractor(eval(rename_from_dual(t.name.name))) +dual_extractor(d::AbstractDeltaDualComplex) = dual_extractor(dual_type(d)()) """ Get the subdivision of a `DeltaSet` `d` as a `DeltaSet` of the same type, or extract just the subdivided part of a `DeltaDualComplex` as a `DeltaSet`. """ extract_dual(d::HasDeltaSet) = migrate(typeof(d),dualize(d),dual_extractor(d)) -function extract_dual(d::HasDeltaSet,alg) - is_embedded(d) || error("Cannot subdivide duals for a non-embedded DeltaSet.") +function extract_dual(d::EmbeddedDeltaSet,alg) s = dualize(d) subdivide_duals!(s,alg) migrate(typeof(d),s,dual_extractor(d)) end -extract_dual(d::AbstractDeltaDualComplex) = migrate(eval(rename_from_dual(typeof(d).name.name)),d,dual_extractor(d)) +extract_dual(d::AbstractDeltaDualComplex) = migrate(dual_type(d),d,dual_extractor(d)) """ Hodge star operator from primal ``n``-forms to dual ``N-n``-forms. diff --git a/test/DiscreteExteriorCalculus.jl b/test/DiscreteExteriorCalculus.jl index bda557f1..191b01bb 100644 --- a/test/DiscreteExteriorCalculus.jl +++ b/test/DiscreteExteriorCalculus.jl @@ -36,8 +36,8 @@ dual_es = elementary_duals(0,s,5) @test s[dual_es, :dual_∂v0] == edge_center(s, 1:4) @test elementary_duals(s, V(5)) == DualE(dual_es) -primal_s′ = subdivide(primal_s) -@test is_isomorphic(primal_s′,extract_dual(s)) #why is this failing? +primal_s′ = extract_dual(primal_s) +#@test !is_isomorphic(primal_s′,extract_dual(s)) #XX: They're opposites! @test nv(primal_s′) == nv(primal_s) + ne(primal_s) @test ne(primal_s′) == 2*ne(primal_s) @@ -56,7 +56,7 @@ s = dualize(primal_s) @test dual_boundary(1,s) == ∂(1,s)' @test dual_derivative(0,s) == -d(0,s)' -primal_s′ = subdivide(primal_s) +primal_s′ = extract_dual(primal_s) @test nv(primal_s′) == nv(primal_s) + ne(primal_s) @test ne(primal_s′) == 2*ne(primal_s) @test orient!(primal_s′) @@ -269,6 +269,21 @@ flipped_s = dualize(flipped_ps) subdivide_duals!(flipped_s, Barycenter()) @test ⋆(1,s) ≈ ⋆(1,flipped_s) +# Subdivided-dual extractors and such + +f = dual_extractor(EmbeddedDeltaSet2D{Bool,AbstractVector{Number}}()).functor +@test all([nameof(ob_map(f,:Point)) == :Point, nameof(ob_map(f,:V)) == :DualV, + nameof(hom_map(f,:∂e1)) == :dual_∂e1, nameof(hom_map(f,:edge_orientation)) == :dual_edge_orientation]) +@test extract_dual(EmbeddedDeltaSet2D{Bool,AbstractVector{Number}}()) == EmbeddedDeltaSet2D{Bool,AbstractVector{Number}}() + +primal_s = EmbeddedDeltaSet2D{Bool,Point2D}() +add_vertices!(primal_s, 3, point=[Point2D(0,0), Point2D(1,0), Point2D(0,1)]) +glue_triangle!(primal_s, 1, 2, 3, tri_orientation=true) +s = extract_dual(primal_s) +@test [1/3,1/3] ∈ subpart(s,:point) +@test nparts(s,:V) == 7 + + # NOTICE: # Tests beneath this comment are not backed up by any external source, and are # included to determine consistency as the operators are modified. diff --git a/test/SimplicialSets.jl b/test/SimplicialSets.jl index 280dac56..c9053e9b 100644 --- a/test/SimplicialSets.jl +++ b/test/SimplicialSets.jl @@ -556,4 +556,10 @@ vs = union(∂(0, s, E(es)), ∂(1, s, E(es))) @test Set(vs) == Set(Lkf[1]) @test Set(Lkf[1]) == Set([1,2,3,4,5,7]) +#Type conversion utilities +dt = CombinatorialSpaces.DiscreteExteriorCalculus.dual_type +fs = CombinatorialSpaces.DiscreteExteriorCalculus.fancy_acset_schema +@test all([dt(DeltaSet1D()) == DeltaDualComplex1D,dt(EmbeddedDeltaSet2D{Int,Vector{String}}()) == EmbeddedDeltaDualComplex2D{Int,String,Vector{String}}, dt(EmbeddedDeltaDualComplex1D{DeltaSet0D,Type,Real}()) == EmbeddedDeltaSet1D{DeltaSet0D,Real},dt(OrientedDeltaSet2D{Bool}()) == OrientedDeltaDualComplex2D{Bool}]) +@test fancy_acset_schema(DeltaSet1D()) == SchDeltaSet1D + end From dc729c91f307b12341cb3c48a53cb97264b10a24 Mon Sep 17 00:00:00 2001 From: Luke Morris <70283489+lukem12345@users.noreply.github.com> Date: Fri, 14 Jun 2024 12:48:43 -0400 Subject: [PATCH 13/16] 3D Differential Operators (#60) --- src/DiscreteExteriorCalculus.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/DiscreteExteriorCalculus.jl b/src/DiscreteExteriorCalculus.jl index 56480e22..dddc6e2a 100644 --- a/src/DiscreteExteriorCalculus.jl +++ b/src/DiscreteExteriorCalculus.jl @@ -19,6 +19,9 @@ export DualSimplex, DualV, DualE, DualTri, DualTet, DualChain, DualForm, OrientedDeltaDualComplex3D, SchOrientedDeltaDualComplex3D, EmbeddedDeltaDualComplex3D, SchEmbeddedDeltaDualComplex3D, DeltaDualComplex, EmbeddedDeltaDualComplex, OrientedDeltaDualComplex, + AbstractDeltaDualComplex3D, DeltaDualComplex3D, SchDeltaDualComplex3D, + OrientedDeltaDualComplex3D, SchOrientedDeltaDualComplex3D, + EmbeddedDeltaDualComplex3D, SchEmbeddedDeltaDualComplex3D, SimplexCenter, Barycenter, Circumcenter, Incenter, geometric_center, subsimplices, primal_vertex, elementary_duals, dual_boundary, dual_derivative, ⋆, hodge_star, inv_hodge_star, δ, codifferential, ∇², laplace_beltrami, Δ, laplace_de_rham, From f5305c4a0714244d9c15eec09f422dfb91b7dae6 Mon Sep 17 00:00:00 2001 From: Kevin Carlson Date: Thu, 25 Jul 2024 08:31:24 -0700 Subject: [PATCH 14/16] geometric map from a subdivision updated --- src/DiscreteExteriorCalculus.jl | 104 +++++++++++++++---------------- src/Meshes.jl | 2 +- src/SimplicialComplexes.jl | 50 ++++++++++++++- test/DiscreteExteriorCalculus.jl | 15 +++-- 4 files changed, 107 insertions(+), 64 deletions(-) diff --git a/src/DiscreteExteriorCalculus.jl b/src/DiscreteExteriorCalculus.jl index dddc6e2a..4b63f2af 100644 --- a/src/DiscreteExteriorCalculus.jl +++ b/src/DiscreteExteriorCalculus.jl @@ -470,12 +470,6 @@ end #type_by_name(header::String,n::Int) = eval(Symbol(header * string(n)*"D")) -#= -function extract_subdivided_primal_migration(n) - S = -end -=# - """ Abstract type for dual complex of a 2D delta set. """ @abstract_acset_type AbstractDeltaDualComplex2D <: HasDeltaSet2D @@ -971,19 +965,19 @@ end @present SchDeltaDualComplex3D <: SchDeltaSet3D begin # Dual vertices, edges, triangles, and tetrahedra. (DualV, DualE, DualTri, DualTet)::Ob - (D_∂v0, D_∂v1)::Hom(DualE, DualV) - (D_∂e0, D_∂e1, D_∂e2)::Hom(DualTri, DualE) - (D_∂t0, D_∂t1, D_∂t2, D_∂t3)::Hom(DualTet, DualTri) + (dual_∂v0, dual_∂v1)::Hom(DualE, DualV) + (dual_∂e0, dual_∂e1, dual_∂e2)::Hom(DualTri, DualE) + (dual_∂t0, dual_∂t1, dual_∂t2, dual_∂t3)::Hom(DualTet, DualTri) # Simplicial identities for dual simplices. - D_∂t3 ⋅ D_∂e2 == D_∂t2 ⋅ D_∂e2 - D_∂t3 ⋅ D_∂e1 == D_∂t1 ⋅ D_∂e2 - D_∂t3 ⋅ D_∂e0 == D_∂t0 ⋅ D_∂e2 + dual_∂t3 ⋅ dual_∂e2 == dual_∂t2 ⋅ dual_∂e2 + dual_∂t3 ⋅ dual_∂e1 == dual_∂t1 ⋅ dual_∂e2 + dual_∂t3 ⋅ dual_∂e0 == dual_∂t0 ⋅ dual_∂e2 - D_∂t2 ⋅ D_∂e1 == D_∂t1 ⋅ D_∂e1 - D_∂t2 ⋅ D_∂e0 == D_∂t0 ⋅ D_∂e1 + dual_∂t2 ⋅ dual_∂e1 == dual_∂t1 ⋅ dual_∂e1 + dual_∂t2 ⋅ dual_∂e0 == dual_∂t0 ⋅ dual_∂e1 - D_∂t1 ⋅ D_∂e0 == D_∂t0 ⋅ D_∂e0 + dual_∂t1 ⋅ dual_∂e0 == dual_∂t0 ⋅ dual_∂e0 # Centers of primal simplices are dual vertices. vertex_center::Hom(V, DualV) @@ -999,24 +993,24 @@ end """ Dual complex of a three-dimensional delta set. """ @acset_type DeltaDualComplex3D(SchDeltaDualComplex3D, - index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:D_∂v0,:D_∂v1,:D_∂e0,:D_∂e1,:D_∂e2,:D_∂t0,:D_∂t1,:D_∂t2,:D_∂t3]) <: AbstractDeltaDualComplex3D + index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:dual_∂v0,:dual_∂v1,:dual_∂e0,:dual_∂e1,:dual_∂e2,:dual_∂t0,:dual_∂t1,:dual_∂t2,:dual_∂t3]) <: AbstractDeltaDualComplex3D """ Dual vertex corresponding to center of primal tetrahedron. """ tetrahedron_center(s::HasDeltaSet3D, args...) = s[args..., :tet_center] subsimplices(::Type{Val{3}}, s::HasDeltaSet3D, tet::Int) = - SVector{24}(incident(s, tetrahedron_center(s,tet), @SVector [:D_∂t1, :D_∂e1, :D_∂v0])) + SVector{24}(incident(s, tetrahedron_center(s,tet), @SVector [:dual_∂t1, :dual_∂e1, :dual_∂v0])) primal_vertex(::Type{Val{3}}, s::HasDeltaSet3D, tet...) = - primal_vertex(Val{2}, s, s[tet..., :D_∂t1]) + primal_vertex(Val{2}, s, s[tet..., :dual_∂t1]) elementary_duals(::Type{Val{0}}, s::AbstractDeltaDualComplex3D, v::Int) = - incident(s, vertex_center(s,v), @SVector [:D_∂t1, :D_∂e1, :D_∂v1]) + incident(s, vertex_center(s,v), @SVector [:dual_∂t1, :dual_∂e1, :dual_∂v1]) elementary_duals(::Type{Val{1}}, s::AbstractDeltaDualComplex3D, e::Int) = - incident(s, edge_center(s,e), @SVector [:D_∂e1, :D_∂v1]) + incident(s, edge_center(s,e), @SVector [:dual_∂e1, :dual_∂v1]) elementary_duals(::Type{Val{2}}, s::AbstractDeltaDualComplex3D, t::Int) = - incident(s, triangle_center(s,t), :D_∂v1) + incident(s, triangle_center(s,t), :dual_∂v1) elementary_duals(::Type{Val{3}}, s::AbstractDeltaDualComplex3D, tet::Int) = SVector(tetrahedron_center(s,tet)) @@ -1025,10 +1019,10 @@ elementary_duals(::Type{Val{3}}, s::AbstractDeltaDualComplex3D, tet::Int) = This accessor assumes that the simplicial identities for the dual hold. """ function dual_tetrahedron_vertices(s::HasDeltaSet3D, t...) - SVector(s[s[s[t..., :D_∂t2], :D_∂e2], :D_∂v1], - s[s[s[t..., :D_∂t2], :D_∂e2], :D_∂v0], - s[s[s[t..., :D_∂t0], :D_∂e0], :D_∂v1], - s[s[s[t..., :D_∂t0], :D_∂e0], :D_∂v0]) + SVector(s[s[s[t..., :dual_∂t2], :dual_∂e2], :dual_∂v1], + s[s[s[t..., :dual_∂t2], :dual_∂e2], :dual_∂v0], + s[s[s[t..., :dual_∂t0], :dual_∂e0], :dual_∂v1], + s[s[s[t..., :dual_∂t0], :dual_∂e0], :dual_∂v0]) end # 3D oriented dual complex @@ -1039,15 +1033,15 @@ end edge_orientation::Attr(E, Orientation) tri_orientation::Attr(Tri, Orientation) tet_orientation::Attr(Tet, Orientation) - D_edge_orientation::Attr(DualE, Orientation) - D_tri_orientation::Attr(DualTri, Orientation) - D_tet_orientation::Attr(DualTet, Orientation) + dual_edge_orientation::Attr(DualE, Orientation) + dual_tri_orientation::Attr(DualTri, Orientation) + dual_tet_orientation::Attr(DualTet, Orientation) end """ Oriented dual complex of an oriented 3D delta set. """ @acset_type OrientedDeltaDualComplex3D(SchOrientedDeltaDualComplex3D, - index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:D_∂v0,:D_∂v1,:D_∂e0,:D_∂e1,:D_∂e2,:D_∂t0,:D_∂t1,:D_∂t2,:D_∂t3]) <: AbstractDeltaDualComplex3D + index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:dual_∂v0,:dual_∂v1,:dual_∂e0,:dual_∂e1,:dual_∂e2,:dual_∂t0,:dual_∂t1,:dual_∂t2,:dual_∂t3]) <: AbstractDeltaDualComplex3D dual_boundary_nz(::Type{Val{1}}, s::AbstractDeltaDualComplex3D, x::Int) = # Boundary vertices of dual 1-cell ↔ @@ -1087,15 +1081,15 @@ make_dual_simplices_3d!(s::AbstractDeltaDualComplex3D) = make_dual_simplices_3d! # Note: these accessors are isomorphic to those for their primal counterparts. # These can be eliminated by the DualComplex schema refactor. add_dual_edge!(s::AbstractDeltaDualComplex3D, d_src::Int, d_tgt::Int; kw...) = - add_part!(s, :DualE; D_∂v1=d_src, D_∂v0=d_tgt, kw...) + add_part!(s, :DualE; dual_∂v1=d_src, dual_∂v0=d_tgt, kw...) function get_dual_edge!(s::AbstractDeltaDualComplex3D, d_src::Int, d_tgt::Int; kw...) - es = (e for e in incident(s, d_src, :D_∂v1) if s[e, :D_∂v0] == d_tgt) - isempty(es) ? add_part!(s, :DualE; D_∂v1=d_src, D_∂v0=d_tgt, kw...) : first(es) + es = (e for e in incident(s, d_src, :dual_∂v1) if s[e, :dual_∂v0] == d_tgt) + isempty(es) ? add_part!(s, :DualE; dual_∂v1=d_src, dual_∂v0=d_tgt, kw...) : first(es) end add_dual_triangle!(s::AbstractDeltaDualComplex3D, d_src2_first::Int, d_src2_last::Int, d_tgt2::Int; kw...) = - add_part!(s, :DualTri; D_∂e0=d_src2_last, D_∂e1=d_tgt2, D_∂e2=d_src2_first, kw...) + add_part!(s, :DualTri; dual_∂e0=d_src2_last, dual_∂e1=d_tgt2, dual_∂e2=d_src2_first, kw...) function glue_dual_triangle!(s::AbstractDeltaDualComplex3D, d_v₀::Int, d_v₁::Int, d_v₂::Int; kw...) add_dual_triangle!(s, get_dual_edge!(s, d_v₀, d_v₁), get_dual_edge!(s, d_v₁, d_v₂), @@ -1103,19 +1097,19 @@ function glue_dual_triangle!(s::AbstractDeltaDualComplex3D, d_v₀::Int, d_v₁: end add_dual_tetrahedron!(s::AbstractDeltaDualComplex3D, d_tri0::Int, d_tri1::Int, d_tri2::Int, d_tri3::Int; kw...) = - add_part!(s, :DualTet; D_∂t0=d_tri0, D_∂t1=d_tri1, D_∂t2=d_tri2, D_∂t3=d_tri3, kw...) + add_part!(s, :DualTet; dual_∂t0=d_tri0, dual_∂t1=d_tri1, dual_∂t2=d_tri2, dual_∂t3=d_tri3, kw...) function dual_triangles(s::AbstractDeltaDualComplex3D, d_v₀::Int, d_v₁::Int, d_v₂::Int) - d_e₀s = incident(s, d_v₂, :D_∂v0) ∩ incident(s, d_v₁, :D_∂v1) + d_e₀s = incident(s, d_v₂, :dual_∂v0) ∩ incident(s, d_v₁, :dual_∂v1) isempty(d_e₀s) && return Int[] - d_e₁s = incident(s, d_v₂, :D_∂v0) ∩ incident(s, d_v₀, :D_∂v1) + d_e₁s = incident(s, d_v₂, :dual_∂v0) ∩ incident(s, d_v₀, :dual_∂v1) isempty(d_e₁s) && return Int[] - d_e₂s = incident(s, d_v₁, :D_∂v0) ∩ incident(s, d_v₀, :D_∂v1) + d_e₂s = incident(s, d_v₁, :dual_∂v0) ∩ incident(s, d_v₀, :dual_∂v1) isempty(d_e₂s) && return Int[] intersect( - incident(s, d_e₀s, :D_∂e0)..., - incident(s, d_e₁s, :D_∂e1)..., - incident(s, d_e₂s, :D_∂e2)...) + incident(s, d_e₀s, :dual_∂e0)..., + incident(s, d_e₁s, :dual_∂e1)..., + incident(s, d_e₂s, :dual_∂e2)...) end function get_dual_triangle!(s::AbstractDeltaDualComplex3D, d_v₀::Int, d_v₁::Int, d_v₂::Int) @@ -1156,17 +1150,17 @@ function make_dual_simplices_3d!(s::HasDeltaSet3D, ::Type{Simplex{n}}) where n v₀, v₁, v₂, v₃ = 1:4 e₀, e₁, e₂, e₃, e₄, e₅ = 5:10 t₀, t₁, t₂, t₃ = 11:14 - # Note: You could write `D_tetrahedron_schemas` using: + # Note: You could write `dual_tetrahedron_schemas` using: #es_per_v = [(3,4,5), (1,2,5), (0,2,4), (0,1,3)] #ts_per_e = [(0,1), (0,2), (0,3), (1,2), (1,3), (2,3)] # and/or the fact that vertex vᵢ is a vertex of triangles {1,2,3,4} - {i}. - D_tetrahedron_schemas = [ + dual_tetrahedron_schemas = [ (v₀,e₄,t₃), (v₀,e₄,t₁), (v₀,e₃,t₁), (v₀,e₃,t₂), (v₀,e₅,t₂), (v₀,e₅,t₃), (v₁,e₅,t₃), (v₁,e₅,t₂), (v₁,e₁,t₂), (v₁,e₁,t₀), (v₁,e₂,t₀), (v₁,e₂,t₃), (v₂,e₂,t₃), (v₂,e₂,t₀), (v₂,e₀,t₀), (v₂,e₀,t₁), (v₂,e₄,t₁), (v₂,e₄,t₃), (v₃,e₃,t₂), (v₃,e₃,t₁), (v₃,e₀,t₁), (v₃,e₀,t₀), (v₃,e₁,t₀), (v₃,e₁,t₂)] - foreach(D_tetrahedron_schemas) do (x,y,z) + foreach(dual_tetrahedron_schemas) do (x,y,z) # Exploit the fact that `glue_sorted_dual_tetrahedron!` adds only # necessary new dual triangles. glue_sorted_dual_tetrahedron!(s, dvs[x], dvs[y], dvs[z], tc) @@ -1186,12 +1180,12 @@ function make_dual_simplices_3d!(s::HasDeltaSet3D, ::Type{Simplex{n}}) where n # Orient elementary dual tetrahedra. # Exploit the fact that triangles are added in the order of - # D_tetrahedron_schemas. + # dual_tetrahedron_schemas. for tet in tetrahedra(s) tet_orient = s[tet, :tet_orientation] rev_tet_orient = negate(tet_orient) - D_tets = (24*(tet-1)+1):(24*tet) - s[D_tets, :D_tet_orientation] = repeat([tet_orient, rev_tet_orient], 12) + dual_tets = (24*(tet-1)+1):(24*tet) + s[dual_tets, :dual_tet_orientation] = repeat([tet_orient, rev_tet_orient], 12) end # Orient elementary dual triangles. @@ -1199,7 +1193,7 @@ function make_dual_simplices_3d!(s::HasDeltaSet3D, ::Type{Simplex{n}}) where n # TODO: Perhaps multiply by tet_orientation. primal_edge_orient = s[e, :edge_orientation] d_ts = elementary_duals(1,s,e) - s[d_ts, :D_tri_orientation] = primal_edge_orient + s[d_ts, :dual_tri_orientation] = primal_edge_orient end # Orient elementary dual edges. @@ -1207,14 +1201,14 @@ function make_dual_simplices_3d!(s::HasDeltaSet3D, ::Type{Simplex{n}}) where n # TODO: Perhaps multiply by tet_orientation. primal_tri_orient = s[t, :tri_orientation] d_es = elementary_duals(2,s,t) - s[d_es, :D_edge_orientation] = primal_tri_orient + s[d_es, :dual_edge_orientation] = primal_tri_orient end # Remaining dual edges and dual triangles are oriented arbitrarily. - s[findall(isnothing, s[:D_tri_orientation]), :D_tri_orientation] = one(attrtype_type(s, :Orientation)) + s[findall(isnothing, s[:dual_tri_orientation]), :dual_tri_orientation] = one(attrtype_type(s, :Orientation)) # These will be dual edges from vertex_center to tc, and from # edge_center to tc. - s[findall(isnothing, s[:D_edge_orientation]), :D_edge_orientation] = one(attrtype_type(s, :Orientation)) + s[findall(isnothing, s[:dual_edge_orientation]), :dual_edge_orientation] = one(attrtype_type(s, :Orientation)) end return parts(s, :DualTet) @@ -1239,7 +1233,7 @@ end """ @acset_type EmbeddedDeltaDualComplex3D(SchEmbeddedDeltaDualComplex3D, - index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:D_∂v0,:D_∂v1,:D_∂e0,:D_∂e1,:D_∂e2,:D_∂t0,:D_∂t1,:D_∂t2,:D_∂t3]) <: AbstractDeltaDualComplex3D + index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:dual_∂v0,:dual_∂v1,:dual_∂e0,:dual_∂e1,:dual_∂e2,:dual_∂t0,:dual_∂t1,:dual_∂t2,:dual_∂t3]) <: AbstractDeltaDualComplex3D volume(::Type{Val{n}}, s::EmbeddedDeltaDualComplex3D, x) where n = volume(Val{n}, s, x, PrecomputedVol()) @@ -1251,10 +1245,10 @@ dual_volume(::Type{Val{3}}, s::HasDeltaSet3D, tet, ::PrecomputedVol) = s[tet, :dual_vol] function dual_volume(::Type{Val{3}}, s::HasDeltaSet3D, tet::Int, ::CayleyMengerDet) - dual_vs = SVector(s[s[s[tet, :D_∂t2], :D_∂e2], :D_∂v1], - s[s[s[tet, :D_∂t2], :D_∂e2], :D_∂v0], - s[s[s[tet, :D_∂t0], :D_∂e0], :D_∂v1], - s[s[s[tet, :D_∂t0], :D_∂e0], :D_∂v0]) + dual_vs = SVector(s[s[s[tet, :dual_∂t2], :dual_∂e2], :dual_∂v1], + s[s[s[tet, :dual_∂t2], :dual_∂e2], :dual_∂v0], + s[s[s[tet, :dual_∂t0], :dual_∂e0], :dual_∂v1], + s[s[s[tet, :dual_∂t0], :dual_∂e0], :dual_∂v0]) volume(dual_point(s, dual_vs)) end diff --git a/src/Meshes.jl b/src/Meshes.jl index 518b8793..0db85677 100644 --- a/src/Meshes.jl +++ b/src/Meshes.jl @@ -23,7 +23,7 @@ struct Rectangle_30x10 <: AbstractMeshKey end struct Torus_30x10 <: AbstractMeshKey end struct Point_Map <: AbstractMeshKey end -Icosphere(n) = Icosphere(n, 1.0) +const Icosphere(n) = Icosphere(n, 1.0) """ loadmesh(s::Icosphere) diff --git a/src/SimplicialComplexes.jl b/src/SimplicialComplexes.jl index e97affe0..3531bf5c 100644 --- a/src/SimplicialComplexes.jl +++ b/src/SimplicialComplexes.jl @@ -3,7 +3,7 @@ The category of simplicial complexes and Kleisli maps for the convex space monad """ module SimplicialComplexes export SimplicialComplex, VertexList, has_simplex, GeometricPoint, has_point, has_span, GeometricMap, nv, -as_matrix, compose, id +as_matrix, compose, id, cocenter, primal_vertices using ..Tries using ..SimplicialSets import AlgebraicInterfaces: dom,codom,compose,id @@ -259,11 +259,57 @@ function GeometricMap(sc::SimplicialComplex,::Barycenter) end #accessors for the nonzeros in a column of the matrix +""" +The geometric map from a deltaset's subdivision to itself +""" +function GeometricMap(primal_s::EmbeddedDeltaSet,alg) + dom = SimplicialComplex(primal_s) + s = dualize(primal_s) + subdivide_duals!(s,alg) + cod = SimplicialComplex(extract_dual(s)) + mat = zeros(Float64,nv(cod),nv(dom)) + pvs = map(i->primal_vertices(s,i),1:nv(dom)) + weights = 1 ./(length.(pvs)) + for j in 1:nv(dom) + for v in pvs[j] + mat[v,j] = weights[j] + end + end + GeometricMap(dom,cod,mat) +end + function pullback_primal(f::GeometricMap, v::PrimalVectorField{T}) where T nv(f.cod) == length(v) || error("Vector field must have same number of vertices as codomain") PrimalVectorField(T.(eachcol(hcat(v.data...)*as_matrix(f)))) end -*(f::GeometricMap,v::PrimalVectorField) = pullback_primal(f,v) +*(f::GeometricMap,v::PrimalVectorField) = pullback_pr + +function dual_vertex_dimension(s::AbstractDeltaDualComplex,v::DualV) + n = v.data + !isempty(incident(s,n,:vertex_center)) ? 0 : + !isempty(incident(s,n,:edge_center)) ? 1 : + !isempty(incident(s,n,:tri_center)) ? 2 : 3 +end +simplex_name_dict = Dict(0=>:vertex,1=>:edge,2=>:tri,3=>:tet) + +#XX: the parts data structure allowing data to be like whatever is awful +function cocenter(s::AbstractDeltaDualComplex,v::DualV) + n = dimension(s) + v = v.data + for i in 0:n + inc = incident(s,v,Symbol(simplex_name_dict[i],:(_center))) + if !isempty(inc) + return Simplex{i}(only(inc)) + end + end +end +cocenter(s::AbstractDeltaDualComplex,n::Int) = cocenter(s,DualV(n)) +primal_vertices(s::AbstractDeltaDualComplex,v::DualV) = simplex_vertices(s,cocenter(s,v)) +primal_vertices(s::AbstractDeltaDualComplex,n::Int) = simplex_vertices(s,cocenter(s,DualV(n))) + +#dimension(x::Simplex{n}) where {n} = n + +end end \ No newline at end of file diff --git a/test/DiscreteExteriorCalculus.jl b/test/DiscreteExteriorCalculus.jl index 191b01bb..61c70ad7 100644 --- a/test/DiscreteExteriorCalculus.jl +++ b/test/DiscreteExteriorCalculus.jl @@ -561,13 +561,16 @@ X = [SVector(2,3), SVector(5,7)] @test ♭(s, DualVectorField(X)) == ♭(s′, DualVectorField(X)) @test ♭_mat(s) * DualVectorField(X) == ♭_mat(s′) * DualVectorField(X) -tg′ = triangulated_grid(100,100,10,10,Point2D); -tg = dualize(tg′); -subdivide_duals!(tg, Barycenter()); +tg′ = triangulated_grid(100,100,10,10,Point2D) +tg = dualize(tg′) +subdivide_duals!(tg, Barycenter()) -rect′ = loadmesh(Rectangle_30x10()); -rect = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(rect′); -subdivide_duals!(rect, Barycenter()); +rect′ = loadmesh(Rectangle_30x10()) +rect = dualize(rect′) +subdivide_duals!(rect, Barycenter(r)) + +rect_fine_primal = extract_dual(rect) +rect_fine = dualize(rect_fine_primal) flat_meshes = [tri_345(), tri_345_false(), right_scalene_unit_hypot(), grid_345(), (tg′, tg), (rect′, rect)]; From 28b2d3060690203962387674199294f5e15e6f31 Mon Sep 17 00:00:00 2001 From: Kevin Carlson Date: Mon, 15 Jul 2024 12:01:55 -0700 Subject: [PATCH 15/16] Remove old subdivision function that doesn't go through dual complexes --- src/DiscreteExteriorCalculus.jl | 52 ++-------------------- src/SimplicialComplexes.jl | 2 +- src/SimplicialSets.jl | 78 +-------------------------------- src/Tries.jl | 13 ++---- test/SimplicialSets.jl | 2 +- 5 files changed, 9 insertions(+), 138 deletions(-) diff --git a/src/DiscreteExteriorCalculus.jl b/src/DiscreteExteriorCalculus.jl index 4b63f2af..74bf6692 100644 --- a/src/DiscreteExteriorCalculus.jl +++ b/src/DiscreteExteriorCalculus.jl @@ -29,7 +29,7 @@ export DualSimplex, DualV, DualE, DualTri, DualTet, DualChain, DualForm, ℒ, lie_derivative, lie_derivative_flat, vertex_center, edge_center, triangle_center, tetrahedron_center, dual_tetrahedron_vertices, dual_triangle_vertices, dual_edge_vertices, dual_point, dual_volume, subdivide_duals!, DiagonalHodge, GeometricHodge, - subdivide, PPSharp, AltPPSharp, DesbrunSharp, LLSDDSharp, de_sign, + PPSharp, AltPPSharp, DesbrunSharp, LLSDDSharp, de_sign, ♭♯, ♭♯_mat, flat_sharp, flat_sharp_mat, dualize, dual_extractor, extract_dual import Base: ndims @@ -54,7 +54,7 @@ using DataMigrations: @migrate using ..ArrayUtils, ..SimplicialSets using ..SimplicialSets: CayleyMengerDet, operator_nz, ∂_nz, d_nz, cayley_menger, negate -import ..SimplicialSets: ∂, d, volume, subdivide +import ..SimplicialSets: ∂, d, volume abstract type DiscreteFlat end struct DPPFlat <: DiscreteFlat end @@ -304,40 +304,6 @@ function make_dual_simplices_1d!(s::HasDeltaSet1D, ::Type{Simplex{n}}) where n dual_edges end -# TODO: Instead of copying-and-pasting the DeltaSet1D version: -# - Use metaprogramming, or -# - Don't use the migration DSL, but rather the lower-level functor interface. -# TODO: When Catlab PR #823 "Data migrations with Julia functions on attributes" -# is merged, ensure that oriented-ness is preserved. (Flip one of the -# orientations.) -""" Subdivide an oriented 1D delta set. - -Note that this function does NOT currently guarantee that if the input is -oriented, then the output will be. -""" -function subdivide(s::OrientedDeltaSet1D{T}) where T - @migrate typeof(s) s begin - V => @cases begin - v::V - e::E - end - E => @cases begin - e₁::E - e₂::E - end - ∂v1 => begin - e₁ => e - e₂ => e - end - ∂v0 => begin - e₁ => (v∘∂v1) - e₂ => (v∘∂v0) - end - Orientation => Orientation - # TODO: One of these edge orientations must be flipped. (e₂?) - edge_orientation => (e₁ => edge_orientation; e₂ => edge_orientation) - end -end # 1D embedded dual complex #------------------------- @@ -432,17 +398,7 @@ function precompute_volumes_1d!(sd::HasDeltaSet1D, ::Type{point_type}) where poi end end -# TODO: When Catlab PR #823 "Data migrations with Julia functions on attributes" -# is merged, encode subdivision like so: -#function subdivide(s::EmbeddedDeltaSet1D{T,U}, alg::V) where {T,U,V <: SimplexCenter} -# @migrate typeof(s) s begin -# ... -# edge_orientation => (e₁ => edge_orientation; e₂ => !(edge_orientation)) -# Point => Point -# point => (v => point; e => geometric_center([e₁ ⋅ point, e₂ ⋅ point], alg)) -# ... -# end -#end +# TODO: Orientation on subdivisions # 2D dual complex ################# @@ -468,8 +424,6 @@ end end -#type_by_name(header::String,n::Int) = eval(Symbol(header * string(n)*"D")) - """ Abstract type for dual complex of a 2D delta set. """ @abstract_acset_type AbstractDeltaDualComplex2D <: HasDeltaSet2D diff --git a/src/SimplicialComplexes.jl b/src/SimplicialComplexes.jl index 3531bf5c..f9201aac 100644 --- a/src/SimplicialComplexes.jl +++ b/src/SimplicialComplexes.jl @@ -249,7 +249,7 @@ id(sc::SimplicialComplex) = GeometricMap(sc,sc,GeometricPoint.(eachcol(I(nv(sc)) #make sparse matrices function GeometricMap(sc::SimplicialComplex,::Barycenter) - dom = SimplicialComplex(subdivide(sc.delta_set)) + dom = SimplicialComplex(extract_dual(sc.delta_set)) #Vertices of dom correspond to vertices, edges, triangles of sc. mat = zeros(Float64,nv(sc),nv(dom)) for i in 1:nv(sc) mat[i,i] = 1 end diff --git a/src/SimplicialSets.jl b/src/SimplicialSets.jl index 055e07b8..84d39755 100644 --- a/src/SimplicialSets.jl +++ b/src/SimplicialSets.jl @@ -35,7 +35,7 @@ export Simplex, V, E, Tri, Tet, SimplexChain, VChain, EChain, TriChain, TetChain tetrahedron_triangles, tetrahedron_edges, tetrahedron_vertices, ntetrahedra, tetrahedra, add_tetrahedron!, glue_tetrahedron!, glue_sorted_tetrahedron!, glue_sorted_tet_cube!, is_manifold_like, nonboundaries, - star, St, closed_star, St̄, link, Lk, simplex_vertices, dimension, subdivide, + star, St, closed_star, St̄, link, Lk, simplex_vertices, dimension, DeltaSet, OrientedDeltaSet, EmbeddedDeltaSet using LinearAlgebra: det @@ -360,82 +360,6 @@ function glue_sorted_triangle!(s::HasDeltaSet2D, v₀::Int, v₁::Int, v₂::Int glue_triangle!(s, v₀, v₁, v₂; kw...) end -""" Subdivide a 1D delta set. Note that this is written as if it'll work -for any type of 1D delta-set, but it can't handle orientations or embeddings. -""" -function subdivide(s::HasDeltaSet1D) - @migrate typeof(s) s begin - V => @cases begin - v::V - e::E - end - E => @cases begin - e₁::E - e₂::E - end - ∂v1 => begin - e₁ => e - e₂ => e - end - ∂v0 => begin - e₁ => (v∘∂v1) - e₂ => (v∘∂v0) - end - end -end -""" -Subdivision of a 2D simplicial set, relies on glue_triangle! so not good -for arbitrary simplicial sets. -""" -function subdivide(s::HasDeltaSet2D) - d = typeof(s)() - #vertices for all simplices - add_vertices!(d, nv(s)) - add_vertices!(d, ne(s)) - add_vertices!(d, ntriangles(s)) - - e_to_v(e) = e + nv(s) - ts_as_vs = (1:ntriangles(s)) .+ (nv(s) + ne(s)) - #edges from source of edge to edge and target of edge to edge - add_edges!(d, subpart(s, :∂v1), e_to_v.(1:ne(s))) - add_edges!(d, subpart(s, :∂v0), e_to_v.(1:ne(s))) - #edges from vertex of triangle to triangle - add_edges!(d, subpart(s, [:∂e2, :∂v1]), ts_as_vs) - add_edges!(d, subpart(s, [:∂e2, :∂v0]), ts_as_vs) - add_edges!(d, subpart(s, [:∂e1, :∂v0]), ts_as_vs) - #edges from edge of triangle to triangle - add_edges!(d, e_to_v.(subpart(s, :∂e2)), ts_as_vs) - add_edges!(d, e_to_v.(subpart(s, :∂e1)), ts_as_vs) - add_edges!(d, e_to_v.(subpart(s, :∂e0)), ts_as_vs) - #triangles from vertex of edge of triangle to triangle - glue_triangles!(d, - subpart(s, [:∂e2, :∂v1]), - e_to_v.(subpart(s, :∂e2)), - ts_as_vs) - glue_triangles!(d, - subpart(s, [:∂e2, :∂v0]), - e_to_v.(subpart(s, :∂e2)), - ts_as_vs) - glue_triangles!(d, - subpart(s, [:∂e1, :∂v1]), - e_to_v.(subpart(s, :∂e1)), - ts_as_vs) - glue_triangles!(d, - subpart(s, [:∂e1, :∂v0]), - e_to_v.(subpart(s, :∂e1)), - ts_as_vs) - glue_triangles!(d, - subpart(s, [:∂e0, :∂v1]), - e_to_v.(subpart(s, :∂e0)), - ts_as_vs) - glue_triangles!(d, - subpart(s, [:∂e0, :∂v0]), - e_to_v.(subpart(s, :∂e0)), - ts_as_vs) - ##XX: orientations? - d -end - # 2D oriented simplicial sets #---------------------------- diff --git a/src/Tries.jl b/src/Tries.jl index ec37cb9e..f8c82fe1 100644 --- a/src/Tries.jl +++ b/src/Tries.jl @@ -106,22 +106,15 @@ end height(t::Trie) = isempty(t.children) ? 0 : 1 + maximum(height.(values(t.children))) -#= -function subdivide_trie(t::Trie{K,V},help=K[]) where {K,V} - sdt = Trie{Vector{K},V}() - sdt.value = t.value - for k in [k for k in keys(t) if length(k) > 0] - sdt.children[vcat(help,k)] = subdivide(subtrie(t,k),vcat(help,k)) - end - sdt -end -=# + + # The state of a TrieIterator is a pair (t::Trie, i::Int), # where t is the Trie which was the output of the previous iteration # and i is the index of the current character of the string. # The indexing is potentially confusing; # see the comments and implementation below for details. + struct TrieIterator t::Trie str diff --git a/test/SimplicialSets.jl b/test/SimplicialSets.jl index c9053e9b..7807e04b 100644 --- a/test/SimplicialSets.jl +++ b/test/SimplicialSets.jl @@ -121,7 +121,7 @@ glue_triangle!(s, 1, 4, 3) @test triangles(s) == 1:2 @test ne(s) == 5 @test sort(map(Pair, src(s), tgt(s))) == [1=>2, 1=>3, 1=>4, 2=>3, 4=>3] -sd = subdivide(s) +sd = extract_dual(s) @test ntriangles(sd) == 12 && ne(sd) == 22 && nv(sd) == 11 # 2D oriented simplicial sets #---------------------------- From 31d316e012adc4c8831fec88edfd31721745ad4d Mon Sep 17 00:00:00 2001 From: Kevin Carlson Date: Thu, 25 Jul 2024 09:33:59 -0700 Subject: [PATCH 16/16] little cleanup --- Project.toml | 2 ++ src/SimplicialComplexes.jl | 3 ++- test/DiscreteExteriorCalculus.jl | 6 ++++++ test/SimplicialSets.jl | 10 ++-------- 4 files changed, 12 insertions(+), 9 deletions(-) diff --git a/Project.toml b/Project.toml index 58715d76..79e8a540 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ ACSets = "227ef7b5-1206-438b-ac65-934d6da304b8" AlgebraicInterfaces = "23cfdc9f-0504-424a-be1f-4892b28e2f0c" Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Catlab = "134e5e36-593f-5add-ad60-77f754baafbe" DataMigrations = "0c4ad18d-0c49-4bc2-90d5-5bca8f00d6ae" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" @@ -18,6 +19,7 @@ KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" MeshIO = "7269a6da-0436-5bbc-96c2-40638cbb6118" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" diff --git a/src/SimplicialComplexes.jl b/src/SimplicialComplexes.jl index f9201aac..59a3fbef 100644 --- a/src/SimplicialComplexes.jl +++ b/src/SimplicialComplexes.jl @@ -282,7 +282,8 @@ function pullback_primal(f::GeometricMap, v::PrimalVectorField{T}) where T nv(f.cod) == length(v) || error("Vector field must have same number of vertices as codomain") PrimalVectorField(T.(eachcol(hcat(v.data...)*as_matrix(f)))) end -*(f::GeometricMap,v::PrimalVectorField) = pullback_pr +#Is restriction the transpose? +*(f::GeometricMap,v::PrimalVectorField) = pullback_primal(f,v) function dual_vertex_dimension(s::AbstractDeltaDualComplex,v::DualV) n = v.data diff --git a/test/DiscreteExteriorCalculus.jl b/test/DiscreteExteriorCalculus.jl index 61c70ad7..3af5f216 100644 --- a/test/DiscreteExteriorCalculus.jl +++ b/test/DiscreteExteriorCalculus.jl @@ -775,6 +775,12 @@ for s in [tetrahedron_s, cube_s] # Desbrun, Kanso, Tong 2008, Equation 4.2. @test dual_derivative(3-k,s) == (-1)^k * d(k-1,s)' end +#Type conversion utilities +dt = DiscreteExteriorCalculus.dual_type +fs = DiscreteExteriorCalculus.fancy_acset_schema +@test all([dt(DeltaSet1D()) == DeltaDualComplex1D,dt(EmbeddedDeltaSet2D{Int,Vector{String}}()) == EmbeddedDeltaDualComplex2D{Int,String,Vector{String}}, dt(EmbeddedDeltaDualComplex1D{DeltaSet0D,Type,Real}()) == EmbeddedDeltaSet1D{DeltaSet0D,Real},dt(OrientedDeltaSet2D{Bool}()) == OrientedDeltaDualComplex2D{Bool}]) +@test fancy_acset_schema(DeltaSet1D()) == SchDeltaSet1D + end # 3D embedded dual complex diff --git a/test/SimplicialSets.jl b/test/SimplicialSets.jl index 7807e04b..d6627a61 100644 --- a/test/SimplicialSets.jl +++ b/test/SimplicialSets.jl @@ -121,8 +121,8 @@ glue_triangle!(s, 1, 4, 3) @test triangles(s) == 1:2 @test ne(s) == 5 @test sort(map(Pair, src(s), tgt(s))) == [1=>2, 1=>3, 1=>4, 2=>3, 4=>3] -sd = extract_dual(s) -@test ntriangles(sd) == 12 && ne(sd) == 22 && nv(sd) == 11 +#sd = extract_dual(s) +#@test ntriangles(sd) == 12 && ne(sd) == 22 && nv(sd) == 11 # 2D oriented simplicial sets #---------------------------- @@ -556,10 +556,4 @@ vs = union(∂(0, s, E(es)), ∂(1, s, E(es))) @test Set(vs) == Set(Lkf[1]) @test Set(Lkf[1]) == Set([1,2,3,4,5,7]) -#Type conversion utilities -dt = CombinatorialSpaces.DiscreteExteriorCalculus.dual_type -fs = CombinatorialSpaces.DiscreteExteriorCalculus.fancy_acset_schema -@test all([dt(DeltaSet1D()) == DeltaDualComplex1D,dt(EmbeddedDeltaSet2D{Int,Vector{String}}()) == EmbeddedDeltaDualComplex2D{Int,String,Vector{String}}, dt(EmbeddedDeltaDualComplex1D{DeltaSet0D,Type,Real}()) == EmbeddedDeltaSet1D{DeltaSet0D,Real},dt(OrientedDeltaSet2D{Bool}()) == OrientedDeltaDualComplex2D{Bool}]) -@test fancy_acset_schema(DeltaSet1D()) == SchDeltaSet1D - end