Merged in Chris' changes dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Tue, 04 Aug 2009 17:12:47 -0700
branchdev
changeset 158 ce2aa05994f0
parent 149 3d15aca95dfb (diff)
parent 157 700cbac5b23c (current diff)
child 159 92318b22ffa5
Merged in Chris' changes
--- a/bindings/python/CMakeLists.txt	Tue Aug 04 16:58:44 2009 -0700
+++ b/bindings/python/CMakeLists.txt	Tue Aug 04 17:12:47 2009 -0700
@@ -11,8 +11,11 @@
                                                 chain.cpp
                                                 static-persistence.cpp
                                                 simplex.cpp
+                                                birthid.cpp
                                                 zigzag-persistence.cpp
+                                                cohomology-persistence.cpp
                                                 rips.cpp
+                                                distances.cpp
                             )
 set                         (bindings_libraries ${libraries})
 
@@ -26,6 +29,7 @@
     link_libraries          (${CGAL_LIBRARY})                                            
 else                            (CGAL_FOUND)
     message(STATUS "CGAL not found, alphashape bindings will not be built")
+    add_definitions         (-DNO_CGAL)
 endif                       (CGAL_FOUND)
 
 add_library                 (_dionysus SHARED   ${sources})
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/bindings/python/birthid.cpp	Tue Aug 04 17:12:47 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	Tue Aug 04 17:12:47 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	Tue Aug 04 17:12:47 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	Tue Aug 04 17:12:47 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	Tue Aug 04 16:58:44 2009 -0700
+++ b/bindings/python/dionysus.cpp	Tue Aug 04 17:12:47 2009 -0700
@@ -7,11 +7,17 @@
 void export_filtration();
 void export_static_persistence();
 void export_chain();
+void export_birthid();
 void export_zigzag_persistence();
+void export_cohomology_persistence();
 
 void export_rips();
+void export_pairwise_distances();
+
+#ifndef NO_CGAL
 void export_alphashapes2d();
 void export_alphashapes3d();
+#endif
 
 #ifdef LOGGING
 void            enable_log(std::string s)
@@ -27,11 +33,17 @@
     export_static_persistence();
     export_chain();
 
+    export_birthid();
     export_zigzag_persistence();
+    export_cohomology_persistence();
 
     export_rips();
+    export_pairwise_distances();
+
+#ifndef NO_CGAL
     export_alphashapes2d();
     export_alphashapes3d();
+#endif
 
 #ifdef LOGGING
     bp::def("enable_log",           &enable_log);
--- a/bindings/python/dionysus/__init__.py	Tue Aug 04 16:58:44 2009 -0700
+++ b/bindings/python/dionysus/__init__.py	Tue Aug 04 17:12:47 2009 -0700
@@ -1,5 +1,5 @@
 from    _dionysus   import *
-from    distances   import *
+from    distances   import l2, ExplicitDistances, points_file
 from    zigzag      import *
 
 
--- a/bindings/python/dionysus/distances.py	Tue Aug 04 16:58:44 2009 -0700
+++ b/bindings/python/dionysus/distances.py	Tue Aug 04 17:12:47 2009 -0700
@@ -1,9 +1,7 @@
 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
+    return sqrt(sum((x**2 for x in p)))
 
 # Pairwise distances between the elements of `points` with respect to some `norm`
 class PairwiseDistances:
@@ -37,5 +35,6 @@
 def points_file(filename):
     fd = open(filename)
     for line in fd.xreadlines():
