include/topology/cohomology-persistence.hpp
author Christos Mantoulidis <cmad@stanford.edu>
Tue, 04 Aug 2009 13:23:16 -0700
branchdev
changeset 156 f75fb57d2831
parent 137 069596c71902
child 172 a6605dc232f2
permissions -rw-r--r--
Changed implementation of WeightedRips to store simplex values (max distance between simplices' vertices) as an invisible layer on top of each simplex object, so that the data() field of WeightedRips has been freed for use by the users again.

#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);
}