bindings/python/dionysus/adaptor.py
author "Arnur Nigmetov <a.nigmetov@gmail.com>"
Mon, 02 Mar 2015 13:20:51 +0100
branchdev
changeset 281 6e883f004ebe
parent 267 2f02384a4d9b
permissions -rw-r--r--
Added executable to compute Wasserstein distance. Powering to 1/q is done in output. Maybe, should be changed in the function itself?

from _dionysus import CohomologyPersistence, PersistenceDiagram, Cocycle

class StaticCohomologyPersistence(object):
    def __init__(self, filtration, prime = 2, subcomplex = lambda s: True):
        self.filtration = filtration
        self.prime = prime
        self.subcomplex = subcomplex
        self.persistence = CohomologyPersistence(prime)
        self.pairs = []

    def pair_simplices(self):
        indices = []
        for i,s in enumerate(self.filtration):
            sc = self.subcomplex(s)
            boundary = (indices[self.filtration(ss)] for ss in s.boundary)
            idx,d,ccl = self.persistence.add(boundary, i, image = sc)
            indices.append(idx)
            self.pairs.append([i, sc, []])
            if d:                           # Death
                if self.pairs[d][1]:        # Birth was in the subcomplex
                    self.pairs[i][0] = d    # i killed d
                    self.pairs[d][0] = i    # d was killed by i
                    self.pairs[d][2] = self._cocycle_list(ccl)  # record the cocycle at the time of death
            else:
                cocycle = self.persistence.__iter__().next()
                self.pairs[-1][2] = cocycle

        # Replace cocycles with lists
        for i in xrange(len(self.pairs)):
            ccl = self.pairs[i][2]
            if isinstance(ccl, Cocycle):
                self.pairs[i][2] = self._cocycle_list(ccl)

    def _cocycle_list(self, ccl):
        return [(n.coefficient if n.coefficient <= self.prime/2 else n.coefficient - self.prime, n.si.order) for n in ccl]

    def __call__(self, n):
        return n.i

    def __len__(self):
        return len(self.pairs)

    def __iter__(self):
        for i, (pair, subcomplex, cocycle) in enumerate(self.pairs):
            if pair == i:       # unpaired
                if subcomplex:
                    yield APNode(i, self.pairs)
            else:
                if pair > i and subcomplex:
                    yield APNode(i, self.pairs)
                elif pair < i:
                    pair_pair, pair_subcomplex, pair_cocycle = self.pairs[pair]
                    if pair_subcomplex:
                        yield APNode(i, self.pairs)

    def make_simplex_map(self, filtration):
        return APSimplexMap(filtration)

class ImagePersistence(StaticCohomologyPersistence):
    def __init__(self, filtration, subcomplex):
        super(ImagePersistence, self).__init__(filtration, subcomplex = subcomplex)

# Remaps APNodes into Simplices
class APSimplexMap:
    def __init__(self, filtration):
        self.filtration = filtration

    def __getitem__(self, n):
        return self.filtration[n.i]

class APNode:
    def __init__(self, i, pairs):
        self.i = i
        self.pairs = pairs

    def sign(self):
        return self.unpaired() or self.i < self._pair()

    def unpaired(self):
        return self.i == self._pair()

    def _pair(self):
        return self.pairs[self.i][0]

    def pair(self):
        return APNode(self._pair(), self.pairs)

    @property
    def cocycle(self):
        return self.pairs[self.i][2]

def init_diagrams_from_adaptor(p, f, evaluator, data):
    if not evaluator:
        evaluator = lambda s: s.data

    if not data:
        data = lambda i: None

    dgms = []
    smap = p.make_simplex_map(f)
    for n in p:
        if not n.sign(): continue

        dim = smap[n].dimension()
        if dim + 1 > len(dgms):
            dgms.append(PersistenceDiagram(dim))

        b = evaluator(smap[n])
        d = evaluator(smap[n.pair()]) if not n.unpaired() else float('inf')
        if b == d: continue

        dgms[dim].append((b,d, data(n)))

    return dgms