+        if line.startswith('#'): continue
         yield map(float, line.strip().split())
     fd.close()
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/bindings/python/distances.cpp	Tue Aug 04 17:12:47 2009 -0700
@@ -0,0 +1,21 @@
+#include <boost/python.hpp>
+namespace bp = boost::python;
+
+#include "distances.h"
+namespace dp = dionysus::python;
+
+boost::shared_ptr<dp::ListPointPairwiseDistances>       init_from_list(bp::list lst)
+{
+    boost::shared_ptr<dp::ListPointPairwiseDistances>   p(new dp::ListPointPairwiseDistances(lst));
+    return p;
+}
+
+void export_pairwise_distances()
+{
+    bp::class_<dp::ListPointPairwiseDistances>("PairwiseDistances", bp::no_init)
+        .def("__init__",        bp::make_constructor(&init_from_list))
+        .def("__len__",         &dp::ListPointPairwiseDistances::size)
+        .def("__call__",        &dp::ListPointPairwiseDistances::operator())
+    ;
+}
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/bindings/python/distances.h	Tue Aug 04 17:12:47 2009 -0700
@@ -0,0 +1,54 @@
+#include <utilities/log.h>
+
+#include <boost/python.hpp>
+namespace bp = boost::python;
+
+namespace dionysus { 
+namespace python   {
+
+typedef     bp::list            ListPoint;
+
+struct ListPointL2Distance:
+    public std::binary_function<bp::object, bp::object, double>
+{
+    result_type     operator()(bp::object p1, bp::object p2) const
+    {
+        ListPoint lp1 = bp::extract<ListPoint>(p1), lp2 = bp::extract<ListPoint>(p2);
+
+        AssertMsg(bp::len(lp1) == bp::len(lp2), "Points must be in the same dimension (in L2Distance): dim1=%d, dim2=%d", bp::len(lp1), bp::len(lp2));
+        result_type sum = 0;
+        for (size_t i = 0; i < bp::len(lp1); ++i)
+        {
+            double diff = bp::extract<double>(lp1[i]) - bp::extract<double>(lp2[i]);
+            sum += diff*diff;
+        }
+
+        return sqrt(sum);
+    }
+};
+
+class ListPointPairwiseDistances
+{
+    public:
+        typedef             bp::list                                        Container;
+        typedef             ListPointL2Distance                             Distance;
+        typedef             unsigned                                        IndexType;
+        typedef             Distance::result_type                           DistanceType;
+
+
+                            ListPointPairwiseDistances(Container container): 
+                                container_(container)                       {}
+
+        DistanceType        operator()(IndexType a, IndexType b) const      { return distance_(container_[a], container_[b]); }
+
+        size_t              size() const                                    { return bp::len(container_); }
+        IndexType           begin() const                                   { return 0; }
+        IndexType           end() const                                     { return size(); }
+
+    private:
+        Container           container_;
+        Distance            distance_;
+};
+
+} }     // namespace dionysus::python
+
--- a/bindings/python/filtration.cpp	Tue Aug 04 16:58:44 2009 -0700
+++ b/bindings/python/filtration.cpp	Tue Aug 04 17:12:47 2009 -0700
@@ -7,7 +7,7 @@
 namespace bp = boost::python;
 
 
-#include "filtration.h"      // defines ListFiltration, ListTraits, ListRandomAccessIterator
+#include "filtration.h"      // defines ListFiltration, ListTraits
 #include "utils.h"           // defines PythonCmp
 namespace dp = dionysus::python;
 
--- a/bindings/python/filtration.h	Tue Aug 04 16:58:44 2009 -0700
+++ b/bindings/python/filtration.h	Tue Aug 04 17:12:47 2009 -0700
@@ -4,48 +4,19 @@
 #include <topology/filtration.h>
 #include <boost/python.hpp>
 #include "simplex.h"
