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