include/topology/dynamic-persistence.hpp
author Dmitriy Morozov <dmitriy@mrzv.org>
Wed, 25 Mar 2009 09:56:52 -0700
branchdev
changeset 119 505b3795d239
parent 115 a3410b6ba79c
child 125 0a2c2283e4a8
permissions -rw-r--r--
Added ZigzagPersistence::is_alive() + consistency zigzag outputs betti 1 at each stage

#include <utilities/log.h>

#ifdef LOGGING
static rlog::RLogChannel* rlTranspositions =    DEF_CHANNEL("topology/persistence/transpositions", rlog::Log_Debug);
#endif // LOGGING

#ifdef COUNTERS
static Counter*  cTransposition =               GetCounter("persistence/transposition");
static Counter*  cTranspositionDiffDim =        GetCounter("persistence/transposition/diffdim");
static Counter*  cTranspositionCase12 =         GetCounter("persistence/transposition/case/1/2");
static Counter*  cTranspositionCase12s =        GetCounter("persistence/transposition/case/1/2/special");
static Counter*  cTranspositionCase112 =        GetCounter("persistence/transposition/case/1/1/2");
static Counter*  cTranspositionCase111 =        GetCounter("persistence/transposition/case/1/1/1");
static Counter*  cTranspositionCase22 =         GetCounter("persistence/transposition/case/2/2");
static Counter*  cTranspositionCase212 =        GetCounter("persistence/transposition/case/2/1/2");
static Counter*  cTranspositionCase211 =        GetCounter("persistence/transposition/case/2/1/1");
static Counter*  cTranspositionCase32 =         GetCounter("persistence/transposition/case/3/2");
static Counter*  cTranspositionCase31 =         GetCounter("persistence/transposition/case/3/1");
static Counter*  cTranspositionCase4 =          GetCounter("persistence/transposition/case/4");
#endif // COUNTERS


/* Trails */

template<class D, class OD, class CI, class CC>
template<class Filtration>
DynamicPersistenceTrails<D,OD,CI,CC>::
DynamicPersistenceTrails(const Filtration& f, const OrderComparison& ocmp, const CC& ccmp):
    Parent(f, ocmp), ccmp_(ccmp)
{}
        
template<class D, class OD, class CI, class CC>
void
DynamicPersistenceTrails<D,OD,CI,CC>::
pair_simplices()
{ 
    Parent::pair_simplices(begin(), end(), PairingTrailsVisitor(begin(), ccmp_));
}