+#include "utils.h"                      // for ListRandomAccessIterator
 
 namespace bp = boost::python;
 
 namespace dionysus { 
 namespace python   {
 
-// Random access iterator into python's list (using integer indices)
-class ListRandomAccessIterator:
-    public boost::iterator_adaptor<ListRandomAccessIterator,                // Derived
-                                   boost::counting_iterator<unsigned>,      // Base
-                                   SimplexObject,                           // Value
-                                   boost::use_default,
-                                   SimplexObject>
-{
-    public:
-        typedef                 ListRandomAccessIterator                                        Self;
-        typedef                 boost::iterator_adaptor<ListRandomAccessIterator,           
-                                                        boost::counting_iterator<unsigned>,     
-                                                        SimplexObject,
-                                                        boost::use_default,
-                                                        SimplexObject>                          Parent;
-                    
-                                ListRandomAccessIterator()                                      {}
-
-                                ListRandomAccessIterator(bp::list l, unsigned i):
-                                    Parent(i), l_(l)                                            {}
-
-    private:
-        friend class boost::iterator_core_access;
-        friend class FiltrationPythonIterator;
-
-        Parent::reference       dereference() const                                             { return bp::object(l_[*(this->base())]); }
-
-        bp::list                l_;
-};
-
 // ComplexTraits describing complexes of type list
 struct ListTraits
 {
     typedef     bp::list                                        Complex;
     typedef     SimplexObject                                   Simplex;
-    typedef     ListRandomAccessIterator                        Index;
+    typedef     ListRandomAccessIterator<Simplex>               Index;
     typedef     std::less<Index>                                IndexComparison;
 
     typedef     BinarySearchMap<SimplexVD, Index,
--- a/bindings/python/optional.h	Tue Aug 04 16:58:44 2009 -0700
+++ b/bindings/python/optional.h	Tue Aug 04 17:12:47 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/utils.h	Tue Aug 04 16:58:44 2009 -0700
+++ b/bindings/python/utils.h	Tue Aug 04 17:12:47 2009 -0700
@@ -7,6 +7,39 @@
 namespace dionysus {
 namespace python   {
 
+// Random access iterator into python's list (using integer indices)
+template<class Value>
+class ListRandomAccessIterator:
+    public boost::iterator_adaptor<ListRandomAccessIterator<Value>,         // Derived
+                                   boost::counting_iterator<unsigned>,      // Base
+                                   Value,                                   // Value
+                                   boost::use_default,
+                                   Value>
+{
+    public:
+        typedef                 ListRandomAccessIterator                                        Self;
+        typedef                 boost::iterator_adaptor<ListRandomAccessIterator,           
+                                                        boost::counting_iterator<unsigned>,     
+                                                        Value,
+                                                        boost::use_default,
+                                                        Value>                                  Parent;
+                    
+                                ListRandomAccessIterator()                                      {}
+
+                                ListRandomAccessIterator(bp::list l, unsigned i):
+                                    Parent(i), l_(l)                                            {}
+
+    private:
+        friend class boost::iterator_core_access;
+        friend class FiltrationPythonIterator;
+
+        typename Parent::reference       
+                                dereference() const                                             { return bp::object(l_[*(this->base())]); }
+
+        bp::list                l_;
+};
+
+// Adaptor of a Pyhon object to act as a C++-style comparison functor
 struct PythonCmp
 {
     template<class T>
@@ -14,7 +47,7 @@
 
                     PythonCmp(bp::object cmp): cmp_(cmp)    {}
 
-    bp::object cmp_;
+    bp::object      cmp_;
 };
 
 } } // namespace dionysus::python
--- a/bindings/python/zigzag-persistence.cpp	Tue Aug 04 16:58:44 2009 -0700
+++ b/bindings/python/zigzag-persistence.cpp	Tue Aug 04 17:12:47 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	Tue Aug 04 16:58:44 2009 -0700
+++ b/bindings/python/zigzag-persistence.h	Tue Aug 04 17:12:47 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	Tue Aug 04 16:58:44 2009 -0700
+++ b/doc/examples/cohomology.rst	Tue Aug 04 17:12:47 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	Tue Aug 04 17:12:47 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	Tue Aug 04 16:58:44 2009 -0700
+++ b/doc/python/overview.rst	Tue Aug 04 17:12:47 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/rips.rst	Tue Aug 04 16:58:44 2009 -0700
+++ b/doc/python/rips.rst	Tue Aug 04 17:12:47 2009 -0700
@@ -60,9 +60,11 @@
         def __call__(self, x, y):
             return math.fabs(y-x)
 
-The bindings provide a pure Python class :class:`PairwiseDistances` to deal with
-explicit points in a Euclidean space. It is defined in
-:sfile:`bindings/python/dionysus/distances.py`::
+The bindings expose a C++ class as a Python class :class:`PairwiseDistances` to deal with
+explicit points in a Euclidean space. In pure Python it could be defined as
+follows (in fact it used to be a pure Python class, and one may still find it in 
+:sfile:`bindings/python/dionysus/distances.py`; its performance is much slower
+than its pure C++ analog)::
 
     class PairwiseDistances:
         def __init__(self, points, norm = l2):
@@ -84,6 +86,9 @@
     distances = PairwiseDistances(points)
     distances = ExplicitDistances(distances)
 
+With :class:`PairwiseDistances` being a C++ class, and
+:class:`ExplicitDistances` being pure Python, the speead-up seems minor.
+
 
 Example
 -------
--- a/doc/python/zigzag-persistence.rst	Tue Aug 04 16:58:44 2009 -0700
+++ b/doc/python/zigzag-persistence.rst	Tue Aug 04 17:12:47 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/doc/tutorial.rst	Tue Aug 04 16:58:44 2009 -0700
+++ b/doc/tutorial.rst	Tue Aug 04 17:12:47 2009 -0700
@@ -70,7 +70,7 @@
 the entire filtration at once, compute its persistence in one operation, and
 then examine the pairings. The second way is to feed simplices one by one in an
 *online* manner and manually keep track of the pairs that are created.
-:class:`ZigzagPersistence` and `CohomologyPersistence` |cpp-only| accept their
+:class:`ZigzagPersistence` and :class:`CohomologyPersistence` accept their
 input this way,
 
 
@@ -123,7 +123,10 @@
 ^^^^^^
 
 Class :class:`ZigzagPersistence` accepts additions and removals of the simplices
-one by one, and returns an internal representation of the simplices. When one
+one by one, and returns an internal representation of the simplices.
+(:class:`CohomologyPersistence` works in the same way, but accepts only
+additions of the simplices.) 
+When one
 adds a simplex via :meth:`ZigzagPersistence.add`, one must provide its boundary
 in this internal representation together with a *birth value* which
 :class:`ZigzagPersistence` will store in case a class is born as a result of the
@@ -161,7 +164,9 @@
 :func:`add_simplices` and :func:`remove_simplices`.
 
 See the :ref:`bindings reference <python-bindings>` for more details, and
-:ref:`triangle-zigzag-example` for an example.
+:ref:`triangle-zigzag-example` for an example of :class:`ZigzagPersistence`.
+See :ref:`rips-pairwise-cohomology` for an example of
+:class:`CohomologyPersistence`.
 
 
 .. _speed-up-suggestions:
@@ -172,7 +177,9 @@
 Currently, when the choice comes between efficiency and flexibility, the Python
 bindings err on the side of flexibility. There is hope that in the future the
 choice won't really be necessary. Meanwhile, one can use a few techniques that
-speed up computation at the expense of memory:
+speed up computation at the expense of memory. Note, however, that since the
+recent switch of :class:`PairwiseDistances` to C++ rather than pure Python, it
+is not clear whether these deliver a substantial speed-up:
 
 * To avoid (possibly expensive) computation of distances during Rips complex
   generation, store :class:`ExplicitDistances` (see :ref:`distances`)::
--- a/examples/cohomology/cocycle.py	Tue Aug 04 16:58:44 2009 -0700
+++ b/examples/cohomology/cocycle.py	Tue Aug 04 17:12:47 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	Tue Aug 04 17:12:47 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)
--- a/examples/consistency/rips-consistency-zigzag.cpp	Tue Aug 04 16:58:44 2009 -0700
+++ b/examples/consistency/rips-consistency-zigzag.cpp	Tue Aug 04 17:12:47 2009 -0700
@@ -68,7 +68,7 @@
 { 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        report_death(std::ostream& out, const BirthInfo& birth, const BirthInfo& death, unsigned skeleton_dimension);
 void        make_boundary(const Smplx& s, Complex& c, const Zigzag& zz, Boundary& b);
 void        process_command_line_options(int           argc,
                                          char*         argv[],
@@ -78,9 +78,9 @@
                                          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);
+                        Complex& complex, Zigzag& zz, std::ostream& out, Timer& add, unsigned skeleton_dimension);
 void        remove_simplex(const Smplx& s, const BirthInfo& birth, const BirthInfo& death, 
-                           Complex& complex, Zigzag& zz, std::ostream& out, Timer& remove);
+                           Complex& complex, Zigzag& zz, std::ostream& out, Timer& remove, unsigned skeleton_dimension);
 
 int main(int argc, char* argv[])
 {
@@ -179,7 +179,7 @@
         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);
+                        complex, zz, out, add, skeleton_dimension);
 
             // Record cofaces with all other vertices outside of the subcomplex
             coface.start();
