Added Python bindings for CohomologyPersistence (+ example + documentation) dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Thu, 30 Jul 2009 10:23:31 -0700
branchdev
changeset 146 4e27f1f7c169
parent 145 ee096f207dfb
child 147 d39a20acb253
Added Python bindings for CohomologyPersistence (+ example + documentation)
bindings/python/CMakeLists.txt
bindings/python/birthid.cpp
bindings/python/birthid.h
bindings/python/cohomology-persistence.cpp
bindings/python/cohomology-persistence.h
bindings/python/dionysus.cpp
bindings/python/optional.h
bindings/python/zigzag-persistence.cpp
bindings/python/zigzag-persistence.h
doc/examples/cohomology.rst
doc/python/cohomology-persistence.rst
doc/python/overview.rst
doc/python/zigzag-persistence.rst
examples/cohomology/cocycle.py
examples/cohomology/rips-pairwise-cohomology.py
--- a/bindings/python/CMakeLists.txt	Fri Jul 24 14:18:16 2009 -0700
+++ b/bindings/python/CMakeLists.txt	Thu Jul 30 10:23:31 2009 -0700
@@ -11,7 +11,9 @@
                                                 chain.cpp
                                                 static-persistence.cpp
                                                 simplex.cpp
+                                                birthid.cpp
                                                 zigzag-persistence.cpp