template<class D, class OD, class CI, class CC>
template<class Visitor>
bool
DynamicPersistenceTrails<D,OD,CI,CC>::
transpose(OrderIndex i, const Visitor& visitor)
{
#if LOGGING
    typename Traits::OutputMap outmap(order());
#endif

    Count(cTransposition);
    typedef                 OrderIndex                                  Index;
    typedef                 typename TrailData::Trail::iterator         TrailIterator;

    visitor.transpose(i);
    
    Index i_prev = i++;

#if 0       // Persistence no longer has the notion of dimension
    if (i_prev->dimension() != i->dimension())
    {
        swap(i_prev, i);
        rLog(rlTranspositions, "Different dimension");
        Count(cTranspositionDiffDim);
        return false;
    }
#endif
    
    bool si = i_prev->sign(), sii = i->sign();
    if (si && sii)
    {
        rLog(rlTranspositions, "Trail prev: %s", i_prev->trail.tostring(outmap).c_str());

        // Case 1
        boost::optional<TrailIterator> i_in_i_prev = i_prev->trail.contains(i);
        if (i_in_i_prev)
        {
            rLog(rlTranspositions, "Case 1, U[i,i+1] = 1");
            i_prev->trail.erase(*i_in_i_prev);
        }

        Index k = i_prev->pair;
        Index l = i->pair;

        // Explicit treatment of unpaired simplex
        if (l == i)
        {
            swap(i_prev, i);
            rLog(rlTranspositions, "Case 1.2 --- unpaired");
            rLog(rlTranspositions, "%s", outmap(i_prev).c_str());
            Count(cTranspositionCase12);
            return false;
        } else if (k == i_prev)
        {
            if (!(l->cycle.contains(i_prev)))
            {
                // Case 1.2
                swap(i_prev, i);
                rLog(rlTranspositions, "Case 1.2 --- unpaired");
                rLog(rlTranspositions, outmap(i_prev).c_str());
                Count(cTranspositionCase12);
                return false;
            } else
            {
                // Case 1.2 --- special version (plain swap, but pairing switches)
                swap(i_prev, i);
                pairing_switch(i_prev, i);
                visitor.switched(i_prev, Case12);
                rLog(rlTranspositions, "Case 1.2 --- unpaired (pairing switch)");
                rLog(rlTranspositions, outmap(i_prev).c_str());
                Count(cTranspositionCase12s);
                return true;
            }
        }
        
        rLog(rlTranspositions, "l cycle: %s", l->cycle.tostring(outmap).c_str());
        if (!(l->cycle.contains(i_prev)))
        {
            // Case 1.2
            swap(i_prev, i);
            rLog(rlTranspositions, "Case 1.2");
            Count(cTranspositionCase12);
            return false;
        } else
        {
            // Case 1.1
            if (not2(ccmp_)(k,l))
            {
                // Case 1.1.1
                swap(i_prev, i);
                l->cycle.add(k->cycle, ccmp_);        // Add column k to l
                k->trail.add(l->trail, ccmp_);        // Add row l to k
                rLog(rlTranspositions, "Case 1.1.1");
                Count(cTranspositionCase111);
                return false;
            } else
            {
                // Case 1.1.2
                swap(i_prev, i);
                k->cycle.add(l->cycle, ccmp_);        // Add column l to k
                l->trail.add(k->trail, ccmp_);        // Add row k to l
                pairing_switch(i_prev, i);
                visitor.switched(i_prev, Case112);
                rLog(rlTranspositions, "Case 1.1.2");
                Count(cTranspositionCase112);
                return true;
            }
        }
    } else if (!si && !sii)
    {
        // Case 2
        if (!(i_prev->trail.contains(i)))
        {
            // Case 2.2
            swap(i_prev, i);
            rLog(rlTranspositions, "Case 2.2");
            Count(cTranspositionCase22);
            return false;
        } else
        {
            // Case 2.1
            Index low_i = i_prev->pair;
            Index low_ii = i->pair;
            i_prev->trail.add(i->trail, ccmp_);            // Add row i to i_prev
            i->cycle.add(i_prev->cycle, ccmp_);            // Add column i_prev to i
            swap(i_prev, i);    
            if (not2(ccmp_)(low_ii, low_i))
            {
                // Case 2.1.2
                i_prev->cycle.add(i->cycle, ccmp_);        // Add column i to i_prev (after transposition)
                i->trail.add(i_prev->trail, ccmp_);        // Add row i to i_prev
                pairing_switch(i_prev, i);
                visitor.switched(i_prev, Case212);
                rLog(rlTranspositions, "Case 2.1.2");
                Count(cTranspositionCase212);
                return true;
            } 
            
            // Case 2.1.1
            rLog(rlTranspositions, "Case 2.1.1");
            Count(cTranspositionCase211);
            return false;
        }
    } else if (!si && sii)
    {
        // Case 3
        if (!(i_prev->trail.contains(i)))
        {
            // Case 3.2
            swap(i_prev, i);
            rLog(rlTranspositions, "Case 3.2");
            Count(cTranspositionCase32);
            return false;
        } else
        {
            // Case 3.1
            i_prev->trail.add(i->trail, ccmp_);            // Add row i to i_prev
            i->cycle.add(i_prev->cycle, ccmp_);            // Add column i_prev to i
            swap(i_prev, i);
            i_prev->cycle.add(i->cycle, ccmp_);            // Add column i_prev to i (after transposition)
            i->trail.add(i_prev->trail, ccmp_);            // Add row i to i_prev
            pairing_switch(i_prev, i);
            visitor.switched(i_prev, Case31);
            rLog(rlTranspositions, "Case 3.1");
            Count(cTranspositionCase31);
            return true;
        }
    } else if (si && !sii)
    {
        // Case 4
        boost::optional<TrailIterator> i_in_i_prev = i_prev->trail.contains(i);
        if (i_in_i_prev)
        {
            rLog(rlTranspositions, "Case 4, U[i,i+1] = 1");
            i_prev->trail.erase(*i_in_i_prev);
        }
        swap(i_prev, i);
        rLog(rlTranspositions, "Case 4");
        Count(cTranspositionCase4);
        return false;
    }
    
    return false; // to avoid compiler complaints; we should never reach this point
}

template<class D, class OD, class CI, class CC>
void
DynamicPersistenceTrails<D,OD,CI,CC>::
swap(OrderIndex i, OrderIndex j)
{
    std::swap<Data>(*i, *j);
    std::swap(i->pair, j->pair);

    std::swap(i->cycle, j->cycle);          // TODO: double-check that the STL container specializations actually get invoked
    std::swap(i->trail, j->trail);
}

template<class D, class OD, class CI, class CC>
void
DynamicPersistenceTrails<D,OD,CI,CC>::
pairing_switch(OrderIndex i, OrderIndex j)
{
    OrderIndex i_pair = i->pair;
    OrderIndex j_pair = j->pair;

    if (i_pair == i)
        j->pair = j;
    else
    {
        j->pair = i_pair;
        i_pair->pair = j;
    }

    if (j_pair == j)
        i->pair = i;
    else
    {
        i->pair = j_pair;
        j_pair->pair = i;
    }
}


/* Chains */

template<class D, class OD, class CI, class CC>
template<class Filtration>
DynamicPersistenceChains<D,OD,CI,CC>::
DynamicPersistenceChains(const Filtration& f, const OrderComparison& ocmp, const CC& ccmp):
    Parent(f, ocmp), ccmp_(ccmp)
{}
        
template<class D, class OD, class CI, class CC>
void
DynamicPersistenceChains<D,OD,CI,CC>::
pair_simplices()
{ 
    Parent::pair_simplices(begin(), end(), PairingChainsVisitor(begin(), ccmp_, size()));
}