@@ -193,7 +193,7 @@
         // 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);
+                        complex, zz, out, add, skeleton_dimension);
         rInfo("  Cross simplices added");
 
 
@@ -213,12 +213,12 @@
 
         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);
+                           complex, zz, out, remove, skeleton_dimension);
         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);
+                           complex, zz, out, remove, skeleton_dimension);
         rInfo("  Subcomplex simplices removed");
         
         Dimension betti_1 = 0;
@@ -249,14 +249,14 @@
     }
 }
 
-void        report_death(std::ostream& out, const BirthInfo& birth, const BirthInfo& death)
+void        report_death(std::ostream& out, const BirthInfo& birth, const BirthInfo& death, unsigned skeleton_dimension)
 {
-    if (death > birth)
+    if (birth.dimension < skeleton_dimension && 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)
+                        Complex& complex, Zigzag& zz, std::ostream& out, Timer& add, unsigned skeleton_dimension)
 {
     rDebug("Adding simplex: %s", tostring(s).c_str());
 
@@ -264,13 +264,13 @@
     make_boundary(s, complex, zz, b);
     add.start();
     boost::tie(idx, d)  = zz.add(b, birth);
-    if (d) report_death(out, *d, death);
+    if (d) report_death(out, *d, death, skeleton_dimension);
     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)
