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