Merged upstream dev
authorChristos Mantoulidis <cmad@stanford.edu>
Tue, 28 Jul 2009 11:42:03 -0700
branchdev
changeset 154 8af75817ecbc
parent 153 7731c42892de (diff)
parent 142 ae2b1702c936 (current diff)
child 155 1dde8a05fcd7
Merged upstream
--- a/examples/rips/CMakeLists.txt	Fri Jul 24 11:08:48 2009 -0700
+++ b/examples/rips/CMakeLists.txt	Tue Jul 28 11:42:03 2009 -0700
@@ -1,6 +1,7 @@
 set                         (targets                        
                              rips
                              rips-pairwise
+                             rips-weighted
                              rips-image-zigzag
                              rips-zigzag)
                              
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/rips/rips-weighted.cpp	Tue Jul 28 11:42:03 2009 -0700
@@ -0,0 +1,145 @@
+#include <topology/weighted-rips.h>
+#include <topology/filtration.h>
+#include <topology/static-persistence.h>
+#include <topology/dynamic-persistence.h>
+#include <topology/persistence-diagram.h>
+
+//#define RIPS_CLOSURE_CECH_SKELETON
+#ifndef RIPS_CLOSURE_CECH_SKELETON
+#include <geometry/weighted-l2distance.h>
+#else
+#include <geometry/weighted-cechdistance.h>
+#endif
+
+#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") );
+#endif
+
+    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);
+    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(), DataDimensionComparison<Smplx>());
+
+    Timer persistence_timer; persistence_timer.start();
+    Persistence p(f);
+    p.pair_simplices();
+    persistence_timer.stop();
+
+    // 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");
+    hidden.add_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);
+    visible.add_options()
+        ("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");
+#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("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);
+    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("input-file"))
+    { 
+        std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl;
+        std::cout << visible << std::endl; 
+        std::abort();
+    }
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/weighted-cechdistance.h	Tue Jul 28 11:42:03 2009 -0700
@@ -0,0 +1,55 @@
+#ifndef __L2_DISTANCE_H__
+#define __L2_DISTANCE_H__
+
+#include <utilities/types.h>
+
+#include <vector>
+#include <fstream>
+#include <functional>
+#include <cmath>
+#include <string>
+#include <sstream>
+
+
+typedef     std::vector<double>                                     Point;
+typedef     std::vector<Point>                                      PointContainer;
+
+struct WeightedL2Distance:
+    public std::binary_function<const Point&, const Point&, double>
+{
+    result_type     operator()(const Point& p1, const Point& p2) const
+    {
+        AssertMsg(p1.size() == p2.size(), "Points must be in the same dimension (in L2Distance): dim1=%d, dim2=%d", p1.size()-1, p2.size()-1);
+
+        /* the distance of a point to itself is the radius at which the "power distance" ball around it appears:
+           d(p,p) := sqrt(-w_p) */
+        if (p1 == p2)
+            return sqrt(-p1[p1.size()-1]); 
+
+        /* otherwise, the distance is the radius at which the power balls around p, q intersect:
+           d(p,q) := sqrt( ( (w_p-w_q)^2 / 4||p-q||^2 ) - (w_p+w_q) / 2 + ||p-q||^2 / 4 ) */
+        result_type sq_l2dist = 0;
+        result_type wp = p1[p1.size()-1];
+        result_type wq = p2[p2.size()-1];
+        for (size_t i = 0; i < p1.size()-1; ++i)
+            sq_l2dist += (p1[i] - p2[i])*(p1[i] - p2[i]);
+        return sqrt( (wp-wq)*(wp-wq) / (4 * sq_l2dist) - (wp+wq) / 2 + sq_l2dist / 4 );
+    }
+};
+
+void    read_weighted_points(const std::string& infilename, PointContainer& points)
+{
+    std::ifstream in(infilename.c_str());
+    std::string   line;
+    while(std::getline(in, line))
+    {
+        if (line[0] == '#') continue;               // comment line in the file
+        std::stringstream linestream(line);
+        double x;
+        points.push_back(Point());
+        while (linestream >> x)
+            points.back().push_back(x);
+    }
+}
+
+#endif // __L2_DISTANCE_H__
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/weighted-l2distance.h	Tue Jul 28 11:42:03 2009 -0700
@@ -0,0 +1,53 @@
+#ifndef __L2_DISTANCE_H__
+#define __L2_DISTANCE_H__
+
+#include <utilities/types.h>
+
+#include <vector>
+#include <fstream>
+#include <functional>
+#include <cmath>
+#include <string>
+#include <sstream>
+
+
+typedef     std::vector<double>                                     Point;
+typedef     std::vector<Point>                                      PointContainer;
+
+struct WeightedL2Distance:
+    public std::binary_function<const Point&, const Point&, double>
+{
+    result_type     operator()(const Point& p1, const Point& p2) const
+    {
+        AssertMsg(p1.size() == p2.size(), "Points must be in the same dimension (in L2Distance): dim1=%d, dim2=%d", p1.size()-1, p2.size()-1);
+
+        /* the distance of a point to itself is the radius at which the "power distance" ball around it appears:
+           d(p,p) := sqrt(-w_p) */
+        if (p1 == p2)
+            return sqrt(-p1[p1.size()-1]);
+
+        /* otherwise, the distance is the radius at which the power balls around p, q intersect:
+           d(p,q) := sqrt( ||p-q||^2 - w_p - w_q ) */
+        result_type sq_l2dist = 0;
+        for (size_t i = 0; i < p1.size()-1; ++i)
+            sq_l2dist += (p1[i] - p2[i])*(p1[i] - p2[i]);
+        return sqrt(sq_l2dist - p1[p1.size()-1] - p2[p2.size()-1]);
+    }
+};
+
+void    read_weighted_points(const std::string& infilename, PointContainer& points)
+{
+    std::ifstream in(infilename.c_str());
+    std::string   line;
+    while(std::getline(in, line))
+    {
+        if (line[0] == '#') continue;               // comment line in the file
+        std::stringstream linestream(line);
+        double x;
+        points.push_back(Point());
+        while (linestream >> x)
+            points.back().push_back(x);
+    }
+}
+
+#endif // __L2_DISTANCE_H__
--- a/include/topology/rips.h	Fri Jul 24 11:08:48 2009 -0700
+++ b/include/topology/rips.h	Tue Jul 28 11:42:03 2009 -0700
@@ -84,7 +84,7 @@
         DistanceType        distance(const Simplex& s1, const Simplex& s2) const;
 
 
