Added arbitrary coefficients in boundaries (to C++ and Python) + cleaned up CohomologyPersistence.add overloading in Python dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Fri, 06 Nov 2009 14:19:08 -0800
branchdev
changeset 173 5fd3f43e6fbf
parent 172 a6605dc232f2
child 174 3f1034dca432
Added arbitrary coefficients in boundaries (to C++ and Python) + cleaned up CohomologyPersistence.add overloading in Python
bindings/python/cohomology-persistence.cpp
doc/python/cohomology-persistence.rst
examples/cohomology/rips-pairwise-cohomology.py
include/topology/cohomology-persistence.h
include/topology/cohomology-persistence.hpp
include/topology/field-arithmetic.h
--- a/bindings/python/cohomology-persistence.cpp	Wed Oct 28 14:38:36 2009 -0700
+++ b/bindings/python/cohomology-persistence.cpp	Fri Nov 06 14:19:08 2009 -0800
@@ -2,6 +2,7 @@
 
 #include <boost/python.hpp>
 #include <boost/python/stl_iterator.hpp>
+#include <boost/python/overloads.hpp>
 #include <boost/shared_ptr.hpp>
 namespace bp = boost::python;
 
@@ -19,36 +20,33 @@
     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);
-}
 
-
-bp::tuple                                   chp_add_store(dp::CohomPersistence& chp, bp::object bdry, dp::BirthID birth, bool store)
+bp::tuple                                   chp_add(dp::CohomPersistence& chp, 
+                                                    bp::object bdry, 
+                                                    dp::BirthID birth, 
+                                                    bool store, 
+                                                    bool image, 
+                                                    bp::object coefficients)
 {
     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, store); 
+
+    if (coefficients)
+    {
+        boost::tie(i,d)                         = chp.add(bp::stl_input_iterator<int>(coefficients),
+                                                          bp::stl_input_iterator<dp::CohomPersistence::SimplexIndex>(bdry),
+                                                          bp::stl_input_iterator<dp::CohomPersistence::SimplexIndex>(),
+                                                          birth, store, dp::CohomPersistence::SimplexData(), image); 
+    } else
+    {
+        boost::tie(i,d)                         = chp.add(bp::stl_input_iterator<dp::CohomPersistence::SimplexIndex>(bdry),
+                                                          bp::stl_input_iterator<dp::CohomPersistence::SimplexIndex>(),
+                                                          birth, store, dp::CohomPersistence::SimplexData(), image); 
+    }
+
     return bp::make_tuple(i,d);
 }
 
-bp::tuple                                   chp_add_store_image(dp::CohomPersistence& chp, bp::object bdry, dp::BirthID birth, bool store, bool image)
-{
-    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, store, dp::CohomPersistence::SimplexData(), image); 
-    return bp::make_tuple(i,d);
-}
 
 dp::CohomPersistence::ZColumn::const_iterator     
 cocycle_zcolumn_begin(dp::CohomPersistence::Cocycle& ccl)                   
@@ -67,6 +65,8 @@
 }
 
 
+BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(add_overloads, add, 2, 4)
+
 void export_cohomology_persistence()
 {
     bp::class_<dp::CohomPersistence::SimplexIndex>("CHSimplexIndex")
@@ -80,9 +80,7 @@
 
     bp::class_<dp::CohomPersistence>("CohomologyPersistence")
         .def("__init__",        bp::make_constructor(&init_from_prime))
-        .def("add",             &chp_add)
-        .def("add",             &chp_add_store)
-        .def("add",             &chp_add_store_image)
+        .def("add",             &chp_add, (bp::arg("bdry"), bp::arg("birth"), bp::arg("store")=true, bp::arg("image")=true, bp::arg("coefficients")=false))
         
         .def("__iter__",        bp::range(&dp::CohomPersistence::begin, &dp::CohomPersistence::end))
     ;
--- a/doc/python/cohomology-persistence.rst	Wed Oct 28 14:38:36 2009 -0700
+++ b/doc/python/cohomology-persistence.rst	Fri Nov 06 14:19:08 2009 -0800
@@ -11,7 +11,7 @@
         this point on all the computation will be performed with coefficients
         coming from :math:`\mathbb{Z}/prime \mathbb{Z}`.
 
-    .. method:: add(boundary, birth, [store = True, [image = True]])
+    .. method:: add(boundary, birth, [store = True], [image = True], [coefficients = []])
         
         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`.
@@ -30,6 +30,10 @@
         by the inclusion of :math:`L` in :math:`K`). `image` indicates whether
         the simplex added belongs to :math:`L` or not.
 
+        If given, `coefficients` is a list parallel to `boundary` that provides
+        coefficients for the corresponding boundary elements. If empty, it is
+        assumed to be :math:`(-1)^i`.
+
         :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
--- a/examples/cohomology/rips-pairwise-cohomology.py	Wed Oct 28 14:38:36 2009 -0700
+++ b/examples/cohomology/rips-pairwise-cohomology.py	Fri Nov 06 14:19:08 2009 -0800
@@ -19,7 +19,7 @@
     # 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 '#', time.asctime(), simplices[0], '...', simplices[-1]
 
     simplices.sort(data_dim_cmp)
     print '#', time.asctime(), "Simplices sorted"
@@ -28,7 +28,7 @@
     complex = {}
 
     for s in simplices:
-        i,d = ch.add([complex[sb] for sb in s.boundary], (s.dimension(), s.data))
+        i,d = ch.add([complex[sb] for sb in s.boundary], (s.dimension(), s.data), store = (s.dimension() < skeleton))
         complex[s] = i
         if d: 
             dimension, birth = d
--- a/include/topology/cohomology-persistence.h	Wed Oct 28 14:38:36 2009 -0700
+++ b/include/topology/cohomology-persistence.h	Fri Nov 06 14:19:08 2009 -0800
@@ -61,9 +61,13 @@
         // BI = BoundaryIterator; it should dereference to a SimplexIndex
         template<class BI>
         IndexDeathPair      add(BI begin, BI end, BirthInfo b, bool store = true, const SimplexData& sd = SimplexData(), bool image = true);
+        
+        // if sign needs to be specified explicitly, provide (parallel) coefficient_iter
+        template<class BI, class CI>
+        IndexDeathPair      add(CI coefficient_iter, BI begin, BI end, BirthInfo b, bool store = true, const SimplexData& sd = SimplexData(), bool image = true);
 
         void                show_cocycles() const;
-        CocycleIndex        begin()                                                     { return cocycles_.begin(); }
+        CocycleIndex        begin()                                                     { return image_begin_; }
         CocycleIndex        end()                                                       { return cocycles_.end(); }
 
     private:
--- a/include/topology/cohomology-persistence.hpp	Wed Oct 28 14:38:36 2009 -0700
+++ b/include/topology/cohomology-persistence.hpp	Fri Nov 06 14:19:08 2009 -0800
@@ -6,6 +6,14 @@
 #include <utilities/log.h>
 #include <utilities/indirect.h>
 
+#include <boost/iterator/transform_iterator.hpp>
+#include <boost/iterator/counting_iterator.hpp>
+#include <boost/foreach.hpp>
+
+#include <boost/lambda/lambda.hpp>
+#include <boost/lambda/bind.hpp>
+namespace bl = boost::lambda;
+
 #ifdef LOGGING
 static rlog::RLogChannel* rlCohomology =                DEF_CHANNEL("topology/cohomology",        rlog::Log_Debug);
 #endif
@@ -16,6 +24,14 @@
     public:
         bool        operator()(const SNode& s1, const SNode& s2) const                  { return s1.si->order < s2.si->order; }
 };
+    
+
+struct Alternator
+{
+    typedef         int         result_type;
+    int             operator()(int x) const                                             { return (x % 2)*2 - 1; }
+};
+    
 
 template<class BirthInfo, class SimplexData, class Field>
 template<class BI>
@@ -23,6 +39,20 @@
 CohomologyPersistence<BirthInfo, SimplexData, Field>::
 add(BI begin, BI end, BirthInfo birth, bool store, const SimplexData& sd, bool image)
 {
+    // Set coefficient to be an iterator over (-1)^i
+    
+    return add(boost::make_transform_iterator(boost::make_counting_iterator(1), Alternator()),
+               begin, end, 
+               birth, store, sd, image);
+}
+
+
+template<class BirthInfo, class SimplexData, class Field>
+template<class BI, class CI>
+typename CohomologyPersistence<BirthInfo, SimplexData, Field>::IndexDeathPair
+CohomologyPersistence<BirthInfo, SimplexData, Field>::
+add(CI coefficient_iter, BI begin, BI end, BirthInfo birth, bool store, const SimplexData& sd, bool image)
+{
     // Create simplex representation
     simplices_.push_back(SHead(sd, simplices_.empty() ? 0 : (simplices_.back().order + 1)));
     SimplexIndex    si = boost::prior(simplices_.end());
@@ -31,14 +61,15 @@
     typedef         std::list<CocycleCoefficientPair>           Candidates;
     Candidates      candidates, candidates_bulk;
     rLog(rlCohomology, "Boundary");
-    bool sign = true;       // TODO: this is very crude, fix later
+
     for (BI cur = begin; cur != end; ++cur)
     {
-        rLog(rlCohomology, "  %d %d", (*cur)->order, sign);
-        for (typename ZRow::const_iterator zcur = (*cur)->row.begin(); zcur != (*cur)->row.end(); ++zcur)
-            candidates_bulk.push_back(std::make_pair(zcur->ci, 
-                                                     (sign ? zcur->coefficient : field_.neg(zcur->coefficient))));
-        sign = !sign;
+        FieldElement coefficient = field_.init(*coefficient_iter++);
+        SimplexIndex cursi = *cur;
+
+        rLog(rlCohomology, "  %d %d", cursi->order, sign);
+        BOOST_FOREACH(const SNode& zcur, std::make_pair(cursi->row.begin(), cursi->row.end()))
+            candidates_bulk.push_back(std::make_pair(zcur.ci, field_.mul(coefficient, zcur.coefficient)));
     }
 
     candidates_bulk.sort(make_first_comparison(make_indirect_comparison(std::less<Cocycle>())));
@@ -78,7 +109,7 @@
         if (!store)
         {
             simplices_.pop_back();
-            return std::make_pair(simplices_.begin(), Death());         // TODO: shouldn't return front
+            return std::make_pair(simplices_.end(), Death());
         }
         
         signed order = 0;
--- a/include/topology/field-arithmetic.h	Wed Oct 28 14:38:36 2009 -0700
+++ b/include/topology/field-arithmetic.h	Fri Nov 06 14:19:08 2009 -0800
@@ -13,6 +13,7 @@
 
         Element     id()  const                                     { return 1; }
         Element     zero()  const                                   { return 0; }
+        Element     init(int a) const                               { return (a % p_ + p_) % p_; }
 
         Element     neg(Element a) const                            { return p_ - a; }
         Element     add(Element a, Element b) const                 { return (a+b) % p_; }