Added consistency zigzag dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Tue, 24 Mar 2009 15:17:39 -0700
branchdev
changeset 118 c4e25fb4082c
parent 117 06a4361bddaa
child 119 505b3795d239
Added consistency zigzag
examples/CMakeLists.txt
examples/consistency/CMakeLists.txt
examples/consistency/rips-consistency-zigzag.cpp
examples/rips/l2distance.h
examples/rips/rips-zigzag.cpp
include/topology/rips.h
include/topology/rips.hpp
--- a/examples/CMakeLists.txt	Mon Mar 23 11:40:29 2009 -0700
+++ b/examples/CMakeLists.txt	Tue Mar 24 15:17:39 2009 -0700
@@ -1,6 +1,7 @@
 add_subdirectory			(alphashapes)
 #add_subdirectory			(ar-vineyard)
 add_subdirectory			(cech-complex)
+add_subdirectory			(consistency)
 add_subdirectory			(cohomology)
 add_subdirectory			(fitness)
 #add_subdirectory			(grid)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/consistency/CMakeLists.txt	Tue Mar 24 15:17:39 2009 -0700
@@ -0,0 +1,7 @@
+set                         (targets                        
+                             rips-consistency-zigzag)
+                             
+foreach                     (t ${targets})
+    add_executable          (${t} ${t}.cpp)
+    target_link_libraries   (${t} ${libraries} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${Boost_SERIALIZATION_LIBRARY})
+endforeach                  (t ${targets})
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/consistency/rips-consistency-zigzag.cpp	Tue Mar 24 15:17:39 2009 -0700
@@ -0,0 +1,338 @@
+#include <topology/rips.h>
+#include <topology/zigzag-persistence.h>
+#include <utilities/types.h>
+#include <utilities/containers.h>
+
+#include <utilities/log.h>
+#include <utilities/memory.h>       // for report_memory()
+#include <utilities/timer.h>
+
+#include <map>
+#include <cmath>
+#include <fstream>
+#include <sstream>
+#include <stack>
+#include <cstdlib>
+
+#include <boost/tuple/tuple.hpp>
+#include <boost/program_options.hpp>
+#include <boost/progress.hpp>
+
+#include "../rips/l2distance.h"         // Point, PointContainer, L2DistanceType
+
+#ifdef COUNTERS
+static Counter*  cComplexSize =                     GetCounter("rips/size");
+static Counter*  cOperations =                      GetCounter("rips/operations");
+#endif // COUNTERS
+
+typedef     PairwiseDistances<PointContainer, L2Distance>           PairDistances;
+typedef     PairDistances::DistanceType                             DistanceType;
+
+typedef     PairDistances::IndexType                                Vertex;
+typedef     std::vector<Vertex>                                     SubsampleVector;
+
+typedef     Simplex<Vertex>                                         Smplx;
+typedef     std::vector<Smplx>                                      SimplexVector;
+typedef     std::set<Smplx, Smplx::VertexDimensionComparison>       SimplexSet;
+
+typedef     std::vector<Vertex>                                     VertexVector;
+typedef     std::vector<DistanceType>                               EpsilonVector;
+typedef     std::vector<std::pair<Vertex, Vertex> >                 EdgeVector;
+
+typedef     Rips<PairDistances, Smplx>                              RipsGenerator;
+typedef     RipsGenerator::Evaluator                                SimplexEvaluator;
+
+struct      BirthInfo;
+typedef     ZigzagPersistence<BirthInfo>                            Zigzag;
+typedef     Zigzag::SimplexIndex                                    Index;
+typedef     Zigzag::Death                                           Death;
+typedef     std::map<Smplx, Index, 
+                            Smplx::VertexDimensionComparison>       Complex;
+typedef     Zigzag::ZColumn                                         Boundary;
+
+// Information we need to know when a class dies
+struct      BirthInfo
+{
+                    BirthInfo(Dimension d = 0, unsigned i = 0, bool u = false):
+                        dimension(d), index(i), un(u)               {}
+    
+    bool            operator<(const BirthInfo& other) const         { if (index == other.index) return (!un && other.un); else return index < other.index; }
+    bool            operator>(const BirthInfo& other) const         { return other.operator<(*this); }
+
+    Dimension       dimension;
+    unsigned        index;
+    bool            un;
+};
+std::ostream&       operator<<(std::ostream& out, const BirthInfo& bi)
+{ out << "(d=" << bi.dimension << ") " << bi.index; if (bi.un) out << " U " << (bi.index + 1); return out; }
+
+// Forward declarations of auxilliary functions
+void        report_death(std::ostream& out, const BirthInfo& birth, const BirthInfo& death);
+void        make_boundary(const Smplx& s, Complex& c, const Zigzag& zz, Boundary& b);
+void        process_command_line_options(int           argc,
+                                         char*         argv[],
+                                         unsigned&     ambient_dimension,
+                                         unsigned&     skeleton_dimension,
+                                         float&        epsilon,
+                                         std::string&  points_filename,
+                                         std::string&  subsample_filename,
+                                         std::string&  outfilename);
+void        add_simplex(const Smplx& s, const BirthInfo& birth, const BirthInfo& death, 
+                        Complex& complex, Zigzag& zz, std::ostream& out, Timer& add);
+void        remove_simplex(const Smplx& s, const BirthInfo& birth, const BirthInfo& death, 
+                           Complex& complex, Zigzag& zz, std::ostream& out, Timer& remove);
+
+int main(int argc, char* argv[])
+{
+#ifdef LOGGING
+    rlog::RLogInit(argc, argv);
+
+    stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+#endif
+
+    Timer total, remove, add, coface, generate;
+    total.start();
+ 
+#if 0
+    SetFrequency(cOperations, 25000);
+    SetTrigger(cOperations, cComplexSize);
+#endif
+
+    unsigned        ambient_dimension;
+    unsigned        skeleton_dimension;
+    float           epsilon;
+    std::string     points_filename, subsample_filename, outfilename;
+    process_command_line_options(argc, argv, ambient_dimension, skeleton_dimension, epsilon, 
+                                             points_filename, subsample_filename, outfilename);
+
+    // Read in points
+    PointContainer      points;
+    read_points(points_filename, points, ambient_dimension);
+    
+    // Read in subsamples
+    std::ifstream subsample_in(subsample_filename.c_str());
+    std::vector<SubsampleVector> subsamples;
+    subsamples.push_back(SubsampleVector());    //  Empty subsample to start from
+    while(subsample_in)
+    {
+        std::string line; std::getline(subsample_in, line);
+        std::istringstream line_in(line);
+        subsamples.push_back(SubsampleVector());
+        unsigned sample;
+        while (line_in >> sample)
+            subsamples.back().push_back(sample);
+    }
+    
+    std::cout << "Subsample size:" << std::endl;
+    for (unsigned i = 0; i < subsamples.size(); ++i)
+        std::cout << "  " << subsamples[i].size() << std::endl;
+
+    // Create output file
+    std::ofstream out(outfilename.c_str());
+
+    // Create pairwise distances
+    PairDistances distances(points);
+    
+
+    // Construct zigzag
+    Complex             complex;
+    Zigzag              zz;
+    RipsGenerator       rips(distances);
+    SimplexEvaluator    size(distances);
+    SimplexVector       subcomplex, across;
+    
+    rInfo("Commencing computation");
+    boost::progress_display show_progress(subsamples.size());
+    for (unsigned i = 0; i < subsamples.size() - 1; ++i)
+    {
+        // Take union of subsamples i and i+1
+        SubsampleVector forward_difference, backward_difference;
+        std::set_difference(subsamples[i+1].begin(), subsamples[i+1].end(),
+                            subsamples[i].begin(),   subsamples[i].end(),
+                            std::back_inserter(forward_difference));
+        std::set_difference(subsamples[i].begin(),   subsamples[i].end(),
+                            subsamples[i+1].begin(), subsamples[i+1].end(),
+                            std::back_inserter(backward_difference));
+        rInfo("Forward difference size:  %d", forward_difference.size());
+        rInfo("Backward difference size: %d", backward_difference.size());
+
+#if LOGGING
+        rDebug("Forward difference:");
+        for (SubsampleVector::const_iterator cur = forward_difference.begin(); cur != forward_difference.end(); ++cur)
+            rDebug("  %d", *cur);
+        
+        rDebug("Backward difference:");
+        for (SubsampleVector::const_iterator cur = backward_difference.begin(); cur != backward_difference.end(); ++cur)
+            rDebug("  %d", *cur);
+#endif
+
+        // Add simplices:       take star of the vertices in subsamples[i+1] - subsamples[i]
+        //                      by first computing the Rips complex on those vertices, and then 
+        //                      taking the star of each individual simplex (perhaps not the most 
+        //                      efficient procedure, but it will do for the first pass)
+        rInfo("Adding");
+        subcomplex.clear(); across.clear();
+        generate.start();
+        rips.generate(skeleton_dimension, epsilon, make_push_back_functor(subcomplex), forward_difference.begin(), forward_difference.end());
+        std::sort(subcomplex.begin(), subcomplex.end(), Smplx::VertexDimensionComparison());
+        generate.stop();
+        rInfo("  Subcomplex size: %d", subcomplex.size());
+        for (SimplexVector::const_iterator cur = subcomplex.begin(); cur != subcomplex.end(); ++cur)
+        {
+            add_simplex(*cur, BirthInfo(cur->dimension(), i, true), BirthInfo(cur->dimension() - 1, i), 
+                        complex, zz, out, add);
+
+            // Record cofaces with all other vertices outside of the subcomplex
+            coface.start();
+            rips.cofaces(*cur, skeleton_dimension, epsilon, make_push_back_functor(across), subsamples[i].begin(), subsamples[i].end());
+            coface.stop();
+        }
+        std::sort(across.begin(), across.end(), Smplx::VertexDimensionComparison());
+        rInfo("  Subcomplex simplices added");
+        rInfo("  Cross simplices size: %d", across.size());
+
+        // Add simplices that cut across VR(K_i) and VR(K_{i+1})
+        for (SimplexVector::const_iterator cur = across.begin(); cur != across.end(); ++cur)
+            add_simplex(*cur, BirthInfo(cur->dimension(), i, true), BirthInfo(cur->dimension() - 1, i), 
+                        complex, zz, out, add);
+        rInfo("  Cross simplices added");
+
+
+        // Remove simplices:    take star of the vertices in subsamples[i] - subsamples[i+1]
+        rInfo("Removing");
+        subcomplex.clear(); across.clear();
+        generate.start();
+        rips.generate(skeleton_dimension, epsilon, make_push_back_functor(subcomplex), backward_difference.begin(), backward_difference.end());
+        std::sort(subcomplex.begin(), subcomplex.end(), Smplx::VertexDimensionComparison());
+        generate.stop();
+        rInfo("  Subcomplex size: %d", subcomplex.size());
+
+        for (SimplexVector::const_iterator cur = subcomplex.begin(); cur != subcomplex.end(); ++cur)
+            rips.cofaces(*cur, skeleton_dimension, epsilon, make_push_back_functor(across), subsamples[i+1].begin(), subsamples[i+1].end());
+        std::sort(across.begin(), across.end(), Smplx::VertexDimensionComparison());
+        rInfo("  Cross simplices size: %d", across.size());
+
+        for (SimplexVector::const_reverse_iterator cur = across.rbegin(); cur != across.rend(); ++cur)
+            remove_simplex(*cur, BirthInfo(cur->dimension() - 1, i+1), BirthInfo(cur->dimension(), i, true),
+                           complex, zz, out, remove);
+        rInfo("  Cross simplices removed");
+
+        for (SimplexVector::const_reverse_iterator cur = subcomplex.rbegin(); cur != subcomplex.rend(); ++cur)
+            remove_simplex(*cur, BirthInfo(cur->dimension() - 1, i+1), BirthInfo(cur->dimension(), i, true),
+                           complex, zz, out, remove);
+        rInfo("  Subcomplex simplices removed");
+        
+        rInfo("Complex size: %d", complex.size());
+        
+        ++show_progress;
+    }
+
+    std::cout << std::endl;
+    total.stop();
+    remove.check("Remove timer          ");
+    add.check   ("Add timer             ");
+    coface.check("Coface timer          ");
+    total.check ("Total timer           ");
+}
+
+
+void        make_boundary(const Smplx& s, Complex& c, const Zigzag& zz, Boundary& b)
+{
+    // rDebug("  Boundary of <%s>", tostring(s).c_str());
+    for (Smplx::BoundaryIterator cur = s.boundary_begin(); cur != s.boundary_end(); ++cur)
+    {
+        b.append(c[*cur], zz.cmp);
+        // rDebug("   %d", c[*cur]->order);
+    }
+}
+
+void        report_death(std::ostream& out, const BirthInfo& birth, const BirthInfo& death)
+{
+    if (death > birth)
+        out << birth << " --- " << death << std::endl;
+}
+
+void        add_simplex(const Smplx& s, const BirthInfo& birth, const BirthInfo& death, 
+                        Complex& complex, Zigzag& zz, std::ostream& out, Timer& add)
+{
+    rDebug("Adding simplex: %s", tostring(s).c_str());
+
+    Index idx; Death d; Boundary b;
+    make_boundary(s, complex, zz, b);
+    add.start();
+    boost::tie(idx, d)  = zz.add(b, birth);
+    if (d) report_death(out, *d, death);
+    add.stop();
+    complex.insert(std::make_pair(s, idx));
+}
+
+void        remove_simplex(const Smplx& s, const BirthInfo& birth, const BirthInfo& death, 
+                           Complex& complex, Zigzag& zz, std::ostream& out, Timer& remove)
+{
+    rDebug("Removing simplex: %s", tostring(s).c_str());
+
+    Complex::iterator si = complex.find(s);
+    remove.start();
+    Death d = zz.remove(si->second, birth);
+    remove.stop();
+    if (d) report_death(out, *d, death);
+    complex.erase(si);
+}
+
+void        process_command_line_options(int           argc,
+                                         char*         argv[],
+                                         unsigned&     ambient_dimension,
+                                         unsigned&     skeleton_dimension,
+                                         float&        epsilon,
+                                         std::string&  points_filename,
+                                         std::string&  subsample_filename,
+                                         std::string&  outfilename)
+{
+    namespace po = boost::program_options;
+
+    po::options_description hidden("Hidden options");
+    hidden.add_options()
+        ("points-file",         po::value<std::string>(&points_filename),   "Point set whose Rips consistency zigzag we want to compute")
+        ("subsample-file",      po::value<std::string>(&subsample_filename),"Subsample list")
+        ("output-file",         po::value<std::string>(&outfilename),       "Location to save persistence pairs");
+    
+    po::options_description visible("Allowed options", 100);
+    visible.add_options()
+        ("help,h",                                                                              "produce help message")
+        ("ambient-dimsnion,a",  po::value<unsigned>(&ambient_dimension)->default_value(3),      "The ambient dimension of the point set")
+        ("skeleton-dimsnion,s", po::value<unsigned>(&skeleton_dimension)->default_value(2),     "Dimension of the Rips complex we want to compute")
+        ("epsilon,e",           po::value<float>(&epsilon)->default_value(0),                   "epsilon to use when computing the Rips complex");
+#if LOGGING
+    std::vector<std::string>    log_channels;
+    visible.add_options()
+        ("log,l",               po::value< std::vector<std::string> >(&log_channels),           "log channels to turn on (info, debug, etc)");
+#endif
+
+    po::positional_options_description pos;
+    pos.add("points-file", 1);
+    pos.add("subsample-file", 1);
+    pos.add("output-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);
+    po::notify(vm);
+
+#if LOGGING
+    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"
+     */
+#endif
+
+    if (vm.count("help") || !vm.count("points-file") || !vm.count("subsample-file") || !vm.count("output-file"))
+    { 
+        std::cout << "Usage: " << argv[0] << " [options] points-file subsample-file output-file" << std::endl;
+        std::cout << visible << std::endl; 
+        std::abort();
+    }
+}
--- a/examples/rips/l2distance.h	Mon Mar 23 11:40:29 2009 -0700
+++ b/examples/rips/l2distance.h	Tue Mar 24 15:17:39 2009 -0700
@@ -17,7 +17,7 @@
 {
     result_type     operator()(const Point& p1, const Point& p2) const
     {
-        AssertMsg(p1.size() == p2.size(), "Points must be in the same dimension (in L2Distance)");
+        AssertMsg(p1.size() == p2.size(), "Points must be in the same dimension (in L2Distance): dim1=%d, dim2=%d", p1.size(), p2.size());
         result_type sum = 0;
         for (size_t i = 0; i < p1.size(); ++i)
             sum += (p1[i] - p2[i])*(p1[i] - p2[i]);
@@ -36,6 +36,7 @@
         {
             double      x;
             in >> x;
+            if (!in) { points.pop_back(); break; }
             points.back().push_back(x);
         }
     }
--- a/examples/rips/rips-zigzag.cpp	Mon Mar 23 11:40:29 2009 -0700
+++ b/examples/rips/rips-zigzag.cpp	Tue Mar 24 15:17:39 2009 -0700
@@ -297,21 +297,6 @@
     }
 }
 
