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.