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