--- a/bindings/python/alphashapes2d.cpp Wed Dec 16 15:39:38 2009 -0800
+++ b/bindings/python/alphashapes2d.cpp Sun Dec 20 07:45:29 2009 -0800
@@ -9,7 +9,7 @@
namespace dp = dionysus::python;
-void fill_alpha2D_complex(bp::object points, bp::list complex)
+void fill_alpha2D_complex(bp::object points, bp::object complex)
{
typedef std::map<AlphaSimplex2D::Vertex, unsigned> ASPointMap;
@@ -32,10 +32,9 @@
for (AlphaSimplex2D::VertexContainer::const_iterator vcur = cur->vertices().begin();
vcur != cur->vertices().end(); ++vcur)
s.add(point_map[*vcur]);
-
- complex.append(s);
- complex[-1].attr("data") = cur->value();
- complex[-1].attr("attached") = cur->attached();
+
+ s.data() = bp::object(std::make_pair(cur->value(), !cur->attached())); // regular/critical rather than attached
+ complex.attr("append")(s);
}
}
--- a/bindings/python/alphashapes3d.cpp Wed Dec 16 15:39:38 2009 -0800
+++ b/bindings/python/alphashapes3d.cpp Sun Dec 20 07:45:29 2009 -0800
@@ -1,15 +1,17 @@
// Wrap includes into namespaces to avoid nameclashes
#include "../../examples/alphashapes/alphashapes3d.h"
+#include <boost/shared_ptr.hpp>
#include <boost/python.hpp>
#include <boost/python/stl_iterator.hpp>
namespace bp = boost::python;
+#include "utils.h"
#include "simplex.h" // defines SimplexVD, Vertex, and Data
namespace dp = dionysus::python;
-void fill_alpha3D_complex(bp::object points, bp::list complex)
+void fill_alpha3D_complex(bp::object points, bp::object complex)
{
typedef std::map<AlphaSimplex3D::Vertex, unsigned> ASPointMap;
@@ -29,14 +31,14 @@
for (AlphaSimplex3D::SimplexSet::const_iterator cur = simplices.begin(); cur != simplices.end(); ++cur)
{
+
dp::SimplexVD s;
for (AlphaSimplex3D::VertexContainer::const_iterator vcur = cur->vertices().begin();
vcur != cur->vertices().end(); ++vcur)
s.add(point_map[*vcur]);
-
- complex.append(s);
- complex[-1].attr("data") = cur->value();
- complex[-1].attr("attached") = cur->attached();
+
+ s.data() = bp::object(std::make_pair(cur->value(), !cur->attached())); // regular/critical rather than attached
+ complex.attr("append")(s);
}
}
--- a/bindings/python/chain.cpp Wed Dec 16 15:39:38 2009 -0800
+++ b/bindings/python/chain.cpp Sun Dec 20 07:45:29 2009 -0800
@@ -1,15 +1,21 @@
+#include <boost/iterator/indirect_iterator.hpp>
+
#include <boost/python.hpp>
#include <boost/python/iterator.hpp>
+#include <boost/python/return_internal_reference.hpp>
namespace bp = boost::python;
#include "chain.h"
namespace dp = dionysus::python;
+boost::indirect_iterator<dp::VChain::const_iterator> chain_begin(dp::VChain& c) { return boost::make_indirect_iterator(c.begin()); }
+boost::indirect_iterator<dp::VChain::const_iterator> chain_end(dp::VChain& c) { return boost::make_indirect_iterator(c.end()); }
+
void export_chain()
{
bp::class_<dp::VChain>("Chain")
- .def("__iter__", bp::iterator<dp::VChain>())
+ .def("__iter__", bp::range<bp::return_internal_reference<1> >(&chain_begin, &chain_end))
.def("__len__", &dp::VChain::size)
;
}
--- a/bindings/python/dionysus.cpp Wed Dec 16 15:39:38 2009 -0800
+++ b/bindings/python/dionysus.cpp Sun Dec 20 07:45:29 2009 -0800
@@ -1,7 +1,9 @@
#include <utilities/log.h>
#include <boost/python.hpp>
+#include "utils.h"
namespace bp = boost::python;
+namespace dp = dionysus::python;
void export_simplex();
void export_filtration();
@@ -28,6 +30,8 @@
BOOST_PYTHON_MODULE(_dionysus)
{
+ bp::to_python_converter<std::pair<double, bool>, dp::PairToTupleConverter<double, bool> >();
+
export_simplex();
export_filtration();
export_static_persistence();
--- a/bindings/python/dionysus/__init__.py Wed Dec 16 15:39:38 2009 -0800
+++ b/bindings/python/dionysus/__init__.py Sun Dec 20 07:45:29 2009 -0800
@@ -3,21 +3,19 @@
from zigzag import *
-def init_with_data(self, v, d = None):
- self._cpp_init_(v)
- if d is not None:
- self.data = d
+def init_with_none(self, iter, data = None): # convenience: data defaults to None
+ self._cpp_init_(iter, data)
def repr_with_data(self):
str = self._cpp_repr_()
- if hasattr(self, 'data'):
+ if type(self.data) == float:
str += ' %f' % self.data
return str
Simplex._cpp_init_ = Simplex.__init__
-Simplex.__init__ = init_with_data
+Simplex.__init__ = init_with_none
Simplex._cpp_repr_ = Simplex.__repr__
-Simplex.__repr__ = repr_with_data
+Simplex.__repr__ = repr_with_data
def data_cmp(s1, s2):
--- a/bindings/python/filtration.cpp Wed Dec 16 15:39:38 2009 -0800
+++ b/bindings/python/filtration.cpp Sun Dec 20 07:45:29 2009 -0800
@@ -1,44 +1,51 @@
#include <topology/filtration.h>
-#include <boost/iterator.hpp>
-#include "simplex.h"
-#include <iostream>
#include <boost/python.hpp>
+#include <boost/iterator.hpp>
+#include <boost/python/return_internal_reference.hpp>
namespace bp = boost::python;
-#include "filtration.h" // defines ListFiltration, ListTraits
+#include "simplex.h"
+#include "filtration.h" // defines PythonFiltration
#include "utils.h" // defines PythonCmp
namespace dp = dionysus::python;
-boost::shared_ptr<dp::ListFiltration> init_from_list(bp::list lst, bp::object cmp)
+boost::shared_ptr<dp::PythonFiltration> init_from_iterator(bp::object iter, bp::object cmp)
{
- boost::shared_ptr<dp::ListFiltration> p(new dp::ListFiltration(dp::ListTraits::begin(lst),
- dp::ListTraits::end(lst),
- dp::PythonCmp(cmp)));
+ typedef dp::PythonFiltration::Simplex Smplx;
+ boost::shared_ptr<dp::PythonFiltration> p(new dp::PythonFiltration(bp::stl_input_iterator<Smplx>(iter),
+ bp::stl_input_iterator<Smplx>(),
+ dp::PythonCmp(cmp)));
return p;
}
-dp::FiltrationPythonIterator
-lf_begin(const dp::ListFiltration& lf)
-{ return lf.begin(); }
+void filtration_sort(dp::PythonFiltration& f, bp::object cmp)
+{ f.sort(dp::PythonCmp(cmp)); }
-dp::FiltrationPythonIterator
-lf_end(const dp::ListFiltration& lf)
-{ return lf.end(); }
+const dp::PythonFiltration::Simplex& f_getitem(const dp::PythonFiltration& f, int i)
+{
+ if (i >= 0)
+ return f.simplex(f.begin() + i);
+ else
+ return f.simplex(f.end() + i);
+}
-unsigned
-lf_getitem(const dp::ListFiltration& lf, unsigned i)
-{ return *(lf_begin(lf) + i); }
+unsigned f_call(const dp::PythonFiltration& f, const dp::PythonFiltration::Simplex& s)
+{ return f.find(s) - f.begin(); }
void export_filtration()
{
- bp::class_<dp::ListFiltration>("Filtration", bp::no_init)
- .def("__init__", bp::make_constructor(&init_from_list))
+ bp::class_<dp::PythonFiltration>("Filtration")
+ .def("__init__", bp::make_constructor(&init_from_iterator))
- .def("__getitem__", &lf_getitem)
- .def("__iter__", bp::range(&lf_begin, &lf_end))
- .def("__len__", &dp::ListFiltration::size)
+ .def("append", &dp::PythonFiltration::push_back)
+ .def("sort", &filtration_sort)
+
+ .def("__getitem__", &f_getitem, bp::return_internal_reference<1>())
+ .def("__call__", &f_call)
+ .def("__iter__", bp::range<bp::return_internal_reference<1> >(&dp::PythonFiltration::begin, &dp::PythonFiltration::end))
+ .def("__len__", &dp::PythonFiltration::size)
;
}
--- a/bindings/python/filtration.h Wed Dec 16 15:39:38 2009 -0800
+++ b/bindings/python/filtration.h Sun Dec 20 07:45:29 2009 -0800
@@ -11,52 +11,7 @@
namespace dionysus {
namespace python {
-// ComplexTraits describing complexes of type list
-struct ListTraits
-{
- typedef bp::list Complex;
- typedef SimplexObject Simplex;
- typedef ListRandomAccessIterator<Simplex> Index;
- typedef std::less<Index> IndexComparison;
-
- typedef BinarySearchMap<SimplexVD, Index,
- SimplexVD::VertexComparison> SimplexIndexMap;
-
- static SimplexIndexMap simplex_index_map(const Complex& l) { return SimplexIndexMap(begin(l), end(l)); }
- static SimplexIndexMap simplex_index_map(Index bg, Index end) { return SimplexIndexMap(bg, end); }
-
- static unsigned size(const Complex& l) { return bp::len(l); }
- static Index begin(const Complex& l) { return Index(l, 0); }
- static Index end(const Complex& l) { return Index(l, size(l)); }
-};
-
-typedef Filtration<bp::list, unsigned, ListTraits> ListFiltration;
-
-
-// Filtration python iterator interface
-class FiltrationPythonIterator:
- public boost::iterator_adaptor<FiltrationPythonIterator, // Derived
- ListFiltration::Index, // Base
- unsigned> // Value
-{
- public:
- typedef FiltrationPythonIterator Self;
- typedef boost::iterator_adaptor<Self,
- ListFiltration::Index,
- unsigned> Parent;
-
- FiltrationPythonIterator(ListFiltration::Index i):
- Parent(i) {}
-
- private:
- friend class boost::iterator_core_access;
-
- Parent::reference dereference() const
- {
- // FIXME: I hate the const_cast here, how do I get rid of it?
- return const_cast<unsigned&>(this->base()->base().base());
- }
-};
+typedef Filtration<SimplexVD> PythonFiltration;
} } // namespace dionysus::python
--- a/bindings/python/simplex.cpp Wed Dec 16 15:39:38 2009 -0800
+++ b/bindings/python/simplex.cpp Sun Dec 20 07:45:29 2009 -0800
@@ -21,14 +21,15 @@
typename Simplex<V,T>::VertexContainer::const_iterator
vertices_end(const Simplex<V,T>& s) { return s.vertices().end(); }
-// Constructor from iterator
+// Constructor from iterator TODO: the default argument is not working yet
template<class V, class T>
-boost::shared_ptr<Simplex<V,T> > init_from_iterator(bp::object iter)
+boost::shared_ptr<Simplex<V,T> > init_from_iterator(bp::object iter, bp::object d)
{
- boost::shared_ptr<Simplex<V,T> > p(new Simplex<V,T>(bp::stl_input_iterator<V>(iter), bp::stl_input_iterator<V>()));
+ boost::shared_ptr<Simplex<V,T> > p(new Simplex<V,T>(bp::stl_input_iterator<V>(iter), bp::stl_input_iterator<V>(), d));
return p;
}
+
// Simplex hash
template<class V, class T>
size_t hash_simplex(const Simplex<V,T>& s)
@@ -42,6 +43,23 @@
return vertex_comparison(a,b) == 0;
}
+template<class S>
+bool contains(const S& s, const S& other)
+{
+ return s.contains(other);
+}
+
+template<class S>
+dp::Data get_data(const S& s)
+{
+ return s.data();
+}
+
+template<class S>
+void set_data(S& s, dp::Data d)
+{
+ s.data() = d;
+}
/* Comparisons */
// VertexComparison
@@ -59,9 +77,10 @@
.def("add", &dp::SimplexVD::add)
.add_property("boundary", bp::range(&dp::SimplexVD::boundary_begin, &dp::SimplexVD::boundary_end))
- .def("contains", &dp::SimplexVD::contains)
+ .def("contains", &contains<dp::SimplexVD>)
.def("join", (void (dp::SimplexVD::*)(const dp::SimplexVD&)) &dp::SimplexVD::join)
.def("dimension", &dp::SimplexVD::dimension)
+ .add_property("data", &get_data<dp::SimplexVD>, &set_data<dp::SimplexVD>)
.add_property("vertices", bp::range(&vertices_begin<dp::Vertex, dp::Data>, &vertices_end<dp::Vertex, dp::Data>))
.def(repr(bp::self))
--- a/bindings/python/simplex.h Wed Dec 16 15:39:38 2009 -0800
+++ b/bindings/python/simplex.h Sun Dec 20 07:45:29 2009 -0800
@@ -17,8 +17,7 @@
* SimplexObject is the representation of Python simplices in C++; i.e. it wraps bp::object and exposes a simplex-like interface.
*/
typedef int Vertex;
-// typedef double Data;
-typedef Empty<> Data;
+typedef bp::object Data;
typedef Simplex<Vertex, Data> SimplexVD;
@@ -28,13 +27,13 @@
public:
typedef SimplexObject Self;
typedef bp::object Parent;
- typedef SimplexVD::BoundaryIterator BoundaryIterator;
+ typedef bp::stl_input_iterator<Self> BoundaryIterator;
SimplexObject(Parent o = Parent()): Parent(o) {}
- BoundaryIterator boundary_begin() const { return bp::extract<const SimplexVD&>(*this)().boundary_begin(); }
- BoundaryIterator boundary_end() const { return bp::extract<const SimplexVD&>(*this)().boundary_end(); }
+ BoundaryIterator boundary_begin() const { return bp::stl_input_iterator<Self>(this->attr("boundary")); }
+ BoundaryIterator boundary_end() const { return bp::stl_input_iterator<Self>(); }
operator SimplexVD() const { return bp::extract<const SimplexVD&>(*this); }
operator bp::object() const { return *this; }
--- a/bindings/python/static-persistence.cpp Wed Dec 16 15:39:38 2009 -0800
+++ b/bindings/python/static-persistence.cpp Sun Dec 20 07:45:29 2009 -0800
@@ -14,32 +14,28 @@
/* SPersistence */
boost::shared_ptr<dp::SPersistence> init_from_filtration(bp::object f)
{
- dp::ListFiltration& lf = bp::extract<dp::ListFiltration&>(f);
- boost::shared_ptr<dp::SPersistence> p(new dp::SPersistence(lf));
+ dp::PythonFiltration& sf = bp::extract<dp::PythonFiltration&>(f);
+ boost::shared_ptr<dp::SPersistence> p(new dp::SPersistence(sf));
return p;
}
-boost::counting_iterator<dp::SPersistenceIndex> py_sp_begin(dp::SPersistence& sp) { return sp.begin(); }
-boost::counting_iterator<dp::SPersistenceIndex> py_sp_end(dp::SPersistence& sp) { return sp.end(); }
-unsigned distance(dp::SPersistence& sp,
- const dp::SPersistenceIndex& i) { return i - sp.begin(); }
-
+unsigned distance(dp::SPersistence& sp,
+ const dp::SPersistenceIndex& i) { return sp.iterator_to(i) - sp.begin(); }
-/* SPersistenceIndex */
-dp::SPersistenceIndex pair(const dp::SPersistenceIndex& i) { return i->pair; }
-bool sign(const dp::SPersistenceIndex& i) { return i->sign(); }
-const dp::VChain& cycle(const dp::SPersistenceIndex& i) { return i->cycle; }
-bool index_eq(const dp::SPersistenceIndex& i,
- const dp::SPersistenceIndex& j) { return i == j; }
+/* SPNode */
+const dp::SPersistenceNode& pair(const dp::SPersistenceNode& sn) { return *sn.pair; }
+/* PersistenceSimplexMap */
+const dp::SimplexVD& psmap_getitem(const dp::PersistenceSimplexMap& psmap,
+ const dp::SPersistenceIndex& i) { return psmap[i]; }
void export_static_persistence()
{
- bp::class_<dp::SPersistenceIndex>("SPNode")
- .add_property("pair", &pair)
- .def("sign", &sign)
- .def("cycle", &cycle, bp::return_internal_reference<1>())
- .def("__eq__", index_eq)
+ bp::class_<dp::SPersistenceNode>("SPNode", bp::no_init)
+ .def("pair", &pair, bp::return_internal_reference<1>())
+ .add_property("cycle", &dp::SPersistenceNode::cycle)
+ .def("sign", &dp::SPersistenceNode::sign)
+ .def("unpaired", &dp::SPersistenceNode::unpaired)
;
bp::class_<dp::SPersistence>("StaticPersistence", bp::no_init)
@@ -47,8 +43,13 @@
.def("pair_simplices", (void (dp::SPersistence::*)()) &dp::SPersistence::pair_simplices)
.def("__call__", &distance)
+ .def("make_simplex_map",&dp::SPersistence::make_simplex_map<dp::PythonFiltration>)
- .def("__iter__", bp::range(&py_sp_begin, &py_sp_end))
+ .def("__iter__", bp::range<bp::return_internal_reference<1> >(&dp::SPersistence::begin, &dp::SPersistence::end))
.def("__len__", &dp::SPersistence::size)
;
+
+ bp::class_<dp::PersistenceSimplexMap>("PersistenceSimplexMap", bp::no_init)
+ .def("__getitem__", &psmap_getitem, bp::return_internal_reference<1>())
+ ;
}
--- a/bindings/python/static-persistence.h Wed Dec 16 15:39:38 2009 -0800
+++ b/bindings/python/static-persistence.h Sun Dec 20 07:45:29 2009 -0800
@@ -3,12 +3,16 @@
#include <topology/static-persistence.h>
+#include "filtration.h"
+
namespace dionysus {
namespace python {
typedef StaticPersistence<> SPersistence;
typedef SPersistence::OrderElement SPersistenceNode;
typedef SPersistence::OrderIndex SPersistenceIndex;
+typedef SPersistence::SimplexMap<PythonFiltration>
+ PersistenceSimplexMap;
} } // namespace dionysus::python
--- a/bindings/python/utils.h Wed Dec 16 15:39:38 2009 -0800
+++ b/bindings/python/utils.h Sun Dec 20 07:45:29 2009 -0800
@@ -2,6 +2,7 @@
#define __PYTHON_UTILS_H__
#include <boost/python.hpp>
+#include <boost/iterator/counting_iterator.hpp>
namespace bp = boost::python;
namespace dionysus {
@@ -50,6 +51,12 @@
bp::object cmp_;
};
+template<class T1, class T2>
+struct PairToTupleConverter
+{
+ static PyObject* convert(const std::pair<T1, T2>& pair) { return bp::incref(bp::make_tuple(pair.first, pair.second).ptr()); }
+};
+
} } // namespace dionysus::python
#endif
--- a/doc/examples/alphashape.rst Wed Dec 16 15:39:38 2009 -0800
+++ b/doc/examples/alphashape.rst Sun Dec 20 07:45:29 2009 -0800
@@ -13,17 +13,18 @@
:language: python
After the points are read into the list ``points``, the functions
-:ref:`fill_alpha*_complex <alphashapes>` fill the list simplices with the
+:ref:`fill_alpha*_complex <alphashapes>` fill the :class:`Filtration` with the
simplices of the Delaunay triangulation. Each one has its :attr:`~Simplex.data`
-attribute set to its alpha shape value (the minimum value of the squared
-distance function on its dual Voronoi cell).
+attribute set to the tuple consisting of its alpha shape value (the minimum value of the squared
+distance function on its dual Voronoi cell) and whether the simplex is regular
+or critical.
-The simplices are put into lexicographic order (required for
-:class:`Filtration`), and then a filtration is created that sorts simplices with
+The filtration then sorts the simplices with
respect to their data and dimension (via :func:`data_dim_cmp`)::
- simplices.sort(vertex_cmp)
- f = Filtration(simplices, data_dim_cmp)
+ f = Filtration()
+ fill_alpha*_complex(points, f)
+ f.sort(data_dim_cmp)
We initialize :class:`StaticPersistence`, and pair the simplices::
@@ -31,6 +32,5 @@
p.pair_simplices()
Iterating over the :class:`StaticPersistence`, we output the points of the
-persistence diagram (dimension, birth, death) in the last for loop. The ``i ==
-i.pair`` condition indicates that the positive simplex is unpaired (i.e. the
-class it creates survives till infinity).
+persistence diagram (dimension, birth, death) in the last for loop. If the
+simplex is unpaired (``i.unpaired()``), the class it creates survives till infinity.
--- a/doc/examples/rips.rst Wed Dec 16 15:39:38 2009 -0800
+++ b/doc/examples/rips.rst Sun Dec 20 07:45:29 2009 -0800
@@ -25,7 +25,7 @@
distances = PairwiseDistances(points)
rips = Rips(distances)
- simplices = []
+ simplices = Filtration()
rips.generate(skeleton, max, simplices.append)
The computation of persistence and output of the persistence diagram are the
--- a/doc/examples/triangle.rst Wed Dec 16 15:39:38 2009 -0800
+++ b/doc/examples/triangle.rst Sun Dec 20 07:45:29 2009 -0800
@@ -14,13 +14,8 @@
simplices. Each :class:`Simplex` constructor takes an iterable sequence of
vertices, and optionally a data value.
-The complex must be sorted lexicographically to compute persistence using
-:class:`StaticPersistence`, and it is accomplished via a helper comparison function :func:`vertex_cmp`, which compares simplices with respect to the lexicographic ordering::
-
- complex.sort(vertex_cmp)
-
A filtration ``f`` is initialized using the :class:`Filtration` class, which
-takes a list of lexicographically sorted simplices and a comparison that defines
+takes a list of simplices (or anything iterable) and a comparison that defines
in what order the simplices should come in the filtration. In this case we use
:func:`data_cmp`, which simply compares simplices' :attr:`~Simplex.data`
attributes.
@@ -33,11 +28,10 @@
p.pair_simplices()
Subsequently, we iterate over ``p`` to access a representation of each simplex
-in the filtration order. We output each simplex, its sign, and its pair. Note
-``complex[f[p(i)]]``: ``p(i)`` gives the integer index of the iterator ``i`` in
-the filtration ``f``; in turn ``f[p(i)]`` gives the index of the simplex in the
-(lexicographically ordered) ``complex``. So the entire expression returns
-individual simplices. Naturally, one could use this to access the
+in the filtration order. We output each simplex, its sign, and its pair. The auxilliary
+``smap = p.make_simplex_map(f)`` remaps the indices of :class:`StaticPersistence` into
+the simplices in the filtration.
+Naturally, one could use this to access the
:attr:`~Simplex.data` attribute of the simplices to output the actual
persistence diagram, as is done in the :ref:`alpha-shape-example` and the
:ref:`rips-example`.
--- a/doc/python/alphashapes.rst Wed Dec 16 15:39:38 2009 -0800
+++ b/doc/python/alphashapes.rst Sun Dec 20 07:45:29 2009 -0800
@@ -17,12 +17,12 @@
.. function:: fill_alpha2D_complex(points, complex)
- Appends to the list `complex` the simplices of the 2D Delaunay triangulation
+ Appends to the `complex` the simplices of the 2D Delaunay triangulation
on the `points`.
.. function:: fill_alpha3D_complex(points, complex)
- Appends to the list `complex` the simplices of the 3D Delaunay triangulation
+ Appends to the `complex` the simplices of the 3D Delaunay triangulation
on the `points`.
@@ -34,9 +34,9 @@
from math import sin, cos, pi
points = [[cos(2*pi*t/10), sin(2*pi*t/10)] for t in xrange(10)]
- complex = []
+ complex = Filtration()
fill_alpha2D_complex(points, complex)
One can extract any given alpha shape with the usual Python list notation::
- alphashape = [s for s in complex if s.data <= .5]
+ alphashape = [s for s in complex if s.data[0] <= .5]
--- a/doc/python/filtration.rst Wed Dec 16 15:39:38 2009 -0800
+++ b/doc/python/filtration.rst Sun Dec 20 07:45:29 2009 -0800
@@ -3,30 +3,40 @@
.. class:: Filtration
- This class serves as a bridge between a complex represented as a
- lexicographically sorted list of simplices, and the
- :class:`StaticPersistence` class which needs to know the order in which the
- simplices appear in the filtration.
+ This class serves as a representation of the simplicial complex. It knows both
+ how to perform a fast lookup of a given simplex, as well as how to
+ iterate over the simplices in a sorted order.
+ .. method:: __init__()
.. method:: __init__(simplices, cmp)
- Initializes :class:`Filtration` by internally arranging the indices of
- the elements in the list `simplices` in the order sorted with respect to
- `cmp`.
+ Initializes :class:`Filtration` by internally storing the elements of the sequence
+ `simplices`, and in the order sorted with respect to `cmp`.
- .. method:: __getitem__()
+ .. method:: append(s)
+
+ Appends the given simplex `s` to the filtration.
+
+ .. method:: sort(cmp)
+
+ Sorts the filtration with respect to the comparison `cmp`.
+
+ .. method:: __getitem__(i)
Random access to the elements of the filtration.
+ .. method:: __call__(s)
+
+ Finds the integer index of the given simplex in the sorted order of the filtration.
+
.. method:: __iter__()
- Iterator over the elements of the filtration, which are simply the
- indices of the simplices in the original list `lst` sorted with respect
+ Iterator over the elements of the filtration sorted with respect
to the comparison `cmp`. E.g.::
simplices = [Simplex([0], 2), ..., Simplex([3,4,5], 3.5)]
f = Filtration(simplices, data_dim_cmp)
- for i in f: print simplices[i]
+ for s in f: print s
.. method:: __len__()
--- a/doc/python/rips.rst Wed Dec 16 15:39:38 2009 -0800
+++ b/doc/python/rips.rst Sun Dec 20 07:45:29 2009 -0800
@@ -101,11 +101,11 @@
points = [for p in points_file('...')]
distances = PairwiseDistances(points)
rips = Rips(distances)
- simplices = []
+ simplices = Filtration()
rips.generate(2, 50, simplices.append)
- f = Filtration(simplices, rips.cmp)
- p = StaticPersistence(f)
+ simplices.sort(rips.cmp)
+ p = StaticPersistence(simplices)
p.pair_simplices()
Essentially the same example is implemented in
--- a/doc/python/static-persistence.rst Wed Dec 16 15:39:38 2009 -0800
+++ b/doc/python/static-persistence.rst Sun Dec 20 07:45:29 2009 -0800
@@ -18,11 +18,17 @@
Given an SPNode in the internal representation, the method returns its
integer offset from the beginning of the filtration. This is useful to
- lookup the actual name of the simplex in the complex. For example, the
- following snippet prints out all the unpaired simplices::
+ lookup the actual name of the simplex in the complex.
+
+ .. method:: make_simplex_map(filtration)
+ Creates an auxilliary :class:`PersistenceSimplexMap` used to lookup the actual
+ simplices from the persistence indices. For example, the following
+ snippet prints out all the unpaired simplices::
+
+ smap = persistence.make_simplex_map(filtration)
for i in persistence:
- if i == i.pair: print complex[filtration[persistence(i)]]
+ if i.unpaired(): print smap[i]
.. method:: __iter__()
@@ -45,7 +51,7 @@
Returns the sign of the simplex: `True` for positive, `False` for
negative.
- .. attribute:: pair
+ .. method:: pair()
Simplex's pair. The pair is set to self if the siplex is unpaired.
@@ -56,12 +62,17 @@
container of :class:`SPNode`. For example, one can print the basis for
the (bounding) cycles::
+ smap = persistence.make_simplex_map(filtration)
for i in persistence:
- for ii in i.cycle(): print complex[filtration[persistence(ii)]]
+ for ii in i.cycle(): print smap[ii]
- .. method:: __eq__(other)
+ .. method:: unpaired()
- Returns true if the two nodes are the same. Useful for determining if
- the node is unpaired (iff ``i == i.pair``), e.g::
+ Indicates whether the simplex is unpaired.
+
+.. class:: PersistenceSimplexMap
- print len([i in persistence if i == i.pair]) # prints the number of unpaired simplices
+ .. method:: __getitem__(i)
+
+ Given a persistence index, i.e. an :class:`SPNode`, returns the
+ :class:`Simplex` it represents.
--- a/doc/tutorial.rst Wed Dec 16 15:39:38 2009 -0800
+++ b/doc/tutorial.rst Sun Dec 20 07:45:29 2009 -0800
@@ -22,23 +22,23 @@
If the points are in :math:`\mathbb{R}^2` or :math:`\mathbb{R}^3`, one can
construct an alphashape filtration::
- simplices = []
+ simplices = Filtration()
fill_alpha2D_complex(points, simplices) # for 2D, or
fill_alpha3D_complex(points, simplices) # for 3D
-Functions :ref:`fill_alpha*_complex <alphashapes>` fill the ``simplices`` list
-with all the :class:`simplices <Simplex>` of the Delaunay triangulation. Each one has attributes
-``data`` and ``attached`` set to its respective value in the alpha shape (the
-smallest value of the squared distance function on the dual Voronoi cell) and
-whether the simplex is attached or not (i.e. whether its dual cell does not or
-does contain a critical point of the distance function). See :ref:`alphashapes`
+Functions :ref:`fill_alpha*_complex <alphashapes>` fill the ``simplices``
+with all the :class:`simplices <Simplex>` of the Delaunay triangulation.
+Each one has its attribute ``data`` set to a pair: the
+smallest value of the squared distance function on the dual Voronoi cell and
+whether the simplex is critical or not (i.e. whether its dual cell does or
+does not contain a critical point of the distance function). See :ref:`alphashapes`
for more details.
As a result, if one wanted only those simplices whose alpha shape value did not
exceed 10, one could obtain them as follows::
- simplices10 = [s for s in simplices if s.data <= 10]
+ simplices10 = [s for s in simplices if s.data[0] <= 10]
If the point set lies in higher dimensions, one may construct a Rips complex on
it. This complex requires only pairwise distances, which makes it very
@@ -78,44 +78,42 @@
^^^^^^^
For the first approach, i.e. to use :class:`StaticPersistence`, one must put the
-simplices into lexicographic ordering (sort them with respect to
-:func:`vertex_cmp`), and construct a filtration with respect to some ordering
+sort the filtration with respect to some ordering
(for example, :func:`data_dim_cmp` for alpha shapes or :meth:`Rips.cmp` for the
Rips complex)::
- simplices.sort(vertex_cmp)
- f = Filtration(simplices, data_dim_cmp) # for the alpha shapes
- f = Filtration(simplices, rips.cmp) # for the rips complex
+ simplices.sort(data_dim_cmp) # for the alpha shapes
+ simplices.sort(rips.cmp) # for the rips complex
Creating an instance of :class:`StaticPersistence` initialized with the
filtration really initializes a boundary matrix. The subsequent call to
:meth:`~StaticPersistence.pair_simplices` reduces the matrix to compute the
persistence of the filtration::
- p = StaticPersistence(f)
+ p = StaticPersistence(simplices)
p.pair_simplices()
Once the simplices are paired, one may examine the pairings by iterating over
-the instance of :class:`StaticPersistence`::
+the instance of :class:`StaticPersistence`. We can use an auxilliary map ``smap``
+to remap the persistence indices into the actual simplices::
+ smap = p.make_simplex_map(simplices)
for i in p:
if i.sign():
- birth = simplices[f[p(i)]]
- if i == i.pair:
+ birth = smap[i]
+ if i.unpaired():
print birth.dimension(), birth.data, "inf"
continue
- death = simplices[f[p(i.pair)]]
+ death = smap[i.pair()]
print birth.dimension(), birth.data, death.data
The iterator ``i`` over the instance ``p`` of :class:`StaticPersistence` is of type
:class:`SPNode`, and represents individual simplices taken in the filtration
-order. It knows about its own :meth:`~SPNode.sign` and :attr:`~SPNode.pair`
-(``i.pair == i`` if the simplex is unpaired). Calling ``p`` with the iterator
-``i`` (``p(i)``) gives its integer index in the filtration, which can then be
-used with the instance ``f`` of :class:`Filtration` to get the integer index in
-the lexicographically sorted list `simplices`: ``simplices[f[p(i)]]``. The
-previous code snippet prints out the persistence diagrams of the given
+order. It knows about its own :meth:`~SPNode.sign` and :meth:`~SPNode.pair` as well as
+whether it is :math:`~SPNode.unpaired`. :meth:`StaticPersistence.make_simplex_map` creates
+a map that we use to get the actual simplices: ``smap[i]`` and ``smap[i.pair()]``.
+The previous code snippet prints out the persistence diagrams of the given
filtration.
@@ -185,14 +183,14 @@
* To avoid the computation of simplex sizes in the Rips complex during the
initialization of a :class:`Filtration`, store them explicitly in
:attr:`Simplex.data` attribute (this is not done by default to save memory);
- then use :func:`data_dim_cmp` in the initialization of the
+ then use :func:`data_dim_cmp` when sorting the
:class:`Filtration`::
rips = Rips(distances)
- simplices = []
+ simplices = Filtration()
rips.generate(..., simplices.append)
for s in simplices: s.data = rips.eval(s)
- f = Filtration(simplices, data_dim_cmp)
+ simplices.sort(data_dim_cmp)
--- a/examples/alphashapes/alphashapes.py Wed Dec 16 15:39:38 2009 -0800
+++ b/examples/alphashapes/alphashapes.py Sun Dec 20 07:45:29 2009 -0800
@@ -12,16 +12,15 @@
exit()
points = [p for p in points_file(argv[1])]
-simplices = []
+f = Filtration()
if len(points[0]) == 2: # 2D
- fill_alpha2D_complex(points, simplices)
+ fill_alpha2D_complex(points, f)
elif len(points[1]) == 3: # 3D
- fill_alpha3D_complex(points, simplices)
+ fill_alpha3D_complex(points, f)
-simplices.sort(vertex_cmp) # Must ensure lexicographic ordering
-print "Total number of simplices:", len(simplices)
+print "Total number of simplices:", len(f)
-f = Filtration(simplices, data_dim_cmp)
+f.sort(data_dim_cmp)
print "Filtration initialized"
p = StaticPersistence(f)
@@ -31,12 +30,13 @@
print "Simplices paired"
print "Outputting persistence diagram"
+smap = p.make_simplex_map(f)
for i in p:
if i.sign():
- b = simplices[f[p(i)]]
- if i == i.pair:
- print b.dimension(), sqrt(b.data), "inf"
+ b = smap[i]
+ if i.unpaired():
+ print b.dimension(), sqrt(b.data[0]), "inf"
continue
- d = simplices[f[p(i.pair)]]
- print b.dimension(), sqrt(b.data), sqrt(d.data)
+ d = smap[i.pair()]
+ print b.dimension(), sqrt(b.data[0]), sqrt(d.data[0])
--- a/examples/cohomology/rips-pairwise-cohomology.py Wed Dec 16 15:39:38 2009 -0800
+++ b/examples/cohomology/rips-pairwise-cohomology.py Sun Dec 20 07:45:29 2009 -0800
@@ -1,6 +1,6 @@
#!/usr/bin/env python
-from dionysus import Simplex, CohomologyPersistence, points_file, PairwiseDistances, ExplicitDistances, Rips
+from dionysus import Simplex, CohomologyPersistence, points_file, PairwiseDistances, ExplicitDistances, Rips, data_dim_cmp
from sys import argv, exit
import time
--- a/examples/lsfiltration.py Wed Dec 16 15:39:38 2009 -0800
+++ b/examples/lsfiltration.py Sun Dec 20 07:45:29 2009 -0800
@@ -21,29 +21,26 @@
vertices.append(float(line.split()[1]))
# Read simplices
- simplices = []
+ fltr = Filtration()
with open(simplices_filename) as f:
for line in f:
if line.startswith('#'): continue
- simplices.append(map(int, line.split()))
-
- # Setup the complex
- complex = [Simplex(s) for s in simplices]
- complex.sort(vertex_cmp)
+ fltr.append(Simplex(map(int, line.split())))
+ fltr.sort(lambda x,y: max_vertex_cmp(x,y,vertices))
# Compute persistence
- f = Filtration(complex, lambda x,y: max_vertex_cmp(x,y,vertices))
- p = StaticPersistence(f)
+ p = StaticPersistence(fltr)
p.pair_simplices()
# Output the persistence diagram
+ smap = p.make_simplex_map(fltr)
for i in p:
if not i.sign(): continue
- b = complex[f[p(i)]]
- d = complex[f[p(i.pair)]]
+ b = smap[i]
+ d = smap[i.pair()]
- if i == i.pair:
+ if i.unpaired():
print b.dimension(), max_vertex(b, vertices), "inf"
continue
--- a/examples/rips/rips-pairwise.py Wed Dec 16 15:39:38 2009 -0800
+++ b/examples/rips/rips-pairwise.py Sun Dec 20 07:45:29 2009 -0800
@@ -13,7 +13,7 @@
rips = Rips(distances)
print time.asctime(), "Rips initialized"
- simplices = []
+ simplices = Filtration()
rips.generate(skeleton, max, simplices.append)
print time.asctime(), "Generated complex: %d simplices" % len(simplices)
@@ -22,27 +22,28 @@
for s in simplices: s.data = rips.eval(s)
print time.asctime(), simplices[0], '...', simplices[-1]
- f = Filtration(simplices, data_dim_cmp) # could be rips.cmp if s.data for s in simplices is not set
+ simplices.sort(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)
+ p = StaticPersistence(simplices)
print time.asctime(), "Initialized StaticPersistence"
p.pair_simplices()
print time.asctime(), "Simplices paired"
print "Outputting persistence diagram"
+ smap = p.make_simplex_map(simplices)
for i in p:
if i.sign():
- b = simplices[f[p(i)]]
+ b = smap[i]
if b.dimension() >= skeleton: continue
- if i == i.pair:
+ if i.unpaired():
print b.dimension(), b.data, "inf"
continue
- d = simplices[f[p(i.pair)]]
+ d = smap[i.pair()]
print b.dimension(), b.data, d.data
if __name__ == '__main__':
--- a/examples/rips/rips.py Wed Dec 16 15:39:38 2009 -0800
+++ b/examples/rips/rips.py Sun Dec 20 07:45:29 2009 -0800
@@ -11,8 +11,8 @@
dist = Distances()
r = Rips(dist)
-lst = []
-lst2 = []
+lst = Filtration()
+lst2 = Filtration()
#enable_log('rips')
@@ -33,8 +33,10 @@
r.vertex_cofaces(2, 1, 3, cofaces.append, [0,2,4])
print "Cofaces of vertex 2 on vertices [0,2,4]:", cofaces
-f = Filtration(lst, r.cmp)
+f = lst
+f.sort(r.cmp)
p = StaticPersistence(f)
p.pair_simplices()
+smap = p.make_simplex_map(f)
for s in p:
- print lst[f[p(s)]], s.sign()
+ print smap[s], s.sign()
--- a/examples/triangle/triangle.py Wed Dec 16 15:39:38 2009 -0800
+++ b/examples/triangle/triangle.py Sun Dec 20 07:45:29 2009 -0800
@@ -14,16 +14,18 @@
print "Data: ", sorted(complex, data_cmp)
print "DataDim:", sorted(complex, data_dim_cmp)
-complex.sort(vertex_cmp) # put complex in the lexicographic order; required for persistence computation
-
f = Filtration(complex, data_cmp)
-print "Complex in the filtration order:", ', '.join((str(complex[i]) for i in f))
+print "Complex in the filtration order:", ', '.join((str(s) for s in f))
p = StaticPersistence(f)
+print "Persistence initialized"
p.pair_simplices()
+print "Simplices paired"
+
+smap = p.make_simplex_map(f)
for i in p:
- print "%s (%d) - %s (%d)" % (complex[f[p(i)]], i.sign(),
- complex[f[p(i.pair)]], i.pair.sign())
- print "Cycle (%d):" % len(i.cycle()), " + ".join((str(complex[f[p(ii)]]) for ii in i.cycle()))
+ print i.sign(), i.pair().sign()
+ print "%s (%d) - %s (%d)" % (smap[i], i.sign(), smap[i.pair()], i.pair().sign())
+ print "Cycle (%d):" % len(i.cycle), " + ".join((str(smap[ii]) for ii in i.cycle))
-print "Number of unpaired simplices:", len([i for i in p if i == i.pair])
+print "Number of unpaired simplices:", len([i for i in p if i.unpaired()])