include/topology/cohomology-persistence.hpp
author mridul@ramrod.Stanford.EDU
Mon, 31 Aug 2009 13:45:28 -0700
branchdev
changeset 168 3b0665eb35bc
parent 137 069596c71902
child 172 a6605dc232f2
permissions -rw-r--r--
The if-else cases in lsqr.py were not exhaustive. They become exhaustive by changing ab > aa to ab >= aa.

#include <boost/utility.hpp>
#include <queue>
#include <vector>

#include <utilities/log.h>
#include <utilities/indirect.h>

#ifdef LOGGING
static rlog::RLogChannel* rlCohomology =                DEF_CHANNEL("topology/cohomology",        rlog::Log_Debug);
#endif

template<class BirthInfo, class SimplexData, class Field>
class CohomologyPersistence<BirthInfo, SimplexData, Field>::CompareSNode
{
    public:
        bool        operator()(const SNode& s1, const SNode& s2) const                  { return s1.si->order < s2.si->order; }
};

template<class BirthInfo, class SimplexData, class Field>
template<class BI>
typename CohomologyPersistence<BirthInfo, SimplexData, Field>::IndexDeathPair
CohomologyPersistence<BirthInfo, SimplexData, Field>::
add(BI begin, BI end, BirthInfo birth, bool store, const SimplexData& sd)
{
    // Create simplex representation
    simplices_.push_back(SHead(sd, simplices_.empty() ? 0 : (simplices_.back().order + 1)));
    SimplexIndex    si = boost::prior(simplices_.end());

    // Find out if there are cocycles that evaluate to non-zero on the new simplex
    typedef         std::list<CocycleCoefficientPair>           Candidates;
    Candidates      candidates, candidates_bulk;
    rLog(rlCohomology, "Boundary");
    bool sign = true;       // TODO: this is very crude, fix later
    for (BI cur = begin; cur != end; ++cur)
    {
        rLog(rlCohomology, "  %d %d", (*cur)->order, sign);
        for (typename ZRow::const_iterator zcur = (*cur)->row.begin(); zcur != (*cur)->row.end(); ++zcur)
            candidates_bulk.push_back(std::make_pair(zcur->ci, 
                                                     (sign ? zcur->coefficient : field_.neg(zcur->coefficient))));
        sign = !sign;
    }

    candidates_bulk.sort(make_first_comparison(make_indirect_comparison(std::less<Cocycle>())));
    
#if LOGGING    
    rLog(rlCohomology,  "  Candidates bulk");
    for (typename Candidates::iterator cur  = candidates_bulk.begin(); 
                                       cur != candidates_bulk.end(); ++cur)
        rLog(rlCohomology, "    %d %d", cur->first->order, cur->second);
#endif

    // Remove duplicates
    {
        typename Candidates::const_iterator cur = candidates_bulk.begin();
        while (cur != candidates_bulk.end())
        {
            typename Candidates::const_iterator next = cur;
            FieldElement sum = field_.zero();
            while (next != candidates_bulk.end() && next->first == cur->first) 
            { 
                sum = field_.add(sum, next->second);
                ++next; 
            }
    
            rLog(rlCohomology, "  Bulk candidate %d sum %d", cur->first->order, sum);
            if (!field_.is_zero(sum))
                candidates.push_back(CocycleCoefficientPair(cur->first, sum));
    
            cur = next;
        }
    }

    // Birth
    if (candidates.empty())
    {
        if (!store)
        {
            simplices_.pop_back();
            return std::make_pair(simplices_.begin(), Death());         // TODO: shouldn't return front
        }
        
        unsigned order = cocycles_.empty() ? 0 : cocycles_.front().order + 1;
        cocycles_.push_front(Cocycle(birth, order));
        
        rLog(rlCohomology,  "Birth: %d", cocycles_.front().order);

        // set up the cocycle
        ZColumn& cocycle = cocycles_.front().zcolumn;
        cocycle.push_back(SNode(si, field_.id(), cocycles_.begin()));
        si->row.push_back(cocycles_.front().zcolumn.front());
        rLog(rlCohomology,  "  Cocyle: %d", si->order);

        return std::make_pair(si, Death());
    }

    // Death
    rLog(rlCohomology,  "Death");

#if LOGGING
    // Debug only, output candidates
    rLog(rlCohomology,  "  Candidates");
    for (typename Candidates::iterator cur  = candidates.begin(); cur != candidates.end(); ++cur)
        rLog(rlCohomology, "    %d %d", cur->first->order, cur->second);
#endif

    CocycleCoefficientPair& z   = candidates.front();
    Death d                     = z.first->birth;

    // add z to everything else in candidates
    for (typename Candidates::iterator cur  = boost::next(candidates.begin()); 
                                       cur != candidates.end(); ++cur)
        add_cocycle(*cur, z);

    for (typename ZColumn::iterator cur = z.first->zcolumn.begin(); cur != z.first->zcolumn.end(); ++cur)
        cur->unlink();
    
    cocycles_.erase(candidates.front().first);

    return std::make_pair(si, d);
}
        
