Merged in Rips complexes dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Mon, 15 Sep 2008 16:41:18 -0700
branchdev
changeset 93 824cf6d63af0
parent 90 dea0e9726c62 (current diff)
parent 92 75b66d8dfa7c (diff)
child 94 900f42b643de
Merged in Rips complexes
--- a/examples/CMakeLists.txt	Mon Sep 15 16:33:34 2008 -0700
+++ b/examples/CMakeLists.txt	Mon Sep 15 16:41:18 2008 -0700
@@ -8,3 +8,4 @@
 add_subdirectory			(grid)
 add_subdirectory			(triangle)
 add_subdirectory			(poincare)
+add_subdirectory			(rips)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/rips/CMakeLists.txt	Mon Sep 15 16:41:18 2008 -0700
@@ -0,0 +1,7 @@
+set							(targets						
+							 rips)
+							 
+foreach 					(t ${targets})
+	add_executable			(${t} ${t}.cpp)
+	target_link_libraries	(${t} ${libraries})
+endforeach 					(t ${targets})
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/rips/rips.cpp	Mon Sep 15 16:41:18 2008 -0700
@@ -0,0 +1,32 @@
+#include <topology/rips.h>
+
+struct Distances
+{
+    typedef         int             IndexType;
+    typedef         double          DistanceType;
+
+    DistanceType    operator()(IndexType a, IndexType b) const      { return std::abs(a - b); }
+
+    size_t          size() const                                    { return 2000; }
+    IndexType       begin() const                                   { return 0; }
+    IndexType       end() const                                     { return size(); }
+};
+
+int main(int argc, char* argv[])
+{
+#ifdef LOGGING
+	rlog::RLogInit(argc, argv);
+
+	stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
+	stdoutLog.subscribeTo( RLOG_CHANNEL("rips/info") );
+#endif
+
+    Distances distances;
+    Rips<Distances> rips(distances);
+
+    //rips.generate(3, distances.size());
+    rips.generate(3, 10);
+    //rips.print();
+    
+    std::cout << "Size: " << rips.size() << std::endl;
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/rips.h	Mon Sep 15 16:41:18 2008 -0700
@@ -0,0 +1,48 @@
+#ifndef __RIPS_H__
+#define __RIPS_H__
+
+#include <vector>
+#include <string>
+#include <topology/simplex.h>
+
+/**
+ * Rips complex class
+ *
+ * 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_ = SimplexWithVertices<typename Distances_::IndexType> >
+class Rips
+{
+    public:
+        typedef             Distances_                                      Distances; 
+        typedef             typename Distances::IndexType                   IndexType;
+        typedef             typename Distances::DistanceType                DistanceType;
+
+        typedef             Simplex_                                        Simplex;
+        typedef             std::vector<Simplex>                            SimplexVector;
+
+    public:
+                            Rips(const Distances& distances): 
+                                distances_(distances)                       {}
+
+        void                generate(Dimension k, DistanceType max);        /// generate k-skeleton of the Rips complex
+
+        void                print() const;
+        size_t              size() const                                    { return simplices.size(); }
+
+    private:
+        struct              ComparePair;
+
+    private:
+        const Distances&    distances_;
+
+        SimplexVector       simplices;
+
+};
+
+#include "rips.hpp"
+
+#endif // __RIPS_H__
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/rips.hpp	Mon Sep 15 16:41:18 2008 -0700
@@ -0,0 +1,101 @@
+#include <algorithm>
+#include <utility>
+#include <boost/utility.hpp>
+#include <iostream>
+#include <utilities/log.h>
+
+#ifdef LOGGING
+static rlog::RLogChannel* rlRips =                  DEF_CHANNEL("rips/info", rlog::Log_Debug);
+#endif // LOGGING
+
+#ifdef COUNTERS
+static Counter*  cClique =                          GetCounter("rips/clique");
+#endif // COUNTERS
+
+template<class Distances_, class Simplex_>
+struct Rips<Distances_, Simplex_>::ComparePair
+{
+                            ComparePair(const Distances& distances): 
+                                distances_(distances)                       {}
+
+        bool                operator()(const std::pair<IndexType, IndexType>& a,
+                                       const std::pair<IndexType, IndexType>& b)            {  return   distances_(a.first, a.second) <
+                                                                                                        distances_(b.first, b.second);  }
+
+        const Distances&    distances_;
+};
+
+template<class DistanceType_, class Simplex_>
+void
+Rips<DistanceType_, Simplex_>::
+generate(Dimension k, DistanceType max)
+{
+    // Order all the edges
+    typedef std::vector< std::pair<IndexType, IndexType> >      EdgeVector;
+    EdgeVector      edges;
+    for (IndexType a = distances_.begin(); a != distances_.end(); ++a)
+    {
+        Simplex ssx; ssx.add(a);
+        simplices.push_back(ssx);
+        for (IndexType b = boost::next(a); b != distances_.end(); ++b)
+        {
+            if (distances_(a,b) <= max)
+                edges.push_back(std::make_pair(a,b));
+        }
+    }
+    std::sort(edges.begin(), edges.end(), ComparePair(distances_));
+
+    // Generate simplices
+    std::vector<std::vector<size_t> >       vertex_star(distances_.size());
+    for(typename EdgeVector::const_iterator cur = edges.begin(); cur != edges.end(); ++cur)
+    {
+        rLog(rlRips, "Current edge: %d %d", cur->first, cur->second);
+
+        // Create the edge
+        Simplex edge; edge.add(cur->first); edge.add(cur->second);
+        simplices.push_back(edge);
+        vertex_star[cur->first].push_back(simplices.size() - 1); 
+        vertex_star[cur->second].push_back(simplices.size() - 1);
+
+        // Go through a star
+        size_t sz = vertex_star[cur->first].size() - 1;
+        for (size_t i = 0; i < sz; ++i)
+        {
+            const Simplex& ssx = simplices[vertex_star[cur->first][i]];
+            rLog(rlRips, "  %s", tostring(ssx).c_str());
+            bool accept = true;
+            for (typename Simplex::VertexContainer::const_iterator v = ssx.vertices().begin(); v != ssx.vertices().end(); ++v)
+            {
+                if (*v == cur->first) continue;
+                
+                if (  distances_(*v, cur->second) >  distances_(cur->first, cur->second) ||
+                    ((distances_(*v, cur->second) == distances_(cur->first, cur->second)) && 
+                     (*v > cur->first)))
+                {
+                    accept = false;
+                    break;
+                }
+            }
+            if (accept)
+            {
+                Simplex tsx(ssx); tsx.add(cur->second);
+                simplices.push_back(tsx);
+                rLog(rlRips, "  Accepting: %s", tostring(tsx).c_str());
+         
+                // Update stars
+                if (tsx.dimension() < k)
+                    for (typename Simplex::VertexContainer::const_iterator v = tsx.vertices().begin(); v != tsx.vertices().end(); ++v)
+                        vertex_star[*v].push_back(simplices.size() - 1);
+            }
+        }
+    }
+}
+
+template<class Distances_, class Simplex_>
+void
+Rips<Distances_, Simplex_>::
+print() const
+{
+    for (typename SimplexVector::const_iterator cur = simplices.begin(); cur != simplices.end(); ++cur)
+        std::cout << *cur << std::endl;
+}
--- a/include/topology/simplex.h	Mon Sep 15 16:33:34 2008 -0700
+++ b/include/topology/simplex.h	Mon Sep 15 16:41:18 2008 -0700
@@ -60,6 +60,7 @@
 		bool					contains(const Vertex& v) const;
 		const VertexContainer&	vertices() const									{ return vertices_; }
 		void					add(const Vertex& v);
+        void                    join(const Self& other);
 		/// @}
 
 		/// \name Assignment and comparison
--- a/include/topology/simplex.hpp	Mon Sep 15 16:33:34 2008 -0700
+++ b/include/topology/simplex.hpp	Mon Sep 15 16:41:18 2008 -0700
@@ -38,6 +38,15 @@
 SimplexWithVertices<V>::
 add(const Vertex& v)
 { vertices_.push_back(v); std::sort(vertices_.begin(), vertices_.end()); }
+	
+template<class V>
+void
+SimplexWithVertices<V>::
+join(const Self& other)
+{ 
+    vertices_.insert(vertices_.end(), other.vertices_.begin(), other.vertices_.end());
+    std::sort(vertices_.begin(), vertices_.end()); 
+}
 
 template<class V>
 std::ostream&