Added Python bindings for the Rips complex + an example dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Sat, 11 Apr 2009 10:29:53 -0700
branchdev
changeset 126 3c3e77ac43d2
parent 125 0a2c2283e4a8
child 127 406c6cc00b9c
Added Python bindings for the Rips complex + an example
bindings/python/CMakeLists.txt
bindings/python/chain.cpp
bindings/python/dionysus.cpp
bindings/python/dionysus/__init__.py
bindings/python/dionysus/distances.py
bindings/python/python-rips.h
bindings/python/rips.cpp
examples/rips/CMakeLists.txt
examples/rips/rips-pairwise.py
examples/rips/rips.py
--- a/bindings/python/CMakeLists.txt	Fri May 01 15:09:38 2009 -0700
+++ b/bindings/python/CMakeLists.txt	Sat Apr 11 10:29:53 2009 -0700
@@ -12,6 +12,16 @@
                                                 static-persistence.cpp
                                                 simplex.cpp
                                                 zigzag-persistence.cpp
+                                                rips.cpp
                             )
+target_link_libraries       (_dionysus ${libraries})
+
 
-target_link_libraries       (_dionysus ${libraries})
+add_custom_target           (dionysus ALL
+                             ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR}/dionysus ${CMAKE_CURRENT_BINARY_DIR}/dionysus
+                             DEPENDS dionysus/__init__.py)
+
+get_target_property         (_dionysus_location _dionysus LOCATION)
+add_custom_target           (dionysus-link ALL 
+                             ${CMAKE_COMMAND} -E create_symlink ${_dionysus_location} ${CMAKE_CURRENT_BINARY_DIR}/dionysus/_dionysus.so
+                             DEPENDS _dionysus)
--- a/bindings/python/chain.cpp	Fri May 01 15:09:38 2009 -0700
+++ b/bindings/python/chain.cpp	Sat Apr 11 10:29:53 2009 -0700
@@ -4,7 +4,7 @@
 
 using namespace boost::python;
 
-typedef     SPersistence::OrderDescriptor::Chains::Chain        VChain;
+typedef     SPersistence::Chain                                 VChain;
 
 
 VChain::const_iterator      begin(const VChain& c)              { return c.begin(); }
--- a/bindings/python/dionysus.cpp	Fri May 01 15:09:38 2009 -0700
+++ b/bindings/python/dionysus.cpp	Sat Apr 11 10:29:53 2009 -0700
@@ -9,6 +9,8 @@
 void export_chain();
 void export_zigzag_persistence();
 
