Rips complex example¶

Todo

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();