Merged upstream dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Sun, 03 Jun 2012 16:39:18 -0700
branchdev
changeset 254 fb0f9b4ca77a
parent 253 6dc64e820477 (diff)
parent 251 870865d25958 (current diff)
child 255 a7e010314a03
Merged upstream
--- 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();
 }