+void export_rips();
+
 #ifdef LOGGING
 void            enable_log(std::string s)
 {
@@ -25,6 +27,8 @@
 
     export_zigzag_persistence();
 
+    export_rips();
+
 #ifdef LOGGING
     bp::def("enable_log",           &enable_log);
 #endif
--- a/bindings/python/dionysus/__init__.py	Fri May 01 15:09:38 2009 -0700
+++ b/bindings/python/dionysus/__init__.py	Sat Apr 11 10:29:53 2009 -0700
@@ -1,4 +1,5 @@
-from _dionysus import *
+from    _dionysus   import *
+from    distances   import *
 
 
 #def init_with_data(s,v, d = None):
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/bindings/python/dionysus/distances.py	Sat Apr 11 10:29:53 2009 -0700
@@ -0,0 +1,41 @@
+from    math        import sqrt
+
+def l2(p):
+    return sqrt(                                        # take the square root
+            reduce(lambda x, y: x + y,                  # add them all up
+                   map(lambda x: x**2, p)))             # square each coordinate
+
+# Pairwise distances between the elements of `points` with respect to some `norm`
+class PairwiseDistances:
+    def __init__(self, points, norm = l2):
+        self.points = points
+        self.norm = norm
+
+    def __len__(self):
+        return len(self.points)
+
+    def __call__(self, p1, p2):
+        return self.norm([x - y for (x,y) in zip(self.points[p1], self.points[p2])])
+
+# Caches all distances specified by `distances`
+class ExplicitDistances:
+    def __init__(self, distances):
+        self.len = len(distances)
+        self.distances = []
+        for i in xrange(self.len): 
+            self.distances.append([])
+            for j in xrange(self.len):
+                self.distances[-1].append(distances(i,j))
+
+    def __len__(self):
+        return self.len
+
+    def __call__(self, p1, p2):
+        return self.distances[p1][p2]
+
+# Generator of all points in a file `filename` with one point per line
+def points_file(filename):
+    fd = open(filename)
+    for line in fd.xreadlines():
+        yield map(float, line.strip().split())
+    fd.close()
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/bindings/python/python-rips.h	Sat Apr 11 10:29:53 2009 -0700
@@ -0,0 +1,108 @@
+#ifndef __PYTHON_RIPS_H__
+#define __PYTHON_RIPS_H__
+
+#include <topology/rips.h>
+#include <utilities/indirect.h>
+
+#include "python-simplex.h"
+
+#include <boost/python.hpp>
+#include <boost/python/stl_iterator.hpp>
+
+
+namespace bp = boost::python;
+
+
+
+// This strange wrapper is necessary because Rips<...> stores only a const reference to the distances. 
+// Something on the C++ side of things must store the actual DistancesWrapper object, and that's the 
+// purpose of this class.
+class RipsWithDistances
+{
+    public:
+        class DistancesWrapper
+        {
+            public:
+                typedef             unsigned                                        IndexType;
+                typedef             double                                          DistanceType;
+        
+                                    DistancesWrapper(bp::object distances):
+                                        distances_(distances)                       {}
+        
+                DistanceType        operator()(IndexType a, IndexType b) const      { return bp::extract<DistanceType>(distances_(a, b)); }
+        
+                IndexType           size() const                                    { return bp::len(distances_); }
+                IndexType           begin() const                                   { return 0; }
+                IndexType           end() const                                     { return size(); }
+        
+            private:
+                bp::object          distances_;
+        };
+
+        typedef             DistancesWrapper::IndexType                             IndexType;
+        typedef             DistancesWrapper::DistanceType                          DistanceType;
+
+        typedef             Rips<DistancesWrapper, SimplexVD>                       RipsDS;
+        typedef             RipsDS::Comparison                                      Comparison;
+        typedef             RipsDS::Evaluator                                       Evaluator;
+
+        class FunctorWrapper
+        {
+            public:
+                                    FunctorWrapper(bp::object functor):
+                                        functor_(functor)                           {}
+        
+                void                operator()(const RipsDS::Simplex& s) const      { functor_(s); }
+        
+            private:
+                bp::object          functor_;
+        };
+
+
+                            RipsWithDistances(bp::object distances):
+                                distances_(distances), rips_(distances_),
+                                cmp_(Comparison(distances_)), eval_(distances_)     {}
+
+        void                generate(Dimension k, DistanceType max, bp::object functor) const
+        { rips_.generate(k, max, FunctorWrapper(functor)); }
+
+        void                vertex_cofaces(IndexType v, Dimension k, DistanceType max, bp::object functor) const
+        { rips_.vertex_cofaces(v, k, max, FunctorWrapper(functor)); }
+        
+        void                edge_cofaces(IndexType u, IndexType v, Dimension k, DistanceType max, bp::object functor) const
+        { rips_.edge_cofaces(u, v, k, max, FunctorWrapper(functor)); }
+
+        void                generate_candidates(Dimension k, DistanceType max, bp::object functor, bp::object seq) const
+        { 
+            rips_.generate(k, max, FunctorWrapper(functor), 
+                           bp::stl_input_iterator<IndexType>(seq), bp::stl_input_iterator<IndexType>());
+        }
+        
+        void                vertex_cofaces_candidate(IndexType v, Dimension k, DistanceType max, 
+                                                     bp::object functor, bp::object seq) const
+        { 
+            rips_.vertex_cofaces(v, k, max, FunctorWrapper(functor), 
+                                 bp::stl_input_iterator<IndexType>(seq), bp::stl_input_iterator<IndexType>());
+        }
+        
+        void                edge_cofaces_candidates(IndexType u, IndexType v, Dimension k, DistanceType max, 
+                                                    bp::object functor, bp::object seq) const
+        { 
+            rips_.edge_cofaces(u, v, k, max, FunctorWrapper(functor), 
+                               bp::stl_input_iterator<IndexType>(seq), bp::stl_input_iterator<IndexType>());
+        }
+
+        int                 cmp(const SimplexVD& s1, const SimplexVD& s2) const     
+        { return cmp_.compare(s1, s2); }
+        
+        DistanceType        eval(const SimplexVD& s) const     
+        { return eval_(s); }
+        
+    private:
+        DistancesWrapper                            distances_;
+        RipsDS                                      rips_;
+        ThreeOutcomeCompare<Comparison>             cmp_;           // in Python, cmp is a three outcome comparison
+        Evaluator                                   eval_;
+};
+
+#endif
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/bindings/python/rips.cpp	Sat Apr 11 10:29:53 2009 -0700
@@ -0,0 +1,31 @@
+#include <topology/rips.h>
+#include <boost/python.hpp>
+#include "python-rips.h"            // defines RipsWithDistances
+
+#include <iostream>
+
+namespace bp = boost::python;
+
+/* Various wrappers for exposing Rips to Python */
+// Constructor from distances
+boost::shared_ptr<RipsWithDistances>        init_from_distances(bp::object distances)
+{ 
+    boost::shared_ptr<RipsWithDistances>    p(new RipsWithDistances(distances));
+    return p;
+}
+
+void export_rips()
+{
+    bp::class_<RipsWithDistances>("Rips", bp::no_init)
+        .def("__init__",            make_constructor(&init_from_distances))
+        .def("generate",            &RipsWithDistances::generate)
+        .def("generate",            &RipsWithDistances::generate_candidates)
+        .def("vertex_cofaces",      &RipsWithDistances::vertex_cofaces)
+        .def("vertex_cofaces",      &RipsWithDistances::vertex_cofaces_candidate)
+        .def("edge_cofaces",        &RipsWithDistances::edge_cofaces)
+        .def("edge_cofaces",        &RipsWithDistances::edge_cofaces_candidates)
+
+        .def("cmp",                 &RipsWithDistances::cmp)
+        .def("eval",                &RipsWithDistances::eval)
+    ;
+}
--- a/examples/rips/CMakeLists.txt	Fri May 01 15:09:38 2009 -0700
+++ b/examples/rips/CMakeLists.txt	Sat Apr 11 10:29:53 2009 -0700
@@ -8,3 +8,8 @@
     add_executable          (${t} ${t}.cpp)
     target_link_libraries   (${t} ${libraries} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${Boost_SERIALIZATION_LIBRARY})
 endforeach                  (t ${targets})
+
+add_custom_target           (rips.py ALL 
+                             ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/rips.py ${CMAKE_CURRENT_BINARY_DIR}/rips.py)
+add_custom_target           (rips-pairwise.py ALL 
+                             ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/rips-pairwise.py ${CMAKE_CURRENT_BINARY_DIR}/rips-pairwise.py)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/rips/rips-pairwise.py	Sat Apr 11 10:29:53 2009 -0700
@@ -0,0 +1,40 @@
+#/usr/bin/env python
+
+
+from    dionysus    import Rips, PairwiseDistances, StaticPersistence, Filtration, points_file, \
+                           ExplicitDistances, data_dim_cmp
+from    sys         import argv, exit
+import  time
+
+
+if len(argv) < 4:
+    print "Usage: %s POINTS SKELETON MAX" % argv[0]
+    exit()
+
+filename = argv[1]
+skeleton = int(argv[2])
+max = float(argv[3])
+
+points = [p for p in points_file(filename)]
+distances = PairwiseDistances(points)
+#distances = ExplicitDistances(distances)           # speeds up generation of the Rips complex at the expense of memory usage
+rips = Rips(distances)
+print time.asctime(), "Rips initialized"
+
+simplices = []
+rips.generate(skeleton, max, simplices.append)
+print time.asctime(), "Generated complex: %d simplices" % len(simplices)
+
+# While this step is unnecessary (Filtration below can be passed rips.cmp), 
+# it greatly speeds up the running times
+for s in simplices: s.data = rips.eval(s)
+print time.asctime(), simplices[0], '...', simplices[-1]
+
+f = Filtration(simplices, data_dim_cmp)             # could be rips.cmp if s.data for s in simplices is not set
+print time.asctime(), "Set up filtration"
+
+p = StaticPersistence(f)
+print time.asctime(), "Initialized StaticPersistence"
+
+p.pair_simplices()
+print time.asctime(), "Simplices paired"
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/rips/rips.py	Sat Apr 11 10:29:53 2009 -0700
@@ -0,0 +1,40 @@
+from    math        import fabs
+from    dionysus    import Rips, Filtration, StaticPersistence #, enable_log
+
+# Simple minded pairwise distance functor distance
+class Distances:
+    def __len__(self):
+        return 5
+
+    def __call__(self, x, y):
+        return fabs(y-x)
+
+dist = Distances()
+r = Rips(dist)
+lst = []
+lst2 = []
+
+#enable_log('rips')
+
+r.generate(1, 3, lst.append)
+r.generate(1, 3, lst2.append, [0,2,4])
+
+print "Rips complex on all vertices:", lst
+print "Rips complex on vertices [0,2,4]):", lst2
+
+print "Values:", [map(r.eval, lst)]
+print "Sorted:", sorted(lst, r.cmp)
+
+cofaces = []
+r.vertex_cofaces(2, 1, 3, cofaces.append)
+print "Cofaces of vertex 2:", cofaces
+
+cofaces = []
+r.vertex_cofaces(2, 1, 3, cofaces.append, [0,2,4])
+print "Cofaces of vertex 2 on vertices [0,2,4]:", cofaces
+
+f = Filtration(lst, r.cmp)
+p = StaticPersistence(f)
+p.pair_simplices()
+for s in p:
+    print s, s.sign()