+                           Complex& complex, Zigzag& zz, std::ostream& out, Timer& remove, unsigned skeleton_dimension)
 {
     rDebug("Removing simplex: %s", tostring(s).c_str());
 
@@ -278,7 +278,7 @@
     remove.start();
     Death d = zz.remove(si->second, birth);
     remove.stop();
-    if (d) report_death(out, *d, death);
+    if (d) report_death(out, *d, death, skeleton_dimension);
     complex.erase(si);
 }
 
--- a/examples/rips/rips-pairwise.py	Tue Aug 04 16:58:44 2009 -0700
+++ b/examples/rips/rips-pairwise.py	Tue Aug 04 17:12:47 2009 -0700
@@ -6,49 +6,52 @@
 from    sys         import argv, exit
 import  time
 
+def main(filename, skeleton, max):
+    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"
 
-if len(argv) < 4:
-    print "Usage: %s POINTS SKELETON MAX" % argv[0]
-    exit()
-
-filename = argv[1]
-skeleton = int(argv[2])
-max = float(argv[3])
+    simplices = []
+    rips.generate(skeleton, max, simplices.append)
+    print time.asctime(), "Generated complex: %d simplices" % len(simplices)
 
-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"
+    # 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 = []
-rips.generate(skeleton, max, simplices.append)
-print time.asctime(), "Generated complex: %d simplices" % len(simplices)
+    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"
 
-# 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]
+    print "Outputting persistence diagram"
+    for i in p:
+        if i.sign():
+            b = simplices[f[p(i)]]
 
-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"
+            if b.dimension() >= skeleton: continue
 
-p = StaticPersistence(f)
-print time.asctime(), "Initialized StaticPersistence"
-
-p.pair_simplices()
-print time.asctime(), "Simplices paired"
+            if i == i.pair:
+                print b.dimension(), b.data, "inf"
+                continue
 
-print "Outputting persistence diagram"
-for i in p:
-    if i.sign():
-        b = simplices[f[p(i)]]
-
-        if b.dimension() >= skeleton: continue
+            d = simplices[f[p(i.pair)]]
+            print b.dimension(), b.data, d.data
 
-        if i == i.pair:
-            print b.dimension(), b.data, "inf"
-            continue
+if __name__ == '__main__':
+    if len(argv) < 4:
+        print "Usage: %s POINTS SKELETON MAX" % argv[0]
+        exit()
 
-        d = simplices[f[p(i.pair)]]
-        print b.dimension(), b.data, d.data
+    filename = argv[1]
+    skeleton = int(argv[2])
+    max = float(argv[3])
+
+    main(filename, skeleton, max)
--- a/examples/triangle/triangle-zigzag.py	Tue Aug 04 16:58:44 2009 -0700
+++ b/examples/triangle/triangle-zigzag.py	Tue Aug 04 17:12:47 2009 -0700
@@ -30,6 +30,6 @@
 for s in sorted(complex.keys(), reverse = True):
     print "%d: Removing %s" % (b, s)
     d = zz.remove(complex[s], b)
-    complex[s] = None
+    del complex[s]
     if d:   print "Interval (%d, %d)" % (d, b-1)
     b += 1
--- a/include/geometry/distances.h	Tue Aug 04 16:58:44 2009 -0700
+++ b/include/geometry/distances.h	Tue Aug 04 17:12:47 2009 -0700
@@ -1,6 +1,8 @@
 #ifndef __DISTANCES_H__
 #define __DISTANCES_H__
 
+#include <vector>
+
 /**
  * Class: ExplicitDistances 
  * Stores the pairwise distances of Distances_ instance passed at construction.