Setup
I am trying to write a function that determines if a set of vertices, edges and faces is a pure simplicial complex.
A pure simplicial complex is a set where all facets have the same degree, a facet is a simplex that is not contained in a simplex of higher degree. So for example, a single vertex is a facet. A set of vertices is a pure simplicial complex because all the facets are of degree one and so on.
In this case a simplicial complex is up to order 2 (triangles) so it is expressed as 2 matrices.
An $|edges|\times|vertices|$ adjacency matrix and a $|face|\times|edges|$ adjacency matrix. If you are curious about the details here's a cool resource on simplicial complexes in the context of CS: https://www.cs.cmu.edu/~kmcrane/Projects/DDG/
I represent a subset of the complex as 3 vectors, a vertex vector, an edge vector and a face vector, where every entry is 0 or 1. 1 means the element is in the subset, 0 means it isn't.
Attempted solution
So to determine if a subset is a pure simplicial complex I thought of the following approach. let's say we want to find all the vertex facets first.
0 out all the edge rows in the edge adjacency matrix that are not in the subset. This is achieved by the product $\text{sub_edge_adjacency} = D(edges) \times \text{edge_adjacency}$ where $D$ is an operator that generates a diagonal matrix from a vector. Then multiply the above by a vector filled with 1's of the same dimension as the number of edges and multiply from the right $ones * \text{subedge_adjacency}$. Which is equivalent to just adding all the rows into a single vector...
From then on...
I am not going to keep explaining the algorithm as it should be obvious this is convoluted. I suspect there is a much nicer way to check if a subset is a pure simplicial complex. Essentially I am trying to code the 8th programming problem of chapter 2 in the resource I linked.
I tested my monstrosity and it works, but I am wondering if someone can come up with a cleaner algorithm.
Constraints
To clarify, the algorithm should be constraint to using mostly, or even exclusively, linear algebra operations. This means that most information should be obtained through the use of matrix products, additions, and queries of the properties of the matrices, like the number of zeroes, the trace of the diagonal etc...
Refined implementation
I managed to simplify the logic quite a bit, but I want to see if it is possible to get a more straighforward implementation:
using SBM = Eigen::SparseVector<bool>;
std::array<SBM, 3> Closure(
const Eigen::SparseMatrix<bool>& edge_matrix,
const Eigen::SparseMatrix<bool>& face_matrix,
const std::array<Eigen::SparseVector<bool>, 3>& simplices)
{
const auto face_vertex_adjacency = face_matrix * edge_matrix;
const auto& vertices = simplices[0];
const auto& edges = simplices[1];
const auto& faces = simplice
s[2];
SBM closure1_edges = edges + face_matrix.transpose() * faces;
SBM closure1_vertices = vertices + edge_matrix.transpose() * edges;
return {closure1_vertices, closure1_edges, faces};
}
bool IsComplex(
const Eigen::SparseMatrix<bool>& edge_matrix,
const Eigen::SparseMatrix<bool>& face_matrix,
const std::array<Eigen::SparseVector<bool>, 3>& simplices)
{
const auto& vertices = simplices[0];
const auto& edges = simplices[1];
const auto& faces = simplices[2];
auto[closure_vertices, closure_edges, closure_faces] =
Closure(edge_matrix, face_matrix, {vertices, edges, faces});
return
closure_vertices.isApprox(vertices) &&
closure_edges.isApprox(edges) &&
closure_faces.isApprox(faces);
}
int IsPureComplex(
const Eigen::SparseMatrix<bool>& edge_matrix,
const Eigen::SparseMatrix<bool>& face_matrix,
const std::array<SBM, 3>& simplices)
{
const SBM& verts = simplices[0];
const SBM& edges = simplices[1];
const SBM& faces = simplices[2];
if(!IsComplex(edge_matrix, face_matrix, simplices)) return -1;
SBM point_faces = edges.transpose() * edge_matrix;
const bool has_point_facets = SBM((verts - point_faces).pruned()).nonZeros();
SBM line_faces = faces.transpose() * face_matrix;
const bool has_line_facets = SBM((edges - line_faces).pruned()).nonZeros();
const bool has_triangle_facets = faces.nonZeros();
if(int(has_point_facets) + int(has_line_facets) + int(has_triangle_facets) != 1)
return -1;
return 3 * has_triangle_facets + 2 * has_line_facets + has_point_facets;
}
I added all the code so that this is copy pasteable, but the only important function is the last one IsPureComplex. The above works and is better than what I originally had, but I wonder if it is possible to make it even simpler. That is to say, shorter and with less cyclomatic complexity.
