#include <utilities/log.h>
#include <boost/utility.hpp>
#include <boost/iterator/filter_iterator.hpp>
#include <algorithm>
#include <utilities/indirect.h>
#include <functional>

#ifdef LOGGING
static rlog::RLogChannel* rlZigzagAdd =                   DEF_CHANNEL("topology/persistence/zigzag/add",        rlog::Log_Debug);
static rlog::RLogChannel* rlZigzagRemove =                DEF_CHANNEL("topology/persistence/zigzag/remove",     rlog::Log_Debug);
static rlog::RLogChannel* rlZigzagAddChain =              DEF_CHANNEL("topology/persistence/zigzag/addchain",   rlog::Log_Debug);
static rlog::RLogChannel* rlZigzagCheckConsistency=       DEF_CHANNEL("topology/persistence/zigzag/check",      rlog::Log_Debug);
#endif // LOGGING

#ifdef COUNTERS
static Counter*  cZigzagAdd =                             GetCounter("zigzag/add");
static Counter*  cZigzagRemove =                          GetCounter("zigzag/remove");
static Counter*  cZigzagConsistency =                     GetCounter("zigzag/consistency");
#endif // COUNTERS

template<class BID, class SD>
template<class Visitor>
typename ZigzagPersistence<BID,SD>::IndexDeathPair
ZigzagPersistence<BID,SD>::
add(ZColumn bdry, const BirthID& birth, Visitor& visitor)
{
    Count(cZigzagAdd);

    rLog(rlZigzagAdd,       "Entered ZigzagPersistence::add()");
    rLog(rlZigzagAdd,       "  Boundary: %s", bdry.tostring(out).c_str());
    rLog(rlZigzagAdd,       "  Boundary size: %d", bdry.size());
    AssertMsg(check_consistency(), "Must be consistent before addition");

    SimplexIndex last_s     = visitor.new_simplex(*this);
    last_s->low             = z_list.end();
#if ZIGZAG_CONSISTENCY
    last_s->boundary        = bdry;     // NB: debug only    
#endif

    rLog(rlZigzagAdd,   "  Reducing among cycles");
    // Reduce bdry among the cycles
    rLog(rlZigzagAdd,       "    Boundary: %s", bdry.tostring(out).c_str());
    BColumn v;                // representation of the boundary in the cycle basis
    while (!bdry.empty())
    {
        SimplexIndex l      = bdry.back();
        ZIndex k            = l->low;
        v.append(k, cmp);
        bdry.add(k->z_column, cmp);
        rLog(rlZigzagAdd,       "    Boundary: %s", bdry.tostring(out).c_str());
    }
    rLog(rlZigzagAdd,   "  Reduced among cycles");

    // Reduce v among boundaries
    BRow u;
    while (!(v.empty()))
    {
        ZIndex l = v.back();
        BIndex k = l->low;

        if (k == b_list.end())
            break;

        u.append(k, cmp);
        v.add(k->b_column, cmp);
    }
    rLog(rlZigzagAdd,   "  Reduced among boundaries");

    if (v.empty())
    {
        rLog(rlZigzagAdd,       "  Birth case in add");

        // Figure out the new cycle z
        ZColumn z;
        std::for_each(u.begin(), u.end(), make_adder(&BNode::c_column, z));
        z.append(last_s, cmp);

        // Birth
        ZIndex new_z                = visitor.new_z_in_add(*this, z, u);
        new_z->birth                = birth;

        // Set s_row
        std::for_each(z.begin(), z.end(), make_appender(&SimplexNode::z_row, new_z));

        // Set z_column
        new_z->z_column.swap(z);

        // Set low
        new_z->low                  = b_list.end();
        last_s->low                 = new_z;

        return std::make_pair(last_s, Death());
    } else
    {
        rLog(rlZigzagAdd,       "  Death case in add");

        // Death
        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.back()->low = last_b;

        // Set b_row
        std::for_each(last_b->b_column.begin(), last_b->b_column.end(), make_appender(&ZNode::b_row, last_b));

        // Set c_column
        CColumn& c                  = last_b->c_column;
        std::for_each(u.begin(), u.end(), make_adder(&BNode::c_column, c));
        c.append(last_s, cmp);

        // Set c_row
        std::for_each(c.begin(), c.end(), make_appender(&SimplexNode::c_row, last_b));

        return std::make_pair(last_s, visitor.death(*this, last_b->b_column.back()));
    }
}


