#!/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)]

     print '#', time.asctime(), "Points read"

     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), store = (s.dimension() < skeleton))

         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)