-bool        face_leaving_subcomplex(Complex::reverse_iterator si, const SimplexEvaluator& size, DistanceType after, DistanceType before)
-{
-    const Smplx& s = si->first;
-    for (Smplx::VertexContainer::const_iterator v1 = s.vertices().begin(); v1 != s.vertices().end(); ++v1)
-        for (Smplx::VertexContainer::const_iterator v2 = boost::next(v1); v2 != s.vertices().end(); ++v2)
-        {
-            Smplx e; e.add(*v1); e.add(*v2);
-            if (size(e) > after && size(e) <= before)
-                return true;
-        }
-
-    return false;
-}
-
-
 std::ostream&   operator<<(std::ostream& out, const BirthInfo& bi)
 { return (out << bi.distance); }
 
--- a/include/topology/rips.h	Mon Mar 23 11:40:29 2009 -0700
+++ b/include/topology/rips.h	Tue Mar 24 15:17:39 2009 -0700
@@ -53,8 +53,14 @@
         void                edge_cofaces(IndexType u, IndexType v, Dimension k, DistanceType max, const Functor& f, 
                                          Iterator candidates_begin, Iterator candidates_end) const;
         
+        // Calls functor f on all the simplices of the Rips complex that contain the given Simplex s
+        // (unlike the previous methods it does not call the functor on the Simplex s itself)
+        template<class Functor, class Iterator>
+        void                cofaces(const Simplex& s, Dimension k, DistanceType max, const Functor& f,
+                                    Iterator candidates_begin, Iterator candidates_end) const;
 