template<class BID, class SD>
template<class Visitor>
typename ZigzagPersistence<BID,SD>::Death
ZigzagPersistence<BID,SD>::
remove(SimplexIndex s, const BirthID& birth, Visitor& visitor)
{
    Count(cZigzagRemove);

    rLog(rlZigzagRemove,        "Entered ZigzagPersistence::remove(%d)", s->order);
    AssertMsg(check_consistency(), "Must be consistent before removal");

    if (s->z_row.empty())
    {
        AssertMsg(!(s->c_row.empty()),  "Birth after removal, row in C must be non-empty");

        // 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));

        ZIndex first_z              = visitor.new_z_in_remove(*this);
        first_z->birth              = birth;
        std::for_each(z.begin(),           z.end(),           make_appender(&SimplexNode::z_row, first_z));
        first_z->z_column.swap(z);
        first_z->low                = b_list.end();

        rLog(rlZigzagRemove,        "  Prepended %d [%s]", first_z->order, z.tostring(out).c_str());
        //AssertMsg(check_consistency(),  "Must be consistent after prepending DC[j] = ZB[j] to Z");

        // 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, 
                                        // 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");

#if ZIGZAG_CONSISTENCY
        {
            ZColumn zz;
            std::for_each(j->b_column.begin(), j->b_column.end(), make_adder(&ZNode::z_column, zz));
            AssertMsg(zz.empty(),       "ZB[j] must be 0 after we prepended the row of s in C to B");
        }
#endif

        typedef         std::not_equal_to<BIndex>       NotEqualBIndex;

        // 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(), 
                                    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()),
                   j, &BNode::c_column, &SimplexNode::c_row);
        add_chain(j, j, &BNode::c_column, &SimplexNode::c_row);
        // TODO: that's how it was done before, now it can be removed
        //       add_chains(first_z->b_row.rbegin(), first_z->b_row.rend(), j, &BNode::c_column, &SimplexNode::c_row);
        //AssertMsg(check_consistency(s_list.end(), z_list.begin(), b_list.end()),  "Must be consistent after subtracting C[j] in remove::birth");

        // 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]", 
                                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()),
                   j, &BNode::b_column, &ZNode::b_row);
        j->b_column.back()->low = b_list.end();     // redundant since l will be deleted (here for check_consistency only)
        add_chain(j, j, &BNode::b_column, &ZNode::b_row);
        // TODO: investigate why this works for ordinary zigzag, but fails for the image zigzag
        //AssertMsg(check_consistency(s_list.end(), first_z, b_list.end()),  "Must be consistent after subtracting B[j] in remove::birth");


        // 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();
        std::for_each(l->z_column.begin(), l->z_column.end(), make_remover(&SimplexNode::z_row, l));

        //show_all();
        rLog(rlZigzagRemove,        "l=%d has z_column: [%s]", l->order, l->z_column.tostring(out).c_str());

        AssertMsg(l->b_row.empty(),     "b_row of l must be empty before erasing in remove::birth");
        AssertMsg(s->z_row.empty(),     "z_row of s must be empty before erasing in remove::birth");
        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()); 
        }
        AssertMsg(s->c_row.empty(),     "c_row of s must be empty before erasing in remove::birth");
        visitor.erasing_z(*this, l);
        b_list.erase(j);
        z_list.erase(l);
        s_list.erase(s);
        AssertMsg(check_consistency(s_list.end(), first_z, b_list.end()),  "Must be consistent before reducing Z in remove::birth");

        // Reduce Z
        rLog(rlZigzagRemove,        "Reducing Z");
        SimplexIndex ls = first_z->z_column.back();
        while(ls->low != first_z)
        {
            if (ls->low == z_list.end())    { ls->low = first_z; break; }

            // 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);

            ls = first_z->z_column.back();
        }
        AssertMsg(check_consistency(),  "Must be consistent at the end of birth case in remove");

        return Death();
    } else
    {
        // Death
        rLog(rlZigzagRemove,        "Death case in remove");

        ZIndex j                    = s->z_row.front();
        CRow c_row                  = s->c_row;

        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 
        // Compute reducers --- columns that we will be adding to other columns
        ZRow z_row                  = s->z_row;
        typedef typename ZRow::reverse_iterator             ZRowReverseIterator;
        typedef std::list<ZRowReverseIterator>              ReducersContainer;
        ReducersContainer  reducers;                        // list of ZColumns that we will be adding to other columns
        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", 
                                    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))); 
                low = (*(reducers.back()))->z_column.back();
                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", 
                                            (*(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", 
                                    reducers.size(), s->order);

        //show_all();

        // Add each reducer to the columns that follow them until the next reducer
        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, 
                         &ZNode::z_column, &SimplexNode::z_row,
                         &ZNode::b_row,    &BNode::b_column);
            if (cur != reducers.rbegin())
            {
                AssertMsg((*zcur)->z_column.back() == (**cur)->z_column.back(),
                          "The back of the z_columns must be the same.");
                (*zcur)->z_column.back()->low = *zcur;
            }
            else
                (**cur)->z_column.back()->low = z_list.end();

            zcur = *cur++;
            // This makes it inconsistent until the next iteration of this update loop
        }
        
        // Drop j and s
        Death d                     = visitor.death(*this, j);

        if (j->z_column.back()->low == j)
            j->z_column.back()->low = z_list.end();
        std::for_each(j->z_column.begin(), j->z_column.end(), make_remover(&SimplexNode::z_row, j));
        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()); 
        }
        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()");
        AssertMsg(s->c_row.empty(),     "c_row of s must be empty before erasing in remove()");
        visitor.erasing_z(*this, j);
        z_list.erase(j);
        s_list.erase(s);

        //show_all();

        AssertMsg(check_consistency(),  "Must be consistent when done in remove()");
        
        return d;
    }
}
        
