--- a/bindings/python/dionysus/__init__.py Mon May 14 10:47:46 2012 -0700
+++ b/bindings/python/dionysus/__init__.py Sun Jun 03 16:17:27 2012 -0700
@@ -1,6 +1,7 @@
from _dionysus import *
from distances import l2, ExplicitDistances, points_file
from zigzag import *
+from adaptor import *
def init_with_none(self, iter, data = None): # convenience: data defaults to None
@@ -52,3 +53,11 @@
res.add(Simplex(face, s.data))
return list(res)
+
+_init_diagrams = init_diagrams
+
+def init_diagrams(p, f, evaluator = None, data = None):
+ if isinstance(p, StaticCohomologyPersistence):
+ return init_diagrams_from_adaptor(p,f, evaluator, data)
+
+ return _init_diagrams(p,f, evaluator, data)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/bindings/python/dionysus/adaptor.py Sun Jun 03 16:17:27 2012 -0700
@@ -0,0 +1,96 @@
+from _dionysus import CohomologyPersistence, PersistenceDiagram
+
+class StaticCohomologyPersistence(object):
+ def __init__(self, filtration, prime = 2, subcomplex = lambda s: True):
+ self.filtration = filtration
+ 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 = 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
+ else:
+ cocycle = self.persistence.__iter__().next()
+ self.pairs[-1][2] = [(n.coefficient, n.si.order) for n in cocycle]
+
+
+ def __iter__(self):
+ for i, (pair, subcomplex, cocycle) in enumerate(self.pairs):
+ if not pair:
+ 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)
+
+ 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 s: 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')
+
+ dgms[dim].append((b,d, data(smap[n])))
+
+ return dgms
--- a/bindings/python/zigzag-persistence.cpp Mon May 14 10:47:46 2012 -0700
+++ b/bindings/python/zigzag-persistence.cpp Sun Jun 03 16:17:27 2012 -0700
@@ -17,23 +17,23 @@
// ZigzagPersistence
bp::tuple zzp_add(dp::ZZPersistence& zzp, bp::object bdry, dp::BirthID birth)
{
- // Make ZColumn
- // NB: it's extremely weird that I have to do it this way,
+ // Make ZColumn
+ // NB: it's extremely weird that I have to do it this way,
// but for some reason I cannot just create boundary on the stack
- boost::shared_ptr<dp::ZZPersistence::ZColumn>
- boundary(new dp::ZZPersistence::ZColumn(bp::stl_input_iterator<dp::ZZPersistence::SimplexIndex>(bdry),
+ boost::shared_ptr<dp::ZZPersistence::ZColumn>
+ boundary(new dp::ZZPersistence::ZColumn(bp::stl_input_iterator<dp::ZZPersistence::SimplexIndex>(bdry),
bp::stl_input_iterator<dp::ZZPersistence::SimplexIndex>()));
boundary->sort(zzp.cmp);
dp::ZZPersistence::SimplexIndex i;
dp::ZZPersistence::Death d;
- boost::tie(i,d) = zzp.add(*boundary, birth);
+ boost::tie(i,d) = zzp.add(*boundary, birth);
return bp::make_tuple(i,d);
}
dp::ZZPersistence::Death zzp_remove(dp::ZZPersistence& zzp, dp::ZZPersistence::SimplexIndex s, dp::ZZPersistence::BirthID birth)
{
- return zzp.remove(s, birth);
+ return zzp.remove(s, birth);
}
bool zzp_is_alive(dp::ZZPersistence& zzp, const dp::ZZPersistence::ZNode& zn)
@@ -44,23 +44,23 @@
// ImageZigzagPersistence
bp::tuple izzp_add(dp::IZZPersistence& izzp, bp::object bdry, bool subcomplex, dp::BirthID birth)
{
- // Make ZColumn
- // NB: it's extremely weird that I have to do it this way,
+ // Make ZColumn
+ // NB: it's extremely weird that I have to do it this way,
// but for some reason I cannot just create boundary on the stack
- boost::shared_ptr<dp::IZZPersistence::ZColumn>
- boundary(new dp::IZZPersistence::ZColumn(bp::stl_input_iterator<dp::IZZPersistence::SimplexIndex>(bdry),
+ boost::shared_ptr<dp::IZZPersistence::ZColumn>
+ boundary(new dp::IZZPersistence::ZColumn(bp::stl_input_iterator<dp::IZZPersistence::SimplexIndex>(bdry),
bp::stl_input_iterator<dp::IZZPersistence::SimplexIndex>()));
boundary->sort(izzp.cmp);
dp::IZZPersistence::SimplexIndex i;
dp::IZZPersistence::Death d;
- boost::tie(i,d) = izzp.add(*boundary, subcomplex, birth);
+ boost::tie(i,d) = izzp.add(*boundary, subcomplex, birth);
return bp::make_tuple(i,d);
}
dp::IZZPersistence::Death izzp_remove(dp::IZZPersistence& izzp, dp::IZZPersistence::SimplexIndex s, dp::IZZPersistence::BirthID birth)
{
- return izzp.remove(s, birth);
+ return izzp.remove(s, birth);
}
@@ -79,12 +79,14 @@
}
// ZNode
-dp::ZZPersistence::ZColumn::const_iterator
-znode_zcolumn_begin(dp::ZZPersistence::ZNode& zn)
+template<class Persistence>
+typename Persistence::ZColumn::const_iterator
+znode_zcolumn_begin(typename Persistence::ZNode& zn)
{ return zn.z_column.begin(); }
-dp::ZZPersistence::ZColumn::const_iterator
-znode_zcolumn_end(dp::ZZPersistence::ZNode& zn)
+template<class Persistence>
+typename Persistence::ZColumn::const_iterator
+znode_zcolumn_end(typename Persistence::ZNode& zn)
{ return zn.z_column.end(); }
@@ -95,9 +97,11 @@
.def("order", &si_order<dp::ZZPersistence::SimplexIndex>)
.def("__repr__", &si_repr<dp::ZZPersistence::SimplexIndex>)
;
-
+
bp::class_<dp::IZZPersistence::SimplexIndex>("IZZSimplexIndex")
- .def("order", &si_order<dp::IZZPersistence::SimplexIndex>);
+ .def("order", &si_order<dp::IZZPersistence::SimplexIndex>)
+ .def("__repr__", &si_repr<dp::IZZPersistence::SimplexIndex>)
+ ;
bp::class_<dp::ZZPersistence>("ZigzagPersistence")
.def("add", &zzp_add)
@@ -105,14 +109,20 @@
.def("is_alive", &zzp_is_alive)
.def("__iter__", bp::range(&dp::ZZPersistence::begin, &dp::ZZPersistence::end))
;
-
+
bp::class_<dp::IZZPersistence>("ImageZigzagPersistence")
.def("add", &izzp_add)
.def("remove", &izzp_remove)
+ .def("__iter__", bp::range(&dp::IZZPersistence::image_begin, &dp::IZZPersistence::image_end))
;
bp::class_<dp::ZZPersistence::ZNode>("ZNode", bp::no_init)
.add_property("birth", &dp::ZZPersistence::ZNode::birth)
- .def("__iter__", bp::range(&znode_zcolumn_begin, &znode_zcolumn_end))
+ .def("__iter__", bp::range(&znode_zcolumn_begin<dp::ZZPersistence>, &znode_zcolumn_end<dp::ZZPersistence>))
+ ;
+
+ bp::class_<dp::IZZPersistence::ZNode>("IZNode", bp::no_init)
+ .add_property("birth", &dp::IZZPersistence::ZNode::birth)
+ .def("__iter__", bp::range(&znode_zcolumn_begin<dp::IZZPersistence>, &znode_zcolumn_end<dp::IZZPersistence>))
;
}
--- a/include/topology/image-zigzag-persistence.h Mon May 14 10:47:46 2012 -0700
+++ b/include/topology/image-zigzag-persistence.h Sun Jun 03 16:17:27 2012 -0700
@@ -6,7 +6,7 @@
struct SimplexSubcomplexData
{
- SimplexSubcomplexData(bool sc = false):
+ SimplexSubcomplexData(bool sc = false):
subcomplex(sc) {}
bool subcomplex;
@@ -40,10 +40,10 @@
cok_order_begin(std::numeric_limits<int>::max()/2)
{}
- IndexDeathPair add(ZColumn bdry,
- bool subcomplex,
+ IndexDeathPair add(ZColumn bdry,
+ bool subcomplex,
const BirthID& birth = BirthID()) { ImageZZVisitor zzv(subcomplex); return Parent::add(bdry, birth, zzv); }
-
+
Death remove(SimplexIndex s, const BirthID& birth = BirthID()) { ImageZZVisitor zzv(s->subcomplex); return Parent::remove(s, birth, zzv); }
@@ -59,31 +59,31 @@
public:
ImageZZVisitor(bool sc = false):
subcomplex(sc), birth_in_image(false) {}
-
- // Sets the subcomplex property of the new simplex
+
+ // Sets the subcomplex property of the new simplex
SimplexIndex new_simplex(Parent& zz);
-
+
// Decides where to put the new column (image or cokernel)
ZIndex new_z_in_add(Parent& zz, const ZColumn& z, const BRow& u);
-
+
// Checks if there is a boundary entirely in the subcomplex, and sets birth_in_image accordingly
BIndex select_j_in_remove(Parent& zz, const CRow& c_row);
-
+
ZIndex new_z_in_remove(Parent& zz);
-
+
// Updates im_last and cok_begin if necessary
void erasing_z(Parent& zz, ZIndex j);
// Determines if there is a death in the image
Death death(Parent& zz, ZIndex dying_z);
-
+
private:
ZIndex append_in_image(Parent& zz);
ZIndex append_in_cokernel(Parent& zz);
ZIndex prepend_in_image(Parent& zz);
ZIndex prepend_in_cokernel(Parent& zz);
bool in_subcomplex(const ZColumn& z);
-
+
bool subcomplex, birth_in_image;
};
--- a/include/topology/zigzag-persistence.h Mon May 14 10:47:46 2012 -0700
+++ b/include/topology/zigzag-persistence.h Sun Jun 03 16:17:27 2012 -0700
@@ -50,7 +50,7 @@
struct ZNode
{
- ZNode(int o, BIndex l):
+ ZNode(int o, BIndex l):
order(o), low(l) {}
int order;
@@ -58,7 +58,7 @@
BRow b_row;
BIndex low; // which BColumn has this ZIndex as low
- BirthID birth; // TODO: do we need to do empty-member optimization?
+ BirthID birth; // TODO: do we need to do empty-member optimization?
// i.e., does it ever make sense for birth to be empty?
};
@@ -73,7 +73,7 @@
struct SimplexNode: public SimplexData
{
- SimplexNode(unsigned o, ZIndex l):
+ SimplexNode(unsigned o, ZIndex l):
order(o), low(l) {}
unsigned order;
@@ -87,10 +87,10 @@
// Constructor: ZigzagPersistence()
ZigzagPersistence() {}
-
+
// Function: add(bdry, birth)
IndexDeathPair add(ZColumn bdry, const BirthID& birth = BirthID()) { ZigzagVisitor zzv; return add<ZigzagVisitor>(bdry, birth, zzv); }
-
+
// Function: remove(s, birth)
Death remove(SimplexIndex s, const BirthID& birth = BirthID()) { ZigzagVisitor zzv; return remove<ZigzagVisitor>(s, birth, zzv); }
@@ -100,8 +100,8 @@
bool is_alive(ZIndex i) const { return i->low == b_list.end(); }
bool is_alive(ZNode zn) const { return zn.low == b_list.end(); }
-
- protected:
+
+ protected:
// Function: add(s)
template<class Visitor>
IndexDeathPair add(ZColumn bdry, const BirthID& birth, Visitor& visitor);
@@ -112,7 +112,7 @@
// Struct: ZigzagVisitor
// Various methods of an instance of this class are called at different stages of addition and removal algorithm.
- // NB: currently the places where it's called are catered for image zigzags, in the future this could be expanded
+ // NB: currently the places where it's called are catered for image zigzags, in the future this could be expanded
// to provide simple support for other algorithms
// TODO: not obvious that the methods should be const (and therefore the reference passed to add() and remove())
// revisit when working on ImageZigzag
@@ -123,7 +123,7 @@
// Function: new_z_in_add(zz, z, u)
// Called when a new cycle is born after adding a simplex. The method is expected to add an element to z_list, and return its ZIndex.
ZIndex new_z_in_add(ZigzagPersistence& zz, const ZColumn& z, const BRow& u);
-
+
BIndex select_j_in_remove(ZigzagPersistence& zz, const CRow& c_row);
ZIndex new_z_in_remove(ZigzagPersistence& zz);
@@ -134,7 +134,7 @@
};
public:
- // Debug; non-const because Indices are iterators, and not const_iterators
+ // Debug; non-const because Indices are iterators, and not const_iterators
void show_all();
bool check_consistency(SimplexIndex s_skip, ZIndex z_skip, BIndex b_skip);
bool check_consistency() { return check_consistency(s_list.end(), z_list.end(), b_list.end()); }
@@ -159,31 +159,31 @@
template<class Index, class IndexFrom, class PrimaryMember, class SecondaryMember>
void add_chains(Index bg, Index end, IndexFrom j, PrimaryMember pm, SecondaryMember sm);
template<class IndexTo, class IndexFrom, class PrimaryMemberTo, class SecondaryMemberTo, class PrimaryMemberFrom>
- void add_chains(IndexTo bg, IndexTo end, IndexFrom j,
+ void add_chains(IndexTo bg, IndexTo end, IndexFrom j,
PrimaryMemberTo pmt, SecondaryMemberTo smt,
PrimaryMemberFrom pmf);
template<class Index, class PrimaryMember, class SecondaryMember>
- void add_chain(Index to, Index from,
+ void add_chain(Index to, Index from,
PrimaryMember pmt, SecondaryMember smt);
template<class IndexTo, class IndexFrom, class PrimaryMemberTo, class SecondaryMemberTo, class PrimaryMemberFrom>
- void add_chain(IndexTo to, IndexFrom from,
+ void add_chain(IndexTo to, IndexFrom from,
PrimaryMemberTo pmt, SecondaryMemberTo smt,
PrimaryMemberFrom pmf);
template<class IndexTo, class IndexFrom, class PrimaryMember, class SecondaryMember, class DualPrimaryMember, class DualSecondaryMember>
- void change_basis(IndexTo bg, IndexTo end, IndexFrom j,
- PrimaryMember pm, SecondaryMember sm,
+ void change_basis(IndexTo bg, IndexTo end, IndexFrom j,
+ PrimaryMember pm, SecondaryMember sm,
DualPrimaryMember dpm, DualSecondaryMember dsm);
public:
struct OrderComparison
{
- template<class T>
+ template<class T>
bool operator()(T a, T b) const { return a->order < b->order; }
} cmp;
-
+
struct OrderOutput
{
- template<class T>
+ template<class T>
std::string operator()(T a) const { std::stringstream s; s << a->order; return s.str(); }
} out;
};