--- a/bindings/python/dionysus/__init__.py Fri May 11 17:06:55 2012 -0700
+++ b/bindings/python/dionysus/__init__.py Sun Jun 03 16:39:18 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:39:18 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 Fri May 11 17:06:55 2012 -0700
+++ b/bindings/python/zigzag-persistence.cpp Sun Jun 03 16:39:18 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 Fri May 11 17:06:55 2012 -0700
+++ b/include/topology/image-zigzag-persistence.h Sun Jun 03 16:39:18 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 Fri May 11 17:06:55 2012 -0700
+++ b/include/topology/zigzag-persistence.h Sun Jun 03 16:39:18 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;
};
--- a/include/topology/zigzag-persistence.hpp Fri May 11 17:06:55 2012 -0700
+++ b/include/topology/zigzag-persistence.hpp Sun Jun 03 16:39:18 2012 -0700
@@ -34,7 +34,7 @@
SimplexIndex last_s = visitor.new_simplex(*this);
last_s->low = z_list.end();
#if ZIGZAG_CONSISTENCY
- last_s->boundary = bdry; // NB: debug only
+ last_s->boundary = bdry; // NB: debug only
#endif
rLog(rlZigzagAdd, " Reducing among cycles");
@@ -98,9 +98,9 @@
unsigned order = b_list.empty() ? 0 : boost::prior(b_list.end())->order + 1;
b_list.push_back(BNode(order));
BIndex last_b = boost::prior(b_list.end());
-
+
// Set b_column and low
- last_b->b_column.swap(v);
+ last_b->b_column.swap(v);
last_b->b_column.back()->low = last_b;
// Set b_row
@@ -137,12 +137,12 @@
// Birth
//show_all();
rLog(rlZigzagRemove, "Birth case in remove");
-
+
// Prepend DC[j] = ZB[j] to Z
rLog(rlZigzagRemove, "Computing the column DC[j] = ZB[j] to prepend to Z");
BIndex j = visitor.select_j_in_remove(*this, s->c_row);
rLog(rlZigzagRemove, " j = %d", j->order);
-
+
ZColumn z;
std::for_each(j->b_column.begin(), j->b_column.end(), make_adder(&ZNode::z_column, z));
@@ -157,7 +157,7 @@
// Prepend row of s in C to B
rLog(rlZigzagRemove, "Prepending the row of s to B");
- first_z->b_row = s->c_row; // copying instead of swapping is inefficient,
+ first_z->b_row = s->c_row; // copying instead of swapping is inefficient,
// but it simplifies logic when subtracting chains later
std::for_each(first_z->b_row.begin(), first_z->b_row.end(), make_appender(&BNode::b_column, first_z));
//AssertMsg(check_consistency(), "Must be consistent after prepending row of s to B");
@@ -175,7 +175,7 @@
// Subtract C[j] from every column of C that contains s
AssertMsg(s->c_row == first_z->b_row, "s->c_row == first_z->b_row before subtracting C[j]");
rLog(rlZigzagRemove, "Subtracting C[j]=[%s] from every column of C that contains s=%d with row [%s]",
- j->c_column.tostring(out).c_str(),
+ j->c_column.tostring(out).c_str(),
s->order, s->c_row.tostring(out).c_str());
add_chains(boost::make_filter_iterator(std::bind2nd(NotEqualBIndex(), j), first_z->b_row.begin(), first_z->b_row.end()),
boost::make_filter_iterator(std::bind2nd(NotEqualBIndex(), j), first_z->b_row.end(), first_z->b_row.end()),
@@ -188,7 +188,7 @@
// Subtract B[j] from every other column of B that has l
ZIndex l = j->b_column.back();
BRow l_row = l->b_row;
- rLog(rlZigzagRemove, "Subtracting B[j], j is %d, l is %d, l_row: [%s]",
+ rLog(rlZigzagRemove, "Subtracting B[j], j is %d, l is %d, l_row: [%s]",
j->order, l->order, l_row.tostring(out).c_str());
add_chains(boost::make_filter_iterator(std::bind2nd(NotEqualBIndex(), j), l_row.rbegin(), l_row.rend()),
boost::make_filter_iterator(std::bind2nd(NotEqualBIndex(), j), l_row.rend(), l_row.rend()),
@@ -200,7 +200,7 @@
// Drop j, l, and s
- //
+ //
// l->z_column is the only non-empty thing, but we drop it,
// the basis is preserved because we added first_z
l->z_column.back()->low = z_list.end();
@@ -214,8 +214,8 @@
rLog(rlZigzagRemove, "s->c_row: [%s]", s->c_row.tostring(out).c_str());
if (!s->c_row.empty())
{
- rLog(rlZigzagRemove, "s->c_row[0]: [%s]", s->c_row.front()->c_column.tostring(out).c_str());
- rLog(rlZigzagRemove, " b_column: [%s]", s->c_row.front()->b_column.tostring(out).c_str());
+ rLog(rlZigzagRemove, "s->c_row[0]: [%s]", s->c_row.front()->c_column.tostring(out).c_str());
+ rLog(rlZigzagRemove, " b_column: [%s]", s->c_row.front()->b_column.tostring(out).c_str());
}
AssertMsg(s->c_row.empty(), "c_row of s must be empty before erasing in remove::birth");
visitor.erasing_z(*this, l);
@@ -233,7 +233,7 @@
// if ls->low precedes first_z, swap them
if (cmp(ls->low, first_z)) std::swap(ls->low, first_z);
-
+
add_chain(first_z, ls->low, &ZNode::b_row, &BNode::b_column);
add_chain(ls->low, first_z, &ZNode::z_column, &SimplexNode::z_row);
std::swap(ls->low, first_z);
@@ -254,13 +254,13 @@
rLog(rlZigzagRemove, "j=%d, j->b_row=[%s]", j->order, j->b_row.tostring(out).c_str());
rLog(rlZigzagRemove, "s=%d, s->c_row=[%s]", s->order, s->c_row.tostring(out).c_str());
rLog(rlZigzagRemove, "s=%d, s->z_row=[%s]", s->order, s->z_row.tostring(out).c_str());
-
+
// Subtract Z[j] from every chain in C that contains s
// (it's Ok to go in the forward order since we are subtracting a column in Z from C)
add_chains(c_row.begin(), c_row.end(), j, &BNode::c_column, &SimplexNode::c_row, &ZNode::z_column);
AssertMsg(check_consistency(), "Must be consistent after subtracting Z[j] from C");
-
- // Change basis to remove s from Z
+
+ // Change basis to remove s from Z
// Compute reducers --- columns that we will be adding to other columns
ZRow z_row = s->z_row;
typedef typename ZRow::reverse_iterator ZRowReverseIterator;
@@ -269,36 +269,42 @@
reducers.push_back(boost::prior(z_row.rend())); // j is the first reducer
AssertMsg(*(reducers.back()) == j, "The first element of reducers should be j");
SimplexIndex low = j->z_column.back();
- rLog(rlZigzagRemove, " Added reducer %d [%s] with low=%d",
+ rLog(rlZigzagRemove, " Added reducer %d [%s] with low=%d",
j->order, j->z_column.tostring(out).c_str(),
low->order);
for (typename ZRow::iterator cur = z_row.begin(); cur != z_row.end(); ++cur)
if (cmp((*cur)->z_column.back(), low))
- {
- reducers.push_back(ZRowReverseIterator(boost::next(cur)));
+ {
+ reducers.push_back(ZRowReverseIterator(boost::next(cur)));
low = (*(reducers.back()))->z_column.back();
- rLog(rlZigzagRemove, " Added reducer %d [%s] with low=%d",
+ rLog(rlZigzagRemove, " Added reducer %d [%s] with low=%d",
(*cur)->order, (*cur)->z_column.tostring(out).c_str(),
low->order);
- rLog(rlZigzagRemove, " reducers.back(): %d [%s] with low=%d",
+ rLog(rlZigzagRemove, " reducers.back(): %d [%s] with low=%d",
(*(reducers.back()))->order,
(*(reducers.back()))->z_column.tostring(out).c_str(),
(*(reducers.back()))->z_column.back()->order);
}
- rLog(rlZigzagRemove, " Reducers size: %d, s is %d",
+ rLog(rlZigzagRemove, " Reducers size: %d, s is %d",
reducers.size(), s->order);
//show_all();
// Add each reducer to the columns that follow them until the next reducer
+ // NB: processing reducers in the reverse order fixes a bug in the paper,
+ // in step Remove.Death.1.4, where the matrix B is updated incorrectly.
+ // I can't find a mention of this bug in my notes, but the fact
+ // that it's fixed in the code suggests that I knew about it. (Or
+ // most likely I didn't recognize that what the paper said is not
+ // exactly what I meant.)
typename ReducersContainer::reverse_iterator cur = reducers.rbegin();
ZRowReverseIterator zcur = z_row.rbegin();
-
+
while (cur != reducers.rend())
{
rLog(rlZigzagRemove, " Cur reducer: %d [%s]", (**cur)->order,
(**cur)->z_column.tostring(out).c_str());
- change_basis(zcur, *cur, **cur,
+ change_basis(zcur, *cur, **cur,
&ZNode::z_column, &SimplexNode::z_row,
&ZNode::b_row, &BNode::b_column);
if (cur != reducers.rbegin())
@@ -313,7 +319,7 @@
zcur = *cur++;
// This makes it inconsistent until the next iteration of this update loop
}
-
+
// Drop j and s
Death d = visitor.death(*this, j);
@@ -323,8 +329,8 @@
rLog(rlZigzagRemove, "j->b_row: [%s]", j->b_row.tostring(out).c_str());
if (!j->b_row.empty())
{
- rLog(rlZigzagRemove, "j->b_row[0]: [%s]", j->b_row.front()->b_column.tostring(out).c_str());
- rLog(rlZigzagRemove, " c_column: [%s]", j->b_row.front()->c_column.tostring(out).c_str());
+ rLog(rlZigzagRemove, "j->b_row[0]: [%s]", j->b_row.front()->b_column.tostring(out).c_str());
+ rLog(rlZigzagRemove, " c_column: [%s]", j->b_row.front()->c_column.tostring(out).c_str());
}
AssertMsg(j->b_row.empty(), "b_row of j must be empty before erasing in remove(). Most likely you are trying to remove a simplex whose coface is still in the complex.");
AssertMsg(s->z_row.empty(), "z_row of s must be empty before erasing in remove()");
@@ -336,11 +342,11 @@
//show_all();
AssertMsg(check_consistency(), "Must be consistent when done in remove()");
-
+
return d;
}
}
-
+
template<class BID, class SD>
void
ZigzagPersistence<BID,SD>::
@@ -360,7 +366,7 @@
for (typename CRow::const_iterator ccur = cur->c_row.begin(); ccur != cur->c_row.end(); ++ccur)
std::cout << (*ccur)->order << " ";
std::cout << std::endl;
-
+
std::cout << " low: ";
if (cur->low != z_list.end())
std::cout << cur->low->order;
@@ -368,12 +374,12 @@
std::cout << "none";
std::cout << std::endl;
}
-
+
std::cout << "z_list:" << std::endl;
for (ZIndex cur = z_list.begin(); cur != z_list.end(); ++cur)
{
std::cout << " " << cur->order << ":" << std::endl;
-
+
std::cout << " birth: " << cur->birth << std::endl;
std::cout << " z_column: ";
@@ -387,7 +393,7 @@
std::cout << std::endl;
std::cout << " low: ";
- if (cur->low != b_list.end())
+ if (cur->low != b_list.end())
std::cout << cur->low->order;
else
std::cout << "none";
@@ -438,7 +444,7 @@
}
if (cur->low != z_list.end())
AssertMsg(!(cur->low->z_column.empty()), "z_column must not be empty");
- if (cur->low != z_list.end() && cur->low->z_column.back() != cur)
+ if (cur->low != z_list.end() && cur->low->z_column.back() != cur)
{
rError("low of SimplexNode %d is incorrect", cur->order);
return false;
@@ -462,7 +468,7 @@
rError("In check_consistency(): ZNode %d not found in b_column of %d", cur->order, (*bcur)->order);
return false;
}
- if (cur->low != b_list.end() && cur->low->b_column.back() != cur)
+ if (cur->low != b_list.end() && cur->low->b_column.back() != cur)
{
rError("low of ZNode %d is incorrect", cur->order);
return false;
@@ -523,42 +529,42 @@
/* Private */
// Class: Appender
-//
+//
// Functor that appends given element to the given member of whatever parameter it is invoked with
template<class BID, class SD>
template<class Member, class Element>
struct ZigzagPersistence<BID,SD>::Appender
{
- Appender(Member mm, Element ee):
+ Appender(Member mm, Element ee):
m(mm), e(ee) {}
template<class T>
void operator()(T& a) { ((*a).*m).append(e, cmp); }
-
+
Member m;
Element e;
OrderComparison cmp;
};
// Class: Remover
-//
+//
// Functor that removes given element from the given member of whatever parameter it is invoked with
template<class BID, class SD>
template<class Member, class Element>
struct ZigzagPersistence<BID,SD>::Remover
{
- Remover(Member mm, Element ee):
+ Remover(Member mm, Element ee):
m(mm), e(ee) {}
template<class T>
void operator()(T& a) { ((*a).*m).remove(e); }
-
+
Member m;
Element e;
};
// Class: Adder
-//
+//
// Functor that adds the given member of whatever it is invoked with to the given chain
template<class BID, class SD>
template<class Member, class Chain>
@@ -575,10 +581,10 @@
OrderComparison cmp;
};
-
+
// Function: add_chains()
-//
-// Special case of add_chains where all Indexes are the same, and
+//
+// Special case of add_chains where all Indexes are the same, and
// therefore PrimaryMemberFrom and PrimaryMemberTo are the same
template<class BID, class SD>
template<class Index, class IndexFrom, class PrimaryMember, class SecondaryMember>
@@ -590,7 +596,7 @@
}
// Function: add_chains()
-//
+//
// Adds PrimaryMember pm of j to pm of every element in the range [bg,end)
// Fixes SecondaryMembers by adding and removing the corresponding elements.
// For example, if we add a column to a number of other columns, then PrimaryMember is that
@@ -634,7 +640,7 @@
// Function: add_chain()
//
-// Adds PrimaryMemberFrom pmf of `from` to PrimaryMemberTo pmt of `to`.
+// Adds PrimaryMemberFrom pmf of `from` to PrimaryMemberTo pmt of `to`.
// Fixes SecondaryMemberTos. See add_chains().
template<class BID, class SD>
template<class IndexTo, class IndexFrom, class PrimaryMemberTo, class SecondaryMemberTo, class PrimaryMemberFrom>
@@ -642,8 +648,8 @@
ZigzagPersistence<BID,SD>::
add_chain(IndexTo to, IndexFrom from, PrimaryMemberTo pmt, SecondaryMemberTo smt, PrimaryMemberFrom pmf)
{
- rLog(rlZigzagAddChain, "Adding %d [%s] to %d [%s]",
- (*from).order,
+ rLog(rlZigzagAddChain, "Adding %d [%s] to %d [%s]",
+ (*from).order,
((*from).*pmf).tostring(out).c_str(),
(*to).order,
((*to).*pmt).tostring(out).c_str());
@@ -663,7 +669,7 @@
((*to).*pmt).end(), ((*to).*pmt).end(),
cmp),
make_appender(smt, to));
-
+
// Add primaries
((*to).*pmt).add((*from).*pmf, cmp);
rLog(rlZigzagAddChain, "Got %s", ((*to).*pmt).tostring(out).c_str());
@@ -704,7 +710,7 @@
ZigzagPersistence<BID,SD>::ZigzagVisitor::
new_z_in_remove(ZigzagPersistence& zz)
{
- int order = zz.z_list.empty() ? 0 : zz.z_list.begin()->order - 1;
+ int order = zz.z_list.empty() ? 0 : zz.z_list.begin()->order - 1;
zz.z_list.push_front(ZNode(order, zz.b_list.end()));
return zz.z_list.begin();
}