Added a note about reducing matrix Z in i/t/zigzag-persistence.hpp dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Mon, 14 May 2012 10:47:46 -0700
branchdev
changeset 252 4de0bd5b4ae1
parent 250 021030a8f97c
child 253 6dc64e820477
Added a note about reducing matrix Z in i/t/zigzag-persistence.hpp
include/topology/zigzag-persistence.hpp
--- a/include/topology/zigzag-persistence.hpp	Fri May 11 10:25:57 2012 -0700
+++ b/include/topology/zigzag-persistence.hpp	Mon May 14 10:47:46 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();
 }