author Christos Mantoulidis <cmad@stanford.edu>
Tue, 04 Aug 2009 13:23:16 -0700
changeset 156 f75fb57d2831
parent 153 7731c42892de
child 182 451748b3c888
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 <topology/weighted-rips.h>
#include <topology/filtration.h>
#include <topology/static-persistence.h>
#include <topology/dynamic-persistence.h>
#include <topology/persistence-diagram.h>

#include <geometry/weighted-l2distance.h>
#include <geometry/weighted-cechdistance.h>

#include <geometry/distances.h>

#include <utilities/containers.h>           // for BackInsertFunctor
#include <utilities/timer.h>

#include <vector>

#include <boost/program_options.hpp>

typedef         PairwiseDistances<PointContainer, WeightedL2Distance>   PairDistances;
typedef         PairDistances::DistanceType                             DistanceType;
typedef         PairDistances::IndexType                                Vertex;

typedef         WeightedRips<PairDistances>                             Generator;
typedef         Generator::Simplex                                      Smplx;
typedef         std::vector<Smplx>                                      SimplexVector;
typedef         Filtration<SimplexVector, unsigned>                     Fltr;
typedef         StaticPersistence<>                                     Persistence;
//typedef         DynamicPersistenceChains<>                              Persistence;
typedef         PersistenceDiagram<>                                    PDgm;

void            program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance);

int main(int argc, char* argv[])
#ifdef LOGGING
    rlog::RLogInit(argc, argv);

    stderrLog.subscribeTo( RLOG_CHANNEL("error") );

    Dimension               skeleton;
    DistanceType            max_distance;
    std::string             infilename;

    program_options(argc, argv, infilename, skeleton, max_distance);

    PointContainer          points;
    read_weighted_points(infilename, points);

    PairDistances           distances(points);
    Generator               rips(distances);
    Generator::Evaluator    size(distances);
    Generator::Comparison   cmp (distances);
    SimplexVector           complex;
    // Generate skeleton of the weighted Rips complex for epsilon = 50
    rips.generate(skeleton, max_distance, make_push_back_functor(complex));

    std::sort(complex.begin(), complex.end(), Smplx::VertexComparison());       // unnecessary
    std::cout << "# Generated complex of size: " << complex.size() << std::endl;

    // Generate filtration with respect to distance and compute its persistence
    Fltr f(complex.begin(), complex.end(), cmp);

    Timer persistence_timer; persistence_timer.start();
    Persistence p(f);

    // Output cycles
    for (Persistence::OrderIndex cur = p.begin(); cur != p.end(); ++cur)
        Persistence::Cycle& cycle = cur->cycle;

        if (!cur->sign())        // only negative simplices have non-empty cycles
            Persistence::OrderIndex birth = cur->pair;      // the cycle that cur killed was born when we added birth (another simplex)

            const Smplx& b = f.simplex(f.begin() + (birth - p.begin()));        // eventually this will be encapsulated behind an interface
            const Smplx& d = f.simplex(f.begin() + (cur   - p.begin()));
            if (b.dimension() >= skeleton) continue;
            std::cout << b.dimension() << " " << size(b) << " " << size(d) << std::endl;
        } else if (cur->pair == cur)    // positive could be unpaired
            const Smplx& b = f.simplex(f.begin() + (cur - p.begin()));
            if (b.dimension() >= skeleton) continue;
            std::cout << b.dimension() << " " << size(b) << " inf" << std::endl;
            //cycle = cur->chain;
    persistence_timer.check("# Persistence timer");

void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance)
    namespace po = boost::program_options;

    po::options_description     hidden("Hidden options");
        ("input-file",          po::value<std::string>(&infilename),        "Point set whose weighed Rips zigzag we want to compute");
    po::options_description visible("Allowed options", 100);
        ("help,h",                                                                                  "produce help message")
        ("skeleton-dimsnion,s", po::value<Dimension>(&skeleton)->default_value(2),                  "Dimension of the weighted Rips complex we want to compute")
        ("max-distance,m",      po::value<DistanceType>(&max_distance)->default_value(Infinity),    "Maximum value for the weighted Rips complex construction");
    std::vector<std::string>    log_channels;
        ("log,l",               po::value< std::vector<std::string> >(&log_channels),           "log channels to turn on (info, debug, etc)");

    po::positional_options_description pos;
    pos.add("input-file", 1);
    po::options_description all; all.add(visible).add(hidden);

    po::variables_map vm;
    po::store(po::command_line_parser(argc, argv).
                  options(all).positional(pos).run(), vm);

    for (std::vector<std::string>::const_iterator cur = log_channels.begin(); cur != log_channels.end(); ++cur)
        stderrLog.subscribeTo( RLOG_CHANNEL(cur->c_str()) );
     * Interesting channels
     * "info", "debug", "topology/persistence"

    if (vm.count("help") || !vm.count("input-file"))
        std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl;
        std::cout << visible << std::endl;