--- 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)