template<class BirthInfo, class SimplexData, class Field>
void
CohomologyPersistence<BirthInfo, SimplexData, Field>::
show_cocycles() const
{
    std::cout << "Cocycles: " << cocycles_.size() << std::endl;
    for (typename Cocycles::const_iterator cur = cocycles_.begin(); cur != cocycles_.end(); ++cur)
    {
        // std::cout << cur->order << " (" << cur->birth << "): ";
        for (typename ZColumn::const_iterator zcur = cur->zcolumn.begin(); zcur != cur->zcolumn.end(); ++zcur)
            std::cout << zcur->coefficient << " * " << zcur->si->order << ", ";
        std::cout << std::endl;
    }
}

template<class BirthInfo, class SimplexData, class Field>
void
CohomologyPersistence<BirthInfo, SimplexData, Field>::
add_cocycle(CocycleCoefficientPair& to, CocycleCoefficientPair& from)
{
    rLog(rlCohomology,  "Adding cocycle %d to %d", from.first->order, to.first->order);

    ZColumn         nw;
    FieldElement    multiplier = field_.neg(field_.div(to.second, from.second));
    CocycleIndex    ci = to.first;

    typename ZColumn::iterator tcur = to.first->zcolumn.begin();
    typename ZColumn::iterator fcur = from.first->zcolumn.begin();
    CompareSNode cmp;
    while (tcur != to.first->zcolumn.end() && fcur != from.first->zcolumn.end())
    {
        rLog(rlCohomology, "  %d %d", tcur->si->order, fcur->si->order);
        if (cmp(*tcur, *fcur))
        {
            nw.push_back(*tcur);
            ++tcur;
        }
        else if (cmp(*fcur, *tcur))
        {
            nw.push_back(SNode(fcur->si, field_.mul(multiplier, fcur->coefficient), ci));
            ++fcur;
        }
        else        // equality
        {
            FieldElement res = field_.mul(multiplier, tcur->coefficient);
            res = field_.add(fcur->coefficient, res);
            if (!field_.is_zero(res))
                nw.push_back(SNode(fcur->si, res, ci));
            ++tcur; ++fcur;
        }
    }
    for (; tcur != to.first->zcolumn.end(); ++tcur)
    {
        rLog(rlCohomology, "  %d", tcur->si->order);
        nw.push_back(SNode(*tcur));
    }
    for (; fcur != from.first->zcolumn.end(); ++fcur)
    {
        rLog(rlCohomology, "  %d", fcur->si->order);
        nw.push_back(SNode(fcur->si, field_.mul(multiplier, fcur->coefficient), ci));
    }


    for (typename ZColumn::iterator cur = to.first->zcolumn.begin(); cur != to.first->zcolumn.end(); ++cur)
        cur->unlink();

    to.first->zcolumn.swap(nw);

    for (typename ZColumn::iterator cur = to.first->zcolumn.begin(); cur != to.first->zcolumn.end(); ++cur)
        cur->si->row.push_back(*cur);
}