examples/cohomology/rips-pairwise-cohomology.py
 author Dmitriy Morozov Mon, 10 Aug 2009 15:22:32 -0700 branch dev changeset 159 92318b22ffa5 parent 146 4e27f1f7c169 child 173 5fd3f43e6fbf permissions -rw-r--r--
Fixed bugs in rips-pairwise-cohomology.py whose origin is unclear to me
```
#!/usr/bin/env python

from    dionysus        import Simplex, CohomologyPersistence, points_file, PairwiseDistances, ExplicitDistances, Rips, data_dim_cmp
from    sys             import argv, exit
import  time

def main(filename, skeleton, max, prime = 11):
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]

simplices.sort(data_dim_cmp)
print '#', time.asctime(), "Simplices sorted"

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

for ccl in ch:
dimension, birth = ccl.birth
if dimension >= skeleton: continue
print dimension, birth, 'inf'         # dimension, simplex data = birth
print "# Cocycle at (dim=%d, birth=%f)" % ccl.birth
for e in ccl:
print "#  ", e.si.order, normalized(e.coefficient, prime)

def normalized(coefficient, prime):
if coefficient > prime/2:
return coefficient - prime
return coefficient

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

filename = argv[1]
skeleton = int(argv[2])
max = float(argv[3])
prime = (len(argv) > 4 and argv[4]) or 11

main(filename, skeleton, max, prime)
```