author Dmitriy Morozov <>
Thu, 18 Dec 2008 16:43:42 -0800
changeset 97 0a9bd3f34419
parent 90 include/topology/filtration.hpp@dea0e9726c62
child 99 6c0da7931e4d
permissions -rw-r--r--
Intermediate commit while converting to the new architecture * Converted Filtration into StaticPersistence and DynamicPersistenceTrails (both right now work with the underlying std::vector rather than std::list order) * Filtration is now just an auxilliary glue (a map between Complex and Persistence) * Whether chains are vectors or lists can be interchanged * Added PersistenceDiagram with a simple bottleneck_distance() function * Converted triangle, alphashapes3d, and cech-complex examples * Lots of organizational changes (factoring utilities out into containers.h, indirect.h, property-maps.h) * Trying to document along the way with NaturalDocs-type comments

#include <utilities/log.h>
#include <utilities/containers.h>
#include <boost/utility/enable_if.hpp>
#include <boost/utility.hpp>
#include <utilities/property-maps.h>

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

static Counter*  cPersistencePair =                         GetCounter("persistence/pair");
static Counter*  cPersistencePairBoundaries =               GetCounter("persistence/pair/boundaries");
static Counter*  cPersistencePairCycleLength =              GetCounter("persistence/pair/cyclelength");
#endif // COUNTERS

template<class D, class OD>
template<class Filtration>
StaticPersistence<D, OD>::
StaticPersistence(const Filtration& filtration, 
                  const OrderComparison& ocmp):
    OrderIndex                          ocur = begin();
    OffsetMap<size_t, OrderIndex>       om(0, ocur);            // TODO: this is customized for std::vector Order
    for (typename Filtration::Index cur = filtration.begin(); cur != filtration.end(); ++cur, ++ocur)
        // Convert the Filtration::IndexBoundary into a Cycle, and 
        typedef typename OrderDescriptor::Chains::Cycle         Cycle;

        OrderElement& oe = *ocur;
        oe.pair = ocur;

        filtration.boundary(cur, oe.cycle, om);

template<class D, class OD>
template<class Visitor>
StaticPersistence<D, OD>::
pair_simplices(OrderIndex bg, OrderIndex end, const Visitor& visitor)
    typename Traits::OutputMap outmap(order_);

    // FIXME: need sane output for logging
    rLog(rlPersistence, "Entered: pair_simplices");
    for (OrderIndex j = bg; j != end; ++j)
        rLog(rlPersistence, "Simplex %s", outmap(j).c_str());

        OrderElement& oe = *j;
        typedef typename OrderDescriptor::Chains::Cycle         Cycle;
        Cycle& z = oe.cycle;
        rLog(rlPersistence, "  has boundary: %s", z.tostring(outmap).c_str());
        CountNum(cPersistencePairBoundaries, oe.cycle.size());

            OrderIndex i =;            // take the youngest element with respect to the OrderComparison
            rLog(rlPersistence, "  %s: %s", outmap(i).c_str(), outmap(i->pair).c_str());
            // TODO: is this even a meaningful assert?
            AssertMsg(!ocmp_(i, j), 
                      "Simplices in the cycle must precede current simplex: (%s in cycle of %s)",
                      outmap(i).c_str(), outmap(j).c_str());

            // i is not paired, so we pair j with i
            if (i->pair == i)
                rLog(rlPersistence, "  Pairing %s and %s with cycle %s", 
                                   outmap(i).c_str(), outmap(j).c_str(), 
                i->pair = j;
                j->pair = i;
                CountNum(cPersistencePairCycleLength, z.size());
                CountBy(cPersistencePairCycleLength, z.size());

            // update element
            z.add(i->pair->cycle, ocmp_);
            visitor.update(j, i);
            rLog(rlPersistence, "    new cycle: %s", z.tostring(outmap).c_str());
        rLog(rlPersistence, "Finished with %s: %s", 
                            outmap(j).c_str(), outmap(oe.pair).c_str());

/* Private */