-    private:
+    protected:
         class               WithinDistance;
 
         template<class Functor, class NeighborTest>
@@ -96,7 +96,7 @@
                                           const Functor&                            functor,
                                           bool                                      check_initial = true) const;
         
-    private:
+    protected:
         const Distances&    distances_;
 };
         
@@ -111,7 +111,7 @@
 
         bool                operator()(Vertex u, Vertex v) const                        { return distances_(u, v) <= max_; }
 
-    private:
+    protected:
         const Distances&    distances_;  
         DistanceType        max_;
 };
@@ -127,7 +127,7 @@
 
         DistanceType        operator()(const Simplex& s) const;
 
-    private:
+    protected:
         const Distances&    distances_;
 };
 
@@ -150,7 +150,7 @@
             return e1 < e2;
         }
 
-    private:
+    protected:
         Evaluator           eval_;
 };
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/weighted-rips.h	Tue Jul 28 11:42:03 2009 -0700
@@ -0,0 +1,114 @@
+#ifndef __WEIGHTED_RIPS_H__
+#define __WEIGHTED_RIPS_H__
+
+#include <vector>
+#include <string>
+#include "simplex.h"
+#include "rips.h"
+#include <boost/iterator/counting_iterator.hpp>
+
+/**
+ * WeightedRips class
+ *
+ * Class providing basic operations to work with Rips complexes. It implements Bron-Kerbosch algorithm, 
+ * and provides simple wrappers for various functions.
+ *
+ * Distances_ is expected to define types IndexType and DistanceType as well as 
+ *               provide operator()(...) which given two IndexTypes should return 
+ *               the distance between them. There should be methods begin() and end() 
+ *               for iterating over IndexTypes as well as a method size().
+ */
+template<class Distances_, class Simplex_ = Simplex<typename Distances_::IndexType, typename Distances_::DistanceType> >
+class WeightedRips : public Rips<Distances_, Simplex_>
+{
+    public:
+
+        /* redeclaring the typedefs because they cannot be inherited at compile-time */
+        typedef             Distances_                                      Distances; 
+        typedef             typename Distances::IndexType                   IndexType;
+        typedef             typename Distances::DistanceType                DistanceType;
+
+        typedef             Simplex_                                        Simplex;
+        typedef             typename Simplex::Vertex                        Vertex;             // should be the same as IndexType
+        typedef             typename Simplex::VertexContainer               VertexContainer;
+
+        class               Evaluator;
+        class               Comparison;
+
+    public:
+                            WeightedRips(const Distances& distances):
+                                Rips<Distances_, Simplex_>(distances)                             {}
+
+                            template<class Functor>
+                            void generate(Dimension k, DistanceType max, const Functor& f) const;
+
+};
+
+/**
+ * DistanceDataStackingFunctor class
+ * 
+ * Class providing a functor that is to be called by WeightedRips::generate(). This functor
+ * takes as an argument (to its constructor) the original functor passed by the user to
+ * generate(), and a new ``double'' functor is created. Assuming that the functor acts on
+ * simplices, first the value of the simplex is computed (the radius at which the simplex
+ * appears in the weighted Rips complex), the data field of the simplex is populated with
+ * this value, and then the original functor is called (it has no idea that it was
+ * intercepted).
+ */
+
+template<class Rips_, class Functor_>
+class DistanceDataStackingFunctor
+{
+	public:
+		typedef     typename Rips_::Simplex         Simplex_;
+
+		DistanceDataStackingFunctor(const Rips_ &r, const Functor_ &f):
+			rips(r), original_functor(f) { }
+
+		void operator()(const Simplex_ &s) const
+		{
+			Simplex_         s_new(s);
+			s_new.data()   = rips.distance(s_new, s_new);
+			original_functor(s_new);
+		}
+
+	private:
+		const Rips_       &rips;
+		const Functor_    &original_functor;
+};
+
+template<class Distances_, class Simplex_>
+template<class Functor>
+void WeightedRips<Distances_,Simplex_>::generate(Dimension k, DistanceType max, const Functor &f) const
+{
+	Rips<Distances_,Simplex_>::generate(k, max, DistanceDataStackingFunctor<WeightedRips<Distances_,Simplex_>,Functor>(*this, f));
+}
+
+template<class Distances_, class Simplex_>
+class WeightedRips<Distances_, Simplex_>::Evaluator: public Rips<Distances_, Simplex_>::Evaluator
+{
+    public:
+                            Evaluator(const Distances& distances): 
+                                Rips<Distances_, Simplex_>::Evaluator(distances)                       {}
+
+        DistanceType       operator()(const Simplex& s) const { return s.data(); }
+};
+
+template<class Distances_, class Simplex_>
+class WeightedRips<Distances_, Simplex_>::Comparison: public Rips<Distances_, Simplex_>::Comparison
+{
+    public:
+                            Comparison(const Distances& distances):
+                                Rips<Distances_, Simplex_>::Comparison(distances)                            {}
+
+        bool                operator()(const Simplex& s1, const Simplex& s2) const    
+        { 
+                            if (s1.dimension() != s2.dimension())
+                                return s1.dimension() < s2.dimension();
+                            return s1.data() < s2.data();
+        }
+};
+
+#include "weighted-rips.hpp"
+
+#endif // __WEIGHTED_RIPS_H__
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/weighted-rips.hpp	Tue Jul 28 11:42:03 2009 -0700
@@ -0,0 +1,18 @@
+#include <algorithm>
+#include <utility>
+#include <boost/utility.hpp>
+#include <iostream>
+#include <utilities/log.h>
+#include <utilities/counter.h>
+#include <utilities/indirect.h>
+#include <boost/iterator/counting_iterator.hpp>
+#include <functional>
+
+#ifdef LOGGING
+static rlog::RLogChannel* rlWeightedRips =                  DEF_CHANNEL("weightedrips/info", rlog::Log_Debug);
+static rlog::RLogChannel* rlWeightedRipsDebug =             DEF_CHANNEL("weightedrips/debug", rlog::Log_Debug);
+#endif // LOGGING
+
+#ifdef COUNTERS
+static Counter*  cClique =                          GetCounter("weightedrips/clique");
+#endif // COUNTERS