Rips complex example

Todo

Explain Vietoris-Rips complex.

There is an elementary example in examples/rips/rips.py that computes a Rips complex of a few points with integer coordinates on a line. It illustrates the use of Rips complexes and in particular of defining your own notion of a Distance based on which the Rips complex is constructed.

A more useful example is given in examples/rips/rips-pairwise.py (and its C++ counterpart in examples/rips/rips-pairwise.cpp). The example takes on the command line the filename of a file with points in Euclidean space (one point per line), and a cut off parameters for the skeleton and the \epsilon parameter for the Rips complex construction. It then constructs the Rips complex up to these cutoff parameters, computes its persistence, and outputs the persistence diagram (one point per line).

#/usr/bin/env python


from    dionysus    import Rips, PairwiseDistances, StaticPersistence, Filtration, points_file, \
                           ExplicitDistances, data_dim_cmp
from    sys         import argv, exit
import  time

def main(filename, skeleton, max):
    points = [p for p in points_file(filename)]
    distances = PairwiseDistances(points)
    # distances = ExplicitDistances(distances)           # speeds up generation of the Rips complex at the expense of memory usage
    rips = Rips(distances)
    print time.asctime(), "Rips initialized"

    simplices = Filtration()
    rips.generate(skeleton, max, simplices.append)
    print time.asctime(), "Generated complex: %d simplices" % len(simplices)

    # While this step is unnecessary (Filtration below can be passed rips.cmp), 
    # it greatly speeds up the running times
    for s in simplices: s.data = rips.eval(s)
    print time.asctime(), simplices[0], '...', simplices[-1]

    simplices.sort(data_dim_cmp)             # could be rips.cmp if s.data for s in simplices is not set
    print time.asctime(), "Set up filtration"

    p = StaticPersistence(simplices)
    print time.asctime(), "Initialized StaticPersistence"

    p.pair_simplices()
    print time.asctime(), "Simplices paired"

    print "Outputting persistence diagram"
    smap = p.make_simplex_map(simplices)
    for i in p:
        if i.sign():
            b = smap[i]

            if b.dimension() >= skeleton: continue

            if i.unpaired():
                print b.dimension(), b.data, "inf"
                continue

            d = smap[i.pair()]
            print b.dimension(), b.data, d.data

if __name__ == '__main__':
    if len(argv) < 4:
        print "Usage: %s POINTS SKELETON MAX" % argv[0]
        exit()

    filename = argv[1]
    skeleton = int(argv[2])
    max = float(argv[3])

    main(filename, skeleton, max)

The bit that sets up the Rips complex is:

distances = PairwiseDistances(points)
rips = Rips(distances)
simplices = Filtration()
rips.generate(skeleton, max, simplices.append)

The computation of persistence and output of the persistence diagram are the same as in the Alpha shape example. The example also incorporates the Speed-up suggestions given in the Brief Tutorial.

C++ sketch

Warning

This section is not finished.

The example given in examples/rips/rips.cpp illustrates how one can use the library to compute persistence of a Vietoris-Rips complex for a given set of distances. At the top of the file a struct Distances is defined. The particular distances in the example are trivial (they are points with integer coordinates on a real line), however, the struct illustrates the basic requirements of any such class to be passed to the Rips<Distances> class.

The Rips complex itself is generated in the line:

rips.generate(2, 50, make_push_back_functor(complex));

which tells it to generate a 2-skeleton of the Rips complex up to distance value of 50, and insert the simplices into the previously defined vector complex.

Subsequent sort is unnecessary since Bron-Kerbosch algorithm that generates the complex will actually generate the simplices in lexicographic order; it’s there for illustration purposes only (the simplices must be sorted lexicographically).

The following “paragraph” sets up the filtration with respect to simplex sizes (specified by Generator::Comparison(distances)), and computes its persistence:

// Generate filtration with respect to distance and compute its persistence
Fltr f(complex.begin(), complex.end(), Generator::Comparison(distances));
Persistence p(f);
p.pair_simplices();