The procedure described below is explained in detail in [dSVJ09].
One can use examples/cohomology/rips-pairwise-cohomology.cpp to compute persistent pairing of the Rips filtration using the persistent cohomology algorithm. It takes as input a file containing a point set in Euclidean space (one per line) as well as the following command-line flags:
The prime to use in the computation (defaults to 11).
Maximum cutoff parameter up to which to compute the complex.
Skeleton to compute; persistent pairs output will be this number minus 1 (defaults to 2).
Filename where to output the boundary matrix.
Prefix of the filenames where to output the 1-dimensional cocycles.
Filename where to output the simplex vertex mapping.
Filename where to output the persistence diagram.
For example:
rips-pairwise-cohomology points.txt -m 1 -b points.bdry -c points -v points.vrt -d points.dgm
Assuming that at the threshold value of 1 (-m 1 above) Rips complex contains 1-dimensional cocycles, they will be output into filenames of the form points-0.ccl, points-1.ccl, etc.
Subsequently one can use examples/cohomology/cocycle.py to assign to each vertex of the input point set a circle-valued function. It takes the boundary matrix, cocycle, and simplex-vertex map as an input (all produced at the previous step):
cocycle.py points.bdry points-0.ccl points.vrt
The above command outputs a file points-0.val which contains values assigned to the input points (the lines match the lines of the input file points.txt, but also contains the indices).
Two auxilliary tools allow one to visualize the values assigned to the points (using Matplotlib): tools/plot-values/plot.py and tools/plot-values/scatter.py:
plot.py points-0.val points.txt scatter.py points-0.val points-1.val
The Python LSQR code (ported from the Stanford MATLAB implementation to Python by Jeffery Kline) included with Dionysus, and used in examples/cohomology/cocycle.py, requires CVXOPT.
examples/cohomology/rips-pairwise-cohomology.py gives an example of the same computation performed in Python (but with the output in a different format).
After the simplicial complex is computed in a list simplices, and the list is sorted with respect to the Rips filtration order, the simplices are inserted into the CohomologyPersistence one by one:
# list simplices is created
ch = CohomologyPersistence(prime)
complex = {}
for s in simplices:
i,d = ch.add([complex[sb] for sb in s.boundary], (s.dimension(), s.data))
complex[s] = i
if d:
dimension, birth = d
print dimension, birth, s.data
# else birth
Above dictionary complex maintains the map of simplices to indices returned by CohomologyPersistence.add(). The pair (dimension, data) is used as the birth value. Here data is the value associated with the simplex in the Rips filtration. The pair is returned back if a death occurs, and is printed on the standard output. After the for loop finishes, one may output infinite persistence classes with the following for loop:
for ccl in ch:
dimension, birth = ccl.birth
if dimension >= skeleton: continue
print dimension, birth, 'inf' # dimension, simplex data = birth
Naturally one may iterate over ccl which is of type Cocycle and extract more information. For example, this is necessary to get the coefficients that serve as the input for examples/cohomology/cocycle.py.