-        // No Iterator argument means IndexType and distances().begin() - distances().end()
+        
+        /* No Iterator argument means Iterator = IndexType and the range is [distances().begin(), distances().end()) */
         template<class Functor>
         void                generate(Dimension k, DistanceType max, const Functor& f) const
         { generate(k, max, f, boost::make_counting_iterator(distances().begin()), boost::make_counting_iterator(distances().end())); }
@@ -67,6 +73,10 @@
         void                edge_cofaces(IndexType u, IndexType v, Dimension k, DistanceType max, const Functor& f) const
         { edge_cofaces(u, v, k, max, f, boost::make_counting_iterator(distances().begin()), boost::make_counting_iterator(distances().end())); }
 
+        template<class Functor>
+        void                cofaces(const Simplex& s, Dimension k, DistanceType max, const Functor& f) const
+        { cofaces(s, k, max, f, boost::make_counting_iterator(distances().begin()), boost::make_counting_iterator(distances().end())); }
+
         
         const Distances&    distances() const                               { return distances_; }
         DistanceType        max_distance() const;
@@ -83,7 +93,8 @@
                                           typename VertexContainer::const_iterator  excluded,
                                           Dimension                                 max_dim,
                                           const NeighborTest&                       neighbor,
-                                          const Functor&                            functor) const;
+                                          const Functor&                            functor,
+                                          bool                                      check_initial = true) const;
         
     private:
         const Distances&    distances_;