+                                                cohomology-persistence.cpp
                                                 rips.cpp
                             )
 set                         (bindings_libraries ${libraries})
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/bindings/python/birthid.cpp	Thu Jul 30 10:23:31 2009 -0700
@@ -0,0 +1,10 @@
+#include "birthid.h"
+#include "optional.h"
+
+namespace dp = dionysus::python;
+
+void export_birthid()
+{
+    python_optional<dp::BirthID>();   
+}
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/bindings/python/birthid.h	Thu Jul 30 10:23:31 2009 -0700
@@ -0,0 +1,15 @@
+#ifndef __PYTHON_BIRTHID_H__
+#define __PYTHON_BIRTHID_H__
+
+#include <boost/python.hpp>
+
+namespace dionysus {
+namespace python   {
+
+//typedef         int                             BirthID;
+//typedef         boost::python::long_            BirthID;
+typedef         boost::python::object           BirthID;
+
+} } // namespace dionysus::python
+
+#endif // __PYTHON_BIRTHID_H__
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/bindings/python/cohomology-persistence.cpp	Thu Jul 30 10:23:31 2009 -0700
@@ -0,0 +1,71 @@
+#include <topology/cohomology-persistence.h>
+
+#include <boost/python.hpp>
+#include <boost/python/stl_iterator.hpp>
+#include <boost/shared_ptr.hpp>
+namespace bp = boost::python;
+
+#include "cohomology-persistence.h"             // defines CohomPersistence
+#include "optional.h"
+namespace dp = dionysus::python;
+
+
+// CohomPersistence
+boost::shared_ptr<dp::CohomPersistence>     init_from_prime(unsigned p = 11)
+{
+    dp::CohomPersistence::Field field(p);       // Zp
+
+    boost::shared_ptr<dp::CohomPersistence> chp(new dp::CohomPersistence(field));
+    return chp;
+}
+
+bp::tuple                                   chp_add(dp::CohomPersistence& chp, bp::object bdry, dp::BirthID birth)
+{
+    dp::CohomPersistence::SimplexIndex      i;
+    dp::CohomPersistence::Death             d;
+    boost::tie(i,d)                                 = chp.add(bp::stl_input_iterator<dp::CohomPersistence::SimplexIndex>(bdry),
+                                                              bp::stl_input_iterator<dp::CohomPersistence::SimplexIndex>(),
+                                                              birth); 
+    return bp::make_tuple(i,d);
+}
+
+dp::CohomPersistence::ZColumn::const_iterator     
+cocycle_zcolumn_begin(dp::CohomPersistence::Cocycle& ccl)                   
+{ return ccl.zcolumn.begin(); }
+
+dp::CohomPersistence::ZColumn::const_iterator     
+cocycle_zcolumn_end(dp::CohomPersistence::Cocycle& ccl)                   
+{ return ccl.zcolumn.end(); }
+
+
+// SimplexIndex
+template<class T>
+unsigned                            si_order(T& si)
+{
+    return si->order;
+}
+
+
+void export_cohomology_persistence()
+{
+    bp::class_<dp::CohomPersistence::SimplexIndex>("CHSimplexIndex")
+        .add_property("order",          &si_order<dp::CohomPersistence::SimplexIndex>)
+    ;
+    
+    bp::class_<dp::CohomPersistence::SNode>("CHSNode", bp::no_init)
+        .add_property("coefficient",    &dp::CohomPersistence::SNode::coefficient)
+        .add_property("si",             &dp::CohomPersistence::SNode::si)
+    ;
+
+    bp::class_<dp::CohomPersistence>("CohomologyPersistence")
+        .def("__init__",        bp::make_constructor(&init_from_prime))
+        .def("add",             &chp_add)
+        
+        .def("__iter__",        bp::range(&dp::CohomPersistence::begin, &dp::CohomPersistence::end))
+    ;
+
+    bp::class_<dp::CohomPersistence::Cocycle>("Cocycle", bp::no_init)
+        .add_property("birth",  &dp::CohomPersistence::Cocycle::birth)
+        .def("__iter__",        bp::range(&cocycle_zcolumn_begin, &cocycle_zcolumn_end))
+    ;
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/bindings/python/cohomology-persistence.h	Thu Jul 30 10:23:31 2009 -0700
@@ -0,0 +1,16 @@
+#ifndef __PYTHON_ZIGZAG_PERSISTENCE_H__
+#define __PYTHON_ZIGZAG_PERSISTENCE_H__
+
+#include <topology/cohomology-persistence.h>
+#include <boost/python.hpp>
+
+#include "birthid.h"
+
+namespace dionysus {
+namespace python   {
+
+typedef         CohomologyPersistence<BirthID>  CohomPersistence;
+
+} } // namespace dionysus::python
+
+#endif
--- a/bindings/python/dionysus.cpp	Fri Jul 24 14:18:16 2009 -0700
+++ b/bindings/python/dionysus.cpp	Thu Jul 30 10:23:31 2009 -0700
@@ -7,7 +7,9 @@
 void export_filtration();
 void export_static_persistence();
 void export_chain();
+void export_birthid();
 void export_zigzag_persistence();
+void export_cohomology_persistence();
 
 void export_rips();
 #ifndef NO_CGAL
@@ -29,7 +31,9 @@
     export_static_persistence();
     export_chain();
 
+    export_birthid();
     export_zigzag_persistence();
+    export_cohomology_persistence();
 
     export_rips();
 #ifndef NO_CGAL
--- a/bindings/python/optional.h	Fri Jul 24 14:18:16 2009 -0700
+++ b/bindings/python/optional.h	Thu Jul 30 10:23:31 2009 -0700
@@ -2,6 +2,8 @@
 #define __PYTHON_OPTIONAL_H__
 
 #include <boost/python.hpp>
+#include <boost/optional.hpp>
+#include <boost/utility.hpp>
 
 // Taken from an email by John Wiegley; 
 // http://mail.python.org/pipermail/cplusplus-sig/2007-May/012003.html
--- a/bindings/python/zigzag-persistence.cpp	Fri Jul 24 14:18:16 2009 -0700
+++ b/bindings/python/zigzag-persistence.cpp	Thu Jul 30 10:23:31 2009 -0700
@@ -67,12 +67,10 @@
 
 void export_zigzag_persistence()
 {
-    python_optional<dp::BirthID>();   
-
-    bp::class_<dp::ZZPersistence::SimplexIndex>("SimplexIndex")
+    bp::class_<dp::ZZPersistence::SimplexIndex>("ZZSimplexIndex")
         .def("order",           &si_order<dp::ZZPersistence::SimplexIndex>);
     
-    bp::class_<dp::IZZPersistence::SimplexIndex>("ISimplexIndex")
+    bp::class_<dp::IZZPersistence::SimplexIndex>("IZZSimplexIndex")
         .def("order",           &si_order<dp::IZZPersistence::SimplexIndex>);
 
     bp::class_<dp::ZZPersistence>("ZigzagPersistence")
--- a/bindings/python/zigzag-persistence.h	Fri Jul 24 14:18:16 2009 -0700
+++ b/bindings/python/zigzag-persistence.h	Thu Jul 30 10:23:31 2009 -0700
@@ -5,13 +5,11 @@
 #include <topology/image-zigzag-persistence.h>
 #include <boost/python.hpp>
 
+#include "birthid.h"
+
 namespace dionysus {
 namespace python   {
 
-//typedef         int                             BirthID;
-//typedef         boost::python::long_            BirthID;
-typedef         boost::python::object           BirthID;
-
 typedef         ZigzagPersistence<BirthID>      ZZPersistence;
 typedef         ImageZigzagPersistence<BirthID> IZZPersistence;
 
--- a/doc/examples/cohomology.rst	Fri Jul 24 14:18:16 2009 -0700
+++ b/doc/examples/cohomology.rst	Thu Jul 30 10:23:31 2009 -0700
@@ -85,3 +85,45 @@
 .. _CVXOPT:                                 http://abel.ee.ucla.edu/cvxopt/
 .. _`Stanford MATLAB implementation`:       http://www.stanford.edu/group/SOL/software/lsqr.html
 .. _`Jeffery Kline`:                        http://pages.cs.wisc.edu/~kline/
+
+
+.. _rips-pairwise-cohomology:
+
+Python cohomology computation
+-----------------------------
+
+:sfile:`examples/cohomology/rips-pairwise-cohomology.py` gives an example of the
+same computation performed in Python (but with the output in a different format).
+
+After the simplicial complex is computed in a list `simplices`, and the list is
+sorted with respect to the Rips filtration order, the simplices are inserted
+into the :class:`CohomologyPersistence` one by one::
+
+    # list simplices is created
+
+    ch = CohomologyPersistence(prime)
+    complex = {}
+
+    for s in simplices:
+        i,d = ch.add([complex[sb] for sb in s.boundary], (s.dimension(), s.data))
+        complex[s] = i
+        if d: 
+            dimension, birth = d
+            print dimension, birth, s.data
+        # else birth
+
+Above dictionary `complex` maintains the map of simplices to indices returned by
+:meth:`CohomologyPersistence.add`.  The pair `(dimension, data)` is used as the
+birth value. Here `data` is the value associated with the simplex in the Rips
+filtration. The pair is returned back if a death occurs, and is printed on the
+standard output. After the for loop finishes, one may output infinite
+persistence classes with the following for loop::
+
+    for ccl in ch:
+        dimension, birth = ccl.birth
+        if dimension >= skeleton: continue
+        print dimension, birth, 'inf'         # dimension, simplex data = birth
+    
+Naturally one may iterate over `ccl` which is of type :class:`Cocycle` and
+extract more information. For example, this is necessary to get the coefficients
+that serve as the input for :sfile:`examples/cohomology/cocycle.py`.
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/doc/python/cohomology-persistence.rst	Thu Jul 30 10:23:31 2009 -0700
@@ -0,0 +1,68 @@
+:class:`CohomologyPersistence` class
+====================================
+
+The :ref:`rips-pairwise-cohomology` illustrates the use of :class:`CohomologyPersistence`.
+
+.. class:: CohomologyPersistence
+
+    .. method:: __init__(prime)
+
+        Initializes :class:`CohomologyPersistence` with the given `prime`; from
+        this point on all the computation will be performed with coefficients
+        coming from :math:`\mathbb{Z}/prime \mathbb{Z}`.
+
+    .. method:: add(boundary, birth)
+        
+        Adds a simplex with the given `boundary` to the complex, i.e. 
+        :math:`K_{i+1} = K_i \cup \sigma` and `boundary` = :math:`\partial \sigma`.
+        If a new class is born as a result of the addition, `birth` is stored with 
+        it for future reference.
+
+        :returns: a pair (`i`, `d`). The first element is the index `i`. 
+                  It is the internal representation of the newly added simplex,
+                  and should be used later for removal or when constructing the
+                  boundaries of its cofaces. In other words, `boundary` must
+                  consist of these indices.  The second element `d` is the death
+                  element. It is `None` if a birth occurred, otherwise it
+                  contains the value passed as `birth` to
+                  :meth:`~CohomologyPersistence.add` when the class that just
+                  died was born.
+
+    .. method:: __iter__()
+
+        Iterator over the live cocycles stored in
+        :class:`CohomologyPersistence`. The returned elements are of the type
+        :class:`Cocycle` below.
+
+       
+.. class:: Cocycle        
+
+    .. attribute:: birth
+
+        The birth value associated with the cocycle. It is passed to
+        :class:`CohomologyPersistence` in method
+        :meth:`~CohomologyPersistence.add`.
+
+    .. method:: __iter__()
+
+        Iterator over the individual nodes (simplices) of the cocycle, each of type
+        :class:`CHSNode`.
+
+.. class:: CHSNode
+
+    .. attribute:: si
+
+        The index of the simplex, of type :class:`CHSimplexIndex`.
+
+    .. attribute:: coefficient
+
+        Coefficient in :math:`\mathbb{Z}/prime \mathbb{Z}` associated with the
+        simplex.
+
+
+.. class:: CHSimplexIndex
+
+    .. attribute:: order
+
+        The count associated with the simplex when it is inserted into
+        :class:`CohomologyPersistence`.
--- a/doc/python/overview.rst	Fri Jul 24 14:18:16 2009 -0700
+++ b/doc/python/overview.rst	Thu Jul 30 10:23:31 2009 -0700
@@ -17,6 +17,7 @@
     simplex.rst
     filtration.rst
     static-persistence.rst
+    cohomology-persistence.rst
     alphashapes.rst
     rips.rst
     zigzag-persistence.rst
--- a/doc/python/zigzag-persistence.rst	Fri Jul 24 14:18:16 2009 -0700
+++ b/doc/python/zigzag-persistence.rst	Thu Jul 30 10:23:31 2009 -0700
@@ -1,7 +1,8 @@
 :class:`ZigzagPersistence` class
 ================================
 
-The class deals with the setting :math:`K_1 \rightarrow K_2 \leftarrow K_3 \rightarrow \dots`
+The class deals with the setting :math:`K_1 \rightarrow K_2 \leftarrow K_3 \rightarrow \dots`.
+The :ref:`triangle-zigzag-example` illustrates the use of :class:`ZigzagPersistence`.
 
 .. class:: ZigzagPersistence
 
@@ -34,8 +35,6 @@
                   :meth:`~ZigzagPersistence.remove` when the class that just
                   died was born.
 
-The :ref:`triangle-zigzag-example` illustrates the use of :class:`ZigzagPersistence`.
-
 
 Auxilliary functions
 --------------------
--- a/examples/cohomology/cocycle.py	Fri Jul 24 14:18:16 2009 -0700
+++ b/examples/cohomology/cocycle.py	Thu Jul 30 10:23:31 2009 -0700
@@ -6,25 +6,7 @@
 from    sys             import argv, exit
 import  os.path
 
-def smooth(boundary_filename, cocycle_filename, vertexmap_filename):
-    boundary_list = []
-    with open(boundary_filename) as fp:
-        for line in fp.xreadlines():
-            if line.startswith('#'): continue
-            boundary_list.append(map(int, line.split()))
-
-    cocycle_list = []
-    with open(cocycle_filename) as fp:
-        for line in fp.xreadlines():
-            if line.startswith('#'): continue
-            cocycle_list.append(map(int, line.split()))
-
-    vertices = []
-    with open(vertexmap_filename) as fp:
-        for line in fp.xreadlines():
-            if line.startswith('#'): continue
-            vertices.append(map(int, line.split())[1])
-
+def smooth(boundary_list, cocycle_list, vertexmap_list):
     dimension = max((max(d[1], d[2]) for d in boundary_list))
     dimension += 1
 
@@ -69,6 +51,15 @@
     return values
 
 
+def read_list_file(filename):
+    list = []
+    with open(filename) as fp:
+        for line in fp.xreadlines():
+            if line.startswith('#'): continue
+            list.append(map(int, line.split()))
+    return list
+
+
 if __name__ == '__main__':
     if len(argv) < 4:
         print "Usage: %s BOUNDARY COCYCLE VERTEXMAP" % argv[0]
@@ -77,7 +68,12 @@
     boundary_filename = argv[1]
     cocycle_filename = argv[2]
     vertexmap_filename = argv[3]
-    values = smooth(boundary_filename, cocycle_filename, vertexmap_filename)
+
+    boundary_list = read_list_file(boundary_filename)
+    cocycle_list = read_list_file(cocycle_filename)
+    vertexmap_list = read_list_file(vertexmap_filename)
+
+    values = smooth(boundary_list, cocycle_list, vertexmap_list)
 
     outfn = os.path.splitext(cocycle_filename)[0] + '.val'
     with open(outfn, 'w') as f:
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/cohomology/rips-pairwise-cohomology.py	Thu Jul 30 10:23:31 2009 -0700
@@ -0,0 +1,61 @@
+#!/usr/bin/env python
+
+from    dionysus        import Simplex, CohomologyPersistence, points_file, PairwiseDistances, ExplicitDistances, Rips
+from    sys             import argv, exit
+import  time
+
+def main(filename, skeleton, max, prime = 11):
+    points = [p for p in points_file(filename)]
+    print '#', time.asctime(), "Points read"
+    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]
+
+    simplices.sort(data_dim_cmp)
+    print '#', time.asctime(), "Simplices sorted"
+
+    ch = CohomologyPersistence(prime)
+    complex = {}
+
+    for s in simplices:
+        i,d = ch.add([complex[sb] for sb in s.boundary], (s.dimension(), s.data))
+        complex[s] = i
+        if d: 
+            dimension, birth = d
+            print dimension, birth, s.data
+        # else birth
+
+    for ccl in ch:
+        dimension, birth = ccl.birth
+        if dimension >= skeleton: continue
+        print dimension, birth, 'inf'         # dimension, simplex data = birth
+        print "# Cocycle at (dim=%d, birth=%f)" % ccl.birth
+        for e in ccl:
+            print "#  ", e.si.order(), normalized(e.coefficient, prime)
+
+def normalized(coefficient, prime):
+    if coefficient > prime/2:
+        return coefficient - prime
+    return coefficient
+
+if __name__ == '__main__':
+    if len(argv) < 4:
+        print "Usage: %s POINTS SKELETON MAX [PRIME=11]" % argv[0]
+        exit()
+
+    filename = argv[1]
+    skeleton = int(argv[2])
+    max = float(argv[3])
+    prime = (len(argv) > 4 and argv[4]) or 11
+
+    main(filename, skeleton, max, prime)