examples/rips/rips-pairwise.py
author Christos Mantoulidis <cmad@stanford.edu>
Tue, 04 Aug 2009 13:23:16 -0700
branchdev
changeset 156 f75fb57d2831
parent 132 2a737609b8bf
child 145 ee096f207dfb
permissions -rw-r--r--
Changed implementation of WeightedRips to store simplex values (max distance between simplices' vertices) as an invisible layer on top of each simplex object, so that the data() field of WeightedRips has been freed for use by the users again.

#/usr/bin/env python


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


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

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

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 = []
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]

f = Filtration(simplices, 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(f)
print time.asctime(), "Initialized StaticPersistence"

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

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

        if b.dimension() >= skeleton: continue

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

        d = simplices[f[p(i.pair)]]
        print b.dimension(), b.data, d.data