template<class BID, class SD>
void
ZigzagPersistence<BID,SD>::
show_all()
{
    std::cout << "s_list:" << std::endl;
    for (SimplexIndex cur = s_list.begin(); cur != s_list.end(); ++cur)
    {
        std::cout << "  " << cur->order << ":" << std::endl;

        std::cout << "    z_row: ";
        for (typename ZRow::const_iterator zcur = cur->z_row.begin(); zcur != cur->z_row.end(); ++zcur)
            std::cout << (*zcur)->order << " ";
        std::cout << std::endl;

        std::cout << "    c_row: ";
        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;
        else
            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: ";
        for (typename ZColumn::const_iterator zcur = cur->z_column.begin(); zcur != cur->z_column.end(); ++zcur)
            std::cout << (*zcur)->order << " ";
        std::cout << std::endl;

        std::cout << "    b_row: ";
        for (typename BRow::const_iterator bcur = cur->b_row.begin(); bcur != cur->b_row.end(); ++bcur)
            std::cout << (*bcur)->order << " ";
        std::cout << std::endl;

        std::cout << "    low: ";
        if (cur->low != b_list.end()) 
            std::cout << cur->low->order;
        else
            std::cout << "none";
        std::cout << std::endl;
    }

    std::cout << "b_list:" << std::endl;
    for (BIndex cur = b_list.begin(); cur != b_list.end(); ++cur)
    {
        std::cout << "  " << cur->order << ":" << std::endl;

        std::cout << "    b_column: ";
        for (typename BColumn::const_iterator bcur = cur->b_column.begin(); bcur != cur->b_column.end(); ++bcur)
            std::cout << (*bcur)->order << " ";
        std::cout << std::endl;

        std::cout << "    c_column: ";
        for (typename CColumn::const_iterator ccur = cur->c_column.begin(); ccur != cur->c_column.end(); ++ccur)
            std::cout << (*ccur)->order << " ";
        std::cout << std::endl;
    }
}