--- a/include/topology/rips.hpp	Mon Mar 23 11:40:29 2009 -0700
+++ b/include/topology/rips.hpp	Tue Mar 24 15:17:39 2009 -0700
@@ -4,6 +4,7 @@
 #include <iostream>
 #include <utilities/log.h>
 #include <utilities/counter.h>
+#include <utilities/indirect.h>
 #include <boost/iterator/counting_iterator.hpp>
 #include <functional>
 
@@ -77,6 +78,42 @@
     bron_kerbosch(current, candidates, boost::prior(candidates.begin()), k, neighbor, f);
 }
 
+template<class D, class S>
+template<class Functor, class Iterator>
+void
+Rips<D,S>::
+cofaces(const Simplex& s, Dimension k, DistanceType max, const Functor& f, Iterator bg, Iterator end) const
+{
+    rLog(rlRipsDebug,   "In cofaces(%s)", tostring(s).c_str());
+
+    WithinDistance neighbor(distances(), max);
+
+    // current      = s.vertices()
+    VertexContainer current(s.vertices().begin(), s.vertices().end());
+    
+    // candidates   = everything - s.vertices()     that is a neighbor() of every vertex in the simplex
+    VertexContainer candidates;
+    typedef difference_iterator<Iterator, 
+                                typename VertexContainer::const_iterator, 
+                                std::less<Vertex> >                     DifferenceIterator;
+    for (DifferenceIterator cur =  DifferenceIterator(bg, end, s.vertices().begin(), s.vertices().end()); 
+                            cur != DifferenceIterator(end, end, s.vertices().end(), s.vertices().end()); 
+                            ++cur)
+    {
+        bool nghbr = true;
+        for (typename VertexContainer::const_iterator v = s.vertices().begin(); v != s.vertices().end(); ++v)
+            if (!neighbor(*v, *cur))        { nghbr = false; break; }
+
+        if (nghbr)  
+        {
+            candidates.push_back(*cur);
+            rLog(rlRipsDebug,   "  added candidate: %d", *cur);
+        }
+    }
+
+    bron_kerbosch(current, candidates, boost::prior(candidates.begin()), k, neighbor, f, false);
+}
+
 
 template<class D, class S>
 template<class Functor, class NeighborTest>
@@ -87,11 +124,12 @@
               typename VertexContainer::const_iterator  excluded,
               Dimension                                 max_dim,    
               const NeighborTest&                       neighbor,       
-              const Functor&                            functor) const
+              const Functor&                            functor,
+              bool                                      check_initial) const
 {
     rLog(rlRipsDebug,       "Entered bron_kerbosch");
     
-    if (!current.empty())
+    if (check_initial && !current.empty())
     {
         Simplex s(current);
         rLog(rlRipsDebug,   "Reporting simplex: %s", tostring(s).c_str());