Using Aleph to calculate the homology of 2-manifolds

Tags: programming, projects, research

In a previous article, I discussed the distribution of Betti numbers of triangulations of 2-manifolds. This article now discusses the code. Coincidentally, it also demonstrates how to use Aleph, my (forthcoming) library for exploring persistent homology and its uses.

Aleph is header-only, so it can be readily deployed in your own projects. In the following, I am assuming that you have installed Aleph, following the official instructions. There are several tutorials available that cover different parts of the persistent homology calculation process. For our purpose, viz. the calculation of homology of triangulated spaces, no large machinery is needed. We start with some includes and using directives:

#include <aleph/topology/io/LexicographicTriangulation.hh>

#include <aleph/persistentHomology/Calculation.hh>

#include <aleph/topology/Simplex.hh>
#include <aleph/topology/SimplicialComplex.hh>

#include <iostream>
#include <string>
#include <vector>

using DataType          = bool;
using VertexType        = unsigned short;

using Simplex           = aleph::topology::Simplex<DataType, VertexType>;
using SimplicialComplex = aleph::topology::SimplicialComplex<Simplex>;

Nothing fancy happened so far. We set up a simplex class, which is the basic data type in most applications, as well as a simplicial complex, which is the basic storage class for multiple simplices. The only interesting thing is our choice of DataType and VertexType. Since our simplices need not store any additional data except for their vertices, we select bool as a data type in order to make the class smaller. This could potentially be achieved with EBCO as well, but I did not yet have time to test it adequately. In addition to the data type, we use unsigned short as the vertex type—the triangulations that we want to analyse only feature 9 or 10 vertices, so unsigned short is the best solution for storing vertex identifiers.

Next, we need some I/O code to read a lexicographic triangulation:

int main(int argc, char* argv[])
{
  if( argc <= 1 )
    return -1;

  std::string filename = argv[1];
  std::vector<SimplicialComplex> simplicialComplexes;

  aleph::topology::io::LexicographicTriangulationReader reader;
  reader( filename, simplicialComplexes );

  for( auto&& K : simplicialComplexes )
  {
    K.createMissingFaces();
    K.sort();
  }
}

Here, we used the LexicographicTriangulationReader, a reader class that supports reading files in the format defined by Frank H. Lutz. However, this format only provides the top-level simplices of a triangulation. Hence, for a 2-manifold, only the 2-simplices are specified. For calculating homology groups, however, all simplices are required. Luckily, the SimplicialComplex class offers a method for just this purpose. By calling createMissingFaces(), all missing faces of the simplicial complex are calculated and added to the simplicial complex. Afterwards, we use sort() to sort simplices in lexicographical order. This order is required to ensure that the homology groups can be calculated correctly—the calculation routines assume that the complex is being "filtrated", so faces need to precede co-faces.

Having created and stored a list of simplicial complexes, we may now finally calculate their homology groups by adding the following code after the last for-loop:

for( auto&& K : simplicialComplexes )
{
  bool dualize                    = true;
  bool includeAllUnpairedCreators = true;

  auto diagrams
    = aleph::calculatePersistenceDiagrams( K,
                                          dualize,
                                          includeAllUnpairedCreators );
}

This code calls calculatePersistenceDiagrams(), which is usually employed to calculate, well, a set of persistence diagrams. The two flags dualize and includeAllUnpairedCreators also deserve some explanation. The first flag only instructs the convenience function as to whether the boundary matrix that is required for (persistent) homology calculations should be dualized or not. Dualization was shown to result in faster computations, so of course we want to use it as well. The second parameter depends on our particular application scenario. Normally, the persistent homology calculation ignores all creator simplices in the top dimension. The reason for this is simple: if we expand a neighbourhood graph to a Vietoris–Rips complex, the top-level simplices are an artefact of the expansion process. Most of these simplices cannot be paired with higher-dimensional simplices, hence they will appear to create a large number of high-dimensional holes. The persistent homology calculation convenience function hence ignores those simplices for the creation of persistence diagrams by default. For our application, however, keeping those simplices is actually the desired behaviour—the top-level simplices are well-defined and an integral part of the triangulation.

As a last step, we want to print the signature of Betti numbers for the given simplicial complex. This is easily achieved by adding a nested loop at the end of the for-loop:

for( auto&& diagram : diagrams )
{
  auto d = diagram.dimension();
  std::cout << "b_" << d << " = " << diagram.betti() << "\n";
}

This will give us all non-zero Betti numbers of the given triangulation. The output format is obviously not optimal—in a research setting, I like to use JSON for creating human-readable and machine-readable results. If you are interested in how this may look, please take a look at a more involved tool for dealing with triangulations. Be aware, however, that this tool is still a work in progress.

This concludes our little foray into Aleph. I hope that I piqued your interest! By the way, I am always looking for contributors…