# Lower-star filtration
elephant_points, elephant_complex = read_off('data/cgal/elephant.off')
show_complex(elephant_points, elephant_complex)

elephant_complex = closure(elephant_complex, 2)
show_complex(elephant_points, elephant_complex)

def axis(points, axis = 1):
    def value(v):
        return points[v][1]

    return value

value = axis(elephant_points, 1)


def max_vertex_compare(points):
    def max_vertex(s):
        return max(value(v) for v in s.vertices)

    def compare(s1, s2):
        return cmp(s1.dimension(), s2.dimension()) or cmp(max_vertex(s1), max_vertex(s2))

    return compare

f = Filtration(elephant_complex, max_vertex_compare(elephant_points))

##persistence = StaticPersistence(f)
persistence = DynamicPersistenceChains(f)
persistence.pair_simplices()

dgms = init_diagrams(persistence, f, lambda s: max(value(v) for v in s.vertices), lambda i: i)

print "Betti numbers:"
for i,dgm in enumerate(dgms):
    print len([p for p in dgm if p[1] == float('inf')])

while True:
    pt = show_diagram(dgms)
    if not pt: break
    print pt
    i = pt[2]
    smap = persistence.make_simplex_map(f)
    chain = [smap[ii] for ii in i.chain]
    pair_cycle = [smap[ii] for ii in i.pair().cycle]
    pair_chain = [smap[ii] for ii in i.pair().chain]
    show_complex(elephant_points, subcomplex = chain)
    show_complex(elephant_points, subcomplex = pair_cycle + pair_chain)
    #show_complex(elephant_points, [s for s in f if s.dimension() != 1], subcomplex = cycle)