template<class BID, class SD>
bool
ZigzagPersistence<BID,SD>::
check_consistency(SimplexIndex s_skip, ZIndex z_skip, BIndex b_skip)
{
#ifdef ZIGZAG_CONSISTENCY
    #warning "Checking consistency in ZigzagPersistence"

    Count(cZigzagConsistency);
    for (SimplexIndex cur = s_list.begin(); cur != s_list.end(); ++cur)
    {
        if (cur == s_skip) continue;
        //rLog(rlZigzagCheckConsistency,      "SimplexIndex cur: %d", cur->order);
        for (typename ZRow::const_iterator zcur = cur->z_row.begin(); zcur != cur->z_row.end(); ++zcur)
            if (std::find((*zcur)->z_column.begin(), (*zcur)->z_column.end(), cur) == (*zcur)->z_column.end())
            {
                rError("In check_consistency(): SimplexNode %d not found in z_column of %d", cur->order, (*zcur)->order);
                return false;
            }
        for (typename CRow::const_iterator ccur = cur->c_row.begin(); ccur != cur->c_row.end(); ++ccur)
            if (std::find((*ccur)->c_column.begin(), (*ccur)->c_column.end(), cur) == (*ccur)->c_column.end())
            {
                rError("In check_consistency(): SimplexNode %d not found in c_column of %d", cur->order, (*ccur)->order);
                return false;
            }
        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) 
        {
            rError("low of SimplexNode %d is incorrect", cur->order);
            return false;
        }
    }

    for (ZIndex cur = z_list.begin(); cur != z_list.end(); ++cur)
    {
        if (cur == z_skip) continue;

        //rLog(rlZigzagCheckConsistency,      "ZIndex cur: %d", cur->order);
        for (typename ZColumn::const_iterator scur = cur->z_column.begin(); scur != cur->z_column.end(); ++scur)
            if (std::find((*scur)->z_row.begin(), (*scur)->z_row.end(), cur) == (*scur)->z_row.end())
            {
                rError("In check_consistency(): ZNode %d not found in z_row of %d", cur->order, (*scur)->order);
                return false;
            }
        for (typename BRow::const_iterator bcur = cur->b_row.begin(); bcur != cur->b_row.end(); ++bcur)
            if (std::find((*bcur)->b_column.begin(), (*bcur)->b_column.end(), cur) == (*bcur)->b_column.end())
            {
                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) 
        {
            rError("low of ZNode %d is incorrect", cur->order);
            return false;
        }
        if (cur->z_column.back()->low != cur)
        {
            rError("The low of the back of the z_column must be set correctly");
            rError("  %d [%s], its back %d with low=%d", cur->order,
                                                         cur->z_column.tostring(out).c_str(),
                                                         cur->z_column.back()->order,
                                                         (cur->z_column.back()->low == z_list.end()) ? 0 : cur->z_column.back()->low->order);
            return false;
        }
    }

    for (BIndex cur = b_list.begin(); cur != b_list.end(); ++cur)
    {
        if (cur == b_skip) continue;

        //rLog(rlZigzagCheckConsistency,      "BIndex cur: %d", cur->order);
        for (typename BColumn::const_iterator zcur = cur->b_column.begin(); zcur != cur->b_column.end(); ++zcur)
            if (std::find((*zcur)->b_row.begin(), (*zcur)->b_row.end(), cur) == (*zcur)->b_row.end())
            {
                rError("In check_consistency(): BNode %d not found in b_row of %d", cur->order, (*zcur)->order);
                return false;
            }
        for (typename CColumn::const_iterator scur = cur->c_column.begin(); scur != cur->c_column.end(); ++scur)
            if (std::find((*scur)->c_row.begin(), (*scur)->c_row.end(), cur) == (*scur)->c_row.end())
            {
                rError("In check_consistency(): BNode %d not found in c_row of %d", cur->order, (*scur)->order);
                return false;
            }
        if (!(cur->b_column.empty() || cur->b_column.back()->low == cur))
        {
            rError("The low of the back of the b_column must be set correctly");
            return false;
        }

        // ZB == DC
        ZColumn zb, dc;
        std::for_each(cur->b_column.begin(), cur->b_column.end(), make_adder(&ZNode::z_column, zb));
        std::for_each(cur->c_column.begin(), cur->c_column.end(), make_adder(&SimplexNode::boundary, dc));
        zb.add(dc, cmp);
        if (!zb.empty())
        {
            rError("   b_column: [%s]",    cur->b_column.tostring(out).c_str());
            rError("   c_column: [%s]",    cur->c_column.tostring(out).c_str());
            rError("   zb - dc:  [%s]",    zb.tostring(out).c_str());
            rError("ZB = DC");
            return false;
        }
    }
#endif

    return true;
}

