author Christos Mantoulidis <cmad@stanford.edu>
Tue, 04 Aug 2009 13:23:16 -0700
changeset 156 f75fb57d2831
parent 155 1dde8a05fcd7
child 157 700cbac5b23c
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/cohomology-persistence.h>
#include <topology/weighted-rips.h>

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

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

#include <string>

#include <boost/tuple/tuple.hpp>
#include <boost/program_options.hpp>
#include <boost/progress.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     SimplexVector::const_iterator                           SV_const_iterator;

typedef     boost::tuple<Dimension, DistanceType>                   BirthInfo;
typedef     CohomologyPersistence<BirthInfo, SV_const_iterator>     Persistence;
typedef     Persistence::SimplexIndex                               Index;
typedef     Persistence::Death                                      Death;
typedef     std::map<Smplx, Index, Smplx::VertexComparison>         Complex;

#include "output.h"         // for output_*()

void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime, std::string& boundary_name, std::string& cocycle_prefix, std::string& vertices_name, std::string& diagram_name);

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

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

    Dimension               skeleton;
    DistanceType            max_distance;
    ZpField::Element        prime;
    std::string             infilename, boundary_name, cocycle_prefix, vertices_name, diagram_name;

    program_options(argc, argv, infilename, skeleton, max_distance, prime, boundary_name, cocycle_prefix, vertices_name, diagram_name);
    std::ofstream           bdry_out(boundary_name.c_str());
    std::ofstream           vertices_out(vertices_name.c_str());
    std::ofstream           diagram_out(diagram_name.c_str());
    std::cout << "Boundary matrix: " << boundary_name << std::endl;
    std::cout << "Cocycles:        " << cocycle_prefix << "*.ccl" << std::endl;
    std::cout << "Vertices:        " << vertices_name << std::endl;
    std::cout << "Diagram:         " << diagram_name << std::endl;

    Timer total_timer; total_timer.start();
    PointContainer          points;
    read_weighted_points(infilename, points);

    PairDistances           distances(points);
    Generator               rips(distances);
    Generator::Evaluator    size(distances);
    Generator::Comparison   cmp(distances);
    SimplexVector           v;
    Complex                 c;
    Timer rips_timer; rips_timer.start();
    rips.generate(skeleton, max_distance, make_push_back_functor(v));
    std::sort(v.begin(), v.end(), cmp);
    std::cout << "Simplex vector generated, size: " << v.size() << std::endl;

    /*output_boundary_matrix(bdry_out, v, cmp);
    output_vertex_indices(vertices_out, v);*/

    Timer persistence_timer; persistence_timer.start();
    ZpField                 zp(prime);
    Persistence             p(zp);
    boost::progress_display show_progress(v.size());
    for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
        std::vector<Index>      boundary;
        for (Smplx::BoundaryIterator bcur  = cur->boundary_begin(); bcur != cur->boundary_end(); ++bcur)
        Index idx; Death d;
        bool store = cur->dimension() < skeleton;
        boost::tie(idx, d)      = p.add(boundary.begin(), boundary.end(), boost::make_tuple(cur->dimension(), size(*cur)), store, cur);
        // c[*cur] = idx;
        if (store)
            c.insert(std::make_pair(*cur, idx));

        if (d && (size(*cur) - d->get<1>()) > 0)
            AssertMsg(d->get<0>() == cur->dimension() - 1, "Dimensions must match");
            diagram_out << (cur->dimension() - 1) << " " << d->get<1>() << " " << size(*cur) << std::endl;
    // output infinte persistence pairs 
    for (Persistence::CocycleIndex cur = p.begin(); cur != p.end(); ++cur)
        diagram_out << cur->birth.get<0>() << " " << cur->birth.get<1>() << " inf" << std::endl;

    // p.show_cocycles();
    // Output alive cocycles of dimension 1
    unsigned i = 0;
    for (Persistence::Cocycles::const_iterator cur = p.begin(); cur != p.end(); ++cur)
        if (cur->birth.get<0>() != 1) continue;
        output_cocycle(cocycle_prefix, i, v, *cur, prime, cmp);
        // std::cout << "Cocycle of dimension: " << cur->birth.get<0>() << " born at " << cur->birth.get<1>() << std::endl;
    rips_timer.check("Rips timer");
    persistence_timer.check("Persistence timer");
    total_timer.check("Total timer");

void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime, std::string& boundary_name, std::string& cocycle_prefix, std::string& vertices_name, std::string& diagram_name)
    namespace po = boost::program_options;

    po::options_description     hidden("Hidden options");
        ("input-file",          po::value<std::string>(&infilename),        "Point set whose 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 Rips complex we want to compute")
        ("prime,p",             po::value<ZpField::Element>(&prime)->default_value(11),             "Prime p for the field F_p")
        ("max-distance,m",      po::value<DistanceType>(&max_distance)->default_value(Infinity),    "Maximum value for the Rips complex construction")
        ("boundary,b",          po::value<std::string>(&boundary_name),                             "Filename where to output the boundary matrix")
        ("cocycle,c",           po::value<std::string>(&cocycle_prefix),                            "Prefix of the filename where to output the 1-dimensional cocycles")
        ("vertices,v",          po::value<std::string>(&vertices_name),                             "Filename where to output the simplex-vertex mapping")
        ("diagram,d",           po::value<std::string>(&diagram_name),                              "Filename where to output the persistence diagram");
    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()) );

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