/* 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): 
                    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): 
                    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>
struct ZigzagPersistence<BID,SD>::Adder
{
                Adder(Member mm, Chain& cc):
                    m(mm), c(cc)                        {}

    template<class T>
    void        operator()(T& a)                        { c.add((*a).*m, cmp); }

    Member          m;
    Chain&          c;
    OrderComparison cmp;
};

        
// Function: add_chains()
//   
// 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>
void
ZigzagPersistence<BID,SD>::
add_chains(Index bg, Index end, IndexFrom j, PrimaryMember pm, SecondaryMember sm)
{
    add_chains(bg, end, j, pm, sm, pm);
}

// 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
// column member, and SecondaryMember is the corresponding row member.
template<class BID, class SD>
template<class IndexTo, class IndexFrom, class PrimaryMemberTo, class SecondaryMemberTo, class PrimaryMemberFrom>
void
ZigzagPersistence<BID,SD>::
add_chains(IndexTo bg, IndexTo end, IndexFrom j, PrimaryMemberTo pmt, SecondaryMemberTo smt, PrimaryMemberFrom pmf)
{
    for (IndexTo cur = bg; cur != end; ++cur)
        add_chain(*cur, j, pmt, smt, pmf);
}

// Function: change_basis()
//
// Changes basis by adding PrimaryMember pm of j to pm of every element in range [bg, end).
// In parallel it performs the reverse (complementary) update on the dual members, i.e.
// column and row operations are performed in sync, so that the product of the two matrices doesn't change
template<class BID, class SD>
template<class IndexTo, class IndexFrom, class PrimaryMember, class SecondaryMember, class DualPrimaryMember, class DualSecondaryMember>
void
ZigzagPersistence<BID,SD>::
change_basis(IndexTo bg, IndexTo end, IndexFrom j, PrimaryMember pm, SecondaryMember sm, DualPrimaryMember dpm, DualSecondaryMember dsm)
{
    for (IndexTo cur = bg; cur != end; ++cur)
    {
        add_chain(*cur, j,  pm,  sm,  pm);
        add_chain(j, *cur, dpm, dsm, dpm);
    }
}

template<class BID, class SD>
template<class Index, class PrimaryMember, class SecondaryMember>
void
ZigzagPersistence<BID,SD>::
add_chain(Index to, Index from, PrimaryMember pm, SecondaryMember sm)
{
    add_chain(to, from, pm, sm, pm);
}

// Function: add_chain()
//
// 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>
void
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, 
                            ((*from).*pmf).tostring(out).c_str(),
                            (*to).order,
                            ((*to).*pmt).tostring(out).c_str());

    // Fix secondaries
    std::for_each(make_intersection_iterator(((*from).*pmf).begin(),  ((*from).*pmf).end(),
                                               ((*to).*pmt).begin(),    ((*to).*pmt).end(),
                                             cmp),
                  make_intersection_iterator(((*from).*pmf).end(),    ((*from).*pmf).end(),
                                               ((*to).*pmt).end(),      ((*to).*pmt).end(),
                                             cmp),
                  make_remover(smt, to));
    std::for_each(make_difference_iterator(((*from).*pmf).begin(),    ((*from).*pmf).end(),
                                             ((*to).*pmt).begin(),      ((*to).*pmt).end(),
                                           cmp),
                  make_difference_iterator(((*from).*pmf).end(),      ((*from).*pmf).end(),
                                             ((*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());
}


/* ZigzagVisitor */
template<class BID, class SD>
typename ZigzagPersistence<BID,SD>::SimplexIndex
ZigzagPersistence<BID,SD>::ZigzagVisitor::
new_simplex(ZigzagPersistence& zz)
{
    unsigned order      = zz.s_list.empty() ? 0 : boost::prior(zz.s_list.end())->order + 1;
    zz.s_list.push_back(SimplexNode(order, zz.z_list.end()));
    return boost::prior(zz.s_list.end());
}

template<class BID, class SD>
typename ZigzagPersistence<BID,SD>::ZIndex
ZigzagPersistence<BID,SD>::ZigzagVisitor::
new_z_in_add(ZigzagPersistence& zz, const ZColumn& z, const BRow& u)
{
    int order                   = zz.z_list.empty() ? 0 : boost::prior(zz.z_list.end())->order + 1;
    zz.z_list.push_back(ZNode(order, zz.b_list.end()));
    return boost::prior(zz.z_list.end());
}

template<class BID, class SD>
typename ZigzagPersistence<BID,SD>::BIndex
ZigzagPersistence<BID,SD>::ZigzagVisitor::
select_j_in_remove(ZigzagPersistence& zz, const CRow& c_row)
{
    return c_row.front();
}

template<class BID, class SD>
typename ZigzagPersistence<BID,SD>::ZIndex
ZigzagPersistence<BID,SD>::ZigzagVisitor::
new_z_in_remove(ZigzagPersistence& zz)
{
    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();
}

template<class BID, class SD>
typename ZigzagPersistence<BID,SD>::Death
ZigzagPersistence<BID,SD>::ZigzagVisitor::
death(ZigzagPersistence& zz, ZIndex dying_z)
{
    return Death(dying_z->birth);
}