--- a/bindings/python/dionysus/__init__.py Mon Jul 25 23:21:29 2011 -0700
+++ b/bindings/python/dionysus/__init__.py Thu May 10 14:24:40 2012 -0700
@@ -34,15 +34,11 @@
def vertex_dim_cmp(s1, s2):
return cmp(s1.dimension(), s2.dimension()) or vertex_cmp(s1, s2)
-# TBD: Port this into C++
-def point_iterator( point ):
-
- yield point.x
- yield point.y
- if not point.data is None:
- yield point.data
-
-Point.__iter__ = point_iterator
+def fill_alpha_complex(points, simplices):
+ if len(points[0]) == 2: # 2D
+ fill_alpha2D_complex(points, simplices)
+ elif len(points[0]) == 3: # 3D
+ fill_alpha3D_complex(points, simplices)
def closure(simplices, k):
"""Compute the k-skeleton of the closure of the list of simplices."""
--- a/bindings/python/filtration.h Mon Jul 25 23:21:29 2011 -0700
+++ b/bindings/python/filtration.h Thu May 10 14:24:40 2012 -0700
@@ -8,7 +8,7 @@
namespace bp = boost::python;
-namespace dionysus {
+namespace dionysus {
namespace python {
typedef Filtration<SimplexVD> PythonFiltration;
--- a/bindings/python/persistence-diagram.cpp Mon Jul 25 23:21:29 2011 -0700
+++ b/bindings/python/persistence-diagram.cpp Thu May 10 14:24:40 2012 -0700
@@ -1,41 +1,88 @@
#include<topology/persistence-diagram.h>
#include<utilities/types.h>
+#include "filtration.h"
+#include "simplex.h"
+#include "static-persistence.h"
+#include "dynamic-persistence.h"
+
+#include <boost/foreach.hpp>
+
#include<boost/python.hpp>
#include<boost/python/init.hpp>
#include<boost/shared_ptr.hpp>
#include<boost/python/stl_iterator.hpp>
#include<boost/python/def.hpp>
-
-
namespace bp = boost::python;
namespace dionysus{
namespace python{
-typedef bp::object Data;
-typedef PDPoint<Data> PointD;
-typedef PersistenceDiagram<Data> PersistenceDiagramD;
+typedef bp::object Data;
+typedef PersistenceDiagram<Data> PersistenceDiagramD;
+typedef PersistenceDiagramD::Point PointD;
+typedef boost::shared_ptr<PersistenceDiagramD> PDgmPtr;
} } //namespace dionysus::python
namespace dp = dionysus::python;
-// Easy way out of overloads
-RealType get_x (const dp::PointD& p) { return p.x(); }
-RealType get_y (const dp::PointD& p) { return p.y(); }
-bp::object get_data(const dp::PointD& p) { return p.data(); }
+struct PointFromTupleConverter
+{
+ PointFromTupleConverter()
+ {
+ boost::python::converter::registry::push_back(&convertible,
+ &construct,
+ boost::python::type_id<dp::PointD>());
+ }
+
+ static void* convertible(PyObject* obj_ptr)
+ {
+ if (!PyTuple_Check(obj_ptr)) return 0;
+ if (PyTuple_Size(obj_ptr) < 2) return 0;
+ return obj_ptr;
+ }
+
+ static void construct(PyObject* obj_ptr,
+ boost::python::converter::rvalue_from_python_stage1_data* data)
+ {
+ //const char* value = PyString_AsString(obj_ptr);
+ //if (value == 0) boost::python::throw_error_already_set();
+
+ // Grab pointer to memory into which to construct the new T
+ void* storage = ( (boost::python::converter::rvalue_from_python_storage<dp::PointD>*) data)->storage.bytes;
+
+ RealType x = bp::extract<RealType>(PyTuple_GetItem(obj_ptr, 0));
+ RealType y = bp::extract<RealType>(PyTuple_GetItem(obj_ptr, 1));
+
+ // in-place construct the new T using the character data extraced from the python object
+ dp::PointD* p = new (storage) dp::PointD(x,y);
+
+ if (PyTuple_Size(obj_ptr) > 2)
+ p->data() = bp::extract<bp::object>(PyTuple_GetItem(obj_ptr, 2));
+
+ // Stash the memory chunk pointer for later use by boost.python
+ data->convertible = storage;
+ }
+};
+
+struct PointToTupleConverter
+{
+ static PyObject* convert(const dp::PointD& p)
+ {
+ if (p.data().ptr() == bp::object().ptr())
+ return bp::incref(bp::make_tuple(p.x(), p.y()).ptr());
+ else
+ return bp::incref(bp::make_tuple(p.x(), p.y(), p.data()).ptr());
+ }
+};
+
void export_point( )
{
- bp::class_<dp::PointD>("Point")
- .def( bp::init<RealType, RealType, bp::optional<dp::Data> >())
- .add_property("x", &get_x)
- .add_property("y", &get_y)
- .add_property("data", &get_data)
- .def( repr(bp::self))
- ;
+ PointFromTupleConverter();
+ bp::to_python_converter<dp::PointD, PointToTupleConverter>();
}
@@ -55,11 +102,106 @@
RealType bottleneck_distance_adapter(const dp::PersistenceDiagramD& dgm1, const dp::PersistenceDiagramD& dgm2)
{
return bottleneck_distance(dgm1, dgm2);
-}
+}
+
+
+template<class Persistence>
+struct InitDiagrams
+{
+ typedef std::map<int, dp::PDgmPtr> DiagramMap;
+ typedef typename Persistence::template SimplexMap<dp::PythonFiltration> SMap;
+
+ struct DataEvaluator
+ {
+ DataEvaluator(const SMap& smap_):
+ smap(smap_) {}
+
+ template<class Key>
+ RealType operator()(Key k) const { return bp::extract<RealType>(smap[k].data()); }
+
+ const SMap& smap;
+ };
+
+ struct PythonEvaluator
+ {
+ PythonEvaluator(const SMap& smap_, bp::object eval_):
+ smap(smap_), eval(eval_) {}
+
+ template<class Key>
+ RealType operator()(Key k) const { return bp::extract<RealType>(eval(smap[k])); }
+
+ const SMap& smap;
+ bp::object eval;
+ };
+
+ // A hack
+ struct DiagramMapOwner: public DiagramMap
+ {
+ typedef dp::PersistenceDiagramD mapped_type;
+
+ mapped_type& operator[](Dimension d)
+ {
+ if (this->find(d) == this->end())
+ this->insert(std::make_pair(d, dp::PDgmPtr(new dp::PersistenceDiagramD(d))));
+ return *DiagramMap::operator[](d);
+ }
+ };
+
+ static
+ bp::list extract_list(const DiagramMapOwner& dgms)
+ {
+ bp::list result;
+ size_t dim = 0;
+ typedef typename DiagramMapOwner::value_type ValType;
+ BOOST_FOREACH(const ValType& dim_dgm, dgms)
+ {
+ while (dim_dgm.first > dim)
+ {
+ result.append(dp::PDgmPtr(new dp::PersistenceDiagramD));
+ ++dim;
+ }
+
+ // dim_dgm.first == dim
+ result.append(dim_dgm.second);
+ dim++;
+ }
+
+ return result;
+ }
+
+ struct PointDataVisitor
+ {
+ PointDataVisitor(bp::object data_): data(data_) {}
+ void point(const typename Persistence::iterator& i, dp::PointD& p) const { p.data() = data(*i); }
+ bp::object data;
+ };
+
+ static
+ bp::list init(const Persistence& p, const dp::PythonFiltration& f, bp::object eval, bp::object data)
+ {
+
+ DiagramMapOwner dgms;
+ SMap smap = p.make_simplex_map(f);
+ if (eval == bp::object())
+ init_diagrams(dgms, p.begin(), p.end(),
+ DataEvaluator(smap),
+ evaluate_through_map(smap, dp::SimplexVD::DimensionExtractor()));
+ else if (data == bp::object())
+ init_diagrams(dgms, p.begin(), p.end(),
+ PythonEvaluator(smap, eval),
+ evaluate_through_map(smap, dp::SimplexVD::DimensionExtractor()));
+ else
+ init_diagrams(dgms, p.begin(), p.end(),
+ PythonEvaluator(smap, eval),
+ evaluate_through_map(smap, dp::SimplexVD::DimensionExtractor()),
+ PointDataVisitor(data));
+ return extract_list(dgms);
+ }
+};
void export_persistence_diagram()
{
- bp::class_<dp::PersistenceDiagramD>("PersistenceDiagram")
+ bp::class_<dp::PersistenceDiagramD, dp::PDgmPtr>("PersistenceDiagram")
.def("__init__", bp::make_constructor(&init_from_points_sequence))
.def( bp::init<Dimension>())
.def("append", &dp::PersistenceDiagramD::push_back)
@@ -69,5 +211,16 @@
.def("__len__", &dp::PersistenceDiagramD::size)
;
+ bp::def("init_diagrams", &InitDiagrams<dp::SPersistence>::init,
+ (bp::arg("persistence"),
+ bp::arg("filtration"),
+ bp::arg("eval")=bp::object(),
+ bp::arg("data")=bp::object()));
+ bp::def("init_diagrams", &InitDiagrams<dp::DPersistenceChains>::init,
+ (bp::arg("persistence"),
+ bp::arg("filtration"),
+ bp::arg("eval")=bp::object(),
+ bp::arg("data")=bp::object()));
+
bp::def("bottleneck_distance", &bottleneck_distance_adapter);
}
--- a/bindings/python/static-persistence.cpp Mon Jul 25 23:21:29 2011 -0700
+++ b/bindings/python/static-persistence.cpp Thu May 10 14:24:40 2012 -0700
@@ -25,7 +25,7 @@
bp::class_<dp::SPersistence>("StaticPersistence", bp::no_init)
.def("__init__", bp::make_constructor(&dp::init_from_filtration<dp::SPersistence>))
-
+
.def("pair_simplices", &pair_simplices)
.def("__call__", &dp::distance<dp::SPersistence, dp::SPersistenceIndex>)
.def("make_simplex_map",&dp::SPersistence::make_simplex_map<dp::PythonFiltration>)
--- a/doc/python/alphashapes.rst Mon Jul 25 23:21:29 2011 -0700
+++ b/doc/python/alphashapes.rst Thu May 10 14:24:40 2012 -0700
@@ -15,13 +15,18 @@
.. _`Delaunay triangulation`: http://en.wikipedia.org/wiki/Delaunay_triangulation
+.. function:: fill_alpha_complex(points, complex)
+
+ Based on the dimension of the first point, decides which of the two functions
+ below to call.
+
.. function:: fill_alpha2D_complex(points, complex)
-
+
Appends to the `complex` the simplices of the 2D Delaunay triangulation
on the `points`.
.. function:: fill_alpha3D_complex(points, complex)
-
+
Appends to the `complex` the simplices of the 3D Delaunay triangulation
on the `points`.
--- a/doc/python/persistence-diagram.rst Mon Jul 25 23:21:29 2011 -0700
+++ b/doc/python/persistence-diagram.rst Thu May 10 14:24:40 2012 -0700
@@ -4,20 +4,24 @@
.. class:: PersistenceDiagram
.. method:: __init__( dimension )
-
+
Initializes : an empty( no points ) :class:`PersistenceDiagram` object and sets
the :attr:`~PersistenceDiagram.dimension` attribute( must be integer ) e.g.::
-
+
dia = PersistenceDiagram( 1 )
.. method:: __init__( dimension, point_seq )
-
- Initializes :class:`PersistenceDiagram` of specified dimension from the given sequence `seq` of :class:`Point` objects, e.g.::
-
- dia = PersistenceDiagram( 1, [Point(1,2)] )
+
+ Initializes :class:`PersistenceDiagram` of specified dimension from the given sequence `seq` of tuples, e.g.::
+
+ dia = PersistenceDiagram( 1, (1,2) )
+
+ The tuples must have at least 2 elements ``(birth, death)``.
+ If there is a third element, it is stored as extra data associated to
+ a point.
.. method:: append( p )
-
+
Adds point `p` to the persistence diagram.
.. attribute:: dimension
@@ -28,7 +32,7 @@
Iterator over the points in the persistence diagram,
e.g.::
-
+
for p in dia: print p
.. method:: __len__( )
@@ -40,42 +44,19 @@
Utility functions for persistence diagrams
--------------------------------------------
+.. function:: init_diagrams(persistence, filtration[, eval = lambda s: s.data[, data = lambda i: None]])
+
+ Initializes a collection of :class:`PersistenceDiagram`s from `persistence`
+ and `filtration`. Optional `eval` can determine how to extract birth and
+ death values from a simplex. For example, if `filtration` was filled using
+ :fun:`fill_alpha_complex()`, the :attr:`~Simplex.data` contains a pair ``(value, critical)``.
+ We can extract the value from the tuple::
+
+ init_diagrams(persistence, filtration, lambda s: s.data[0])
+
+ Optional `data` argument can return arbitrary data to associate with each point,
+ given an node of `persistence`.
.. function:: bottleneck_distance(dia1, dia2)
-
- Calculates the bottleneck distance between the two persistence diagrams.
-
-
-.. class:: Point
-
- .. method:: __init__( x, y [, data])
-
- Initializes :class:`Point` with the given real-valued coordinates and
- optionally real value `data`, e.g.::
-
- s = Point( 1, 1.1, 1.11 )
-
- .. attribute:: x
-
- x coordinate of the point
-
- .. attribute:: y
-
- y coordinate of the point
-
- .. attribute:: data
-
- Real value stored in the simplex.
-
- .. method:: __iter__( )
-
- Point objects are iterable, returning two or three elements depending on presence of data, e.g.::
-
- p = Point( 1, 1.1, 1.11 )
- for i in p: print p
-
- 1
- 1.1
- 1.11
-
+ Calculates the bottleneck distance between the two persistence diagrams.
--- a/include/topology/persistence-diagram.h Mon Jul 25 23:21:29 2011 -0700
+++ b/include/topology/persistence-diagram.h Thu May 10 14:24:40 2012 -0700
@@ -32,7 +32,7 @@
Data& data() { return point_.second(); }
std::ostream& operator<<(std::ostream& out) const { return (out << x() << " " << y()); } // << " " << data()); }
-
+
struct Visitor
{
template<class Iterator>
@@ -45,17 +45,17 @@
private:
boost::compressed_pair<std::pair<RealType, RealType>, Data> point_;
-
+
private:
- /* Serialization */
- friend class boost::serialization::access;
-
- template<class Archive>
- void serialize(Archive& ar, version_type );
+ /* Serialization */
+ friend class boost::serialization::access;
+
+ template<class Archive>
+ void serialize(Archive& ar, version_type );
};
template<class Data>
-std::ostream& operator<<(std::ostream& out, const PDPoint<Data>& point)
+std::ostream& operator<<(std::ostream& out, const PDPoint<Data>& point)
{ return (point.operator<<(out)); }
template<class Point, class Iterator, class Evaluator, class Visitor>
@@ -71,7 +71,7 @@
/**
* Class: PersistenceDiagram
*
- * Stores birth-death pairs, i.e. points in the extended plane. Each point can also store
+ * Stores birth-death pairs, i.e. points in the extended plane. Each point can also store
* additional information described by `Data_` template parameter.
*/
template<class Data_ = Empty<> >
@@ -85,24 +85,24 @@
PersistenceDiagram() {}
- PersistenceDiagram( Dimension dimension ):
+ PersistenceDiagram( Dimension dimension ):
dimension_( dimension ) {}
template<class OtherData>
PersistenceDiagram(const PersistenceDiagram<OtherData>& other);
template<class Iterator, class Evaluator>
- PersistenceDiagram(Iterator bg, Iterator end,
+ PersistenceDiagram(Iterator bg, Iterator end,
const Evaluator& eval = Evaluator());
template<class Iterator, class Evaluator, class Visitor>
- PersistenceDiagram(Iterator bg, Iterator end,
- const Evaluator& eval = Evaluator(),
+ PersistenceDiagram(Iterator bg, Iterator end,
+ const Evaluator& eval = Evaluator(),
const Visitor& visitor = Visitor());
template<class Iterator, class Evaluator, class Visitor>
- void init(Iterator bg, Iterator end,
- const Evaluator& eval = Evaluator(),
+ void init(Iterator bg, Iterator end,
+ const Evaluator& eval = Evaluator(),
const Visitor& visitor = Visitor());
const_iterator begin() const { return points_.begin(); }
@@ -110,7 +110,7 @@
size_t size() const { return points_.size(); }
void push_back(const Point& point) { points_.push_back(point); }
-
+
std::ostream& operator<<(std::ostream& out) const;
Dimension dimension() const { return dimension_; }
@@ -118,53 +118,53 @@
private:
PointVector points_;
Dimension dimension_;
-
+
private:
- /* Serialization */
- friend class boost::serialization::access;
-
- template<class Archive>
- void serialize(Archive& ar, version_type );
+ /* Serialization */
+ friend class boost::serialization::access;
+
+ template<class Archive>
+ void serialize(Archive& ar, version_type );
};
template<class Data>
-std::ostream& operator<<(std::ostream& out, const PersistenceDiagram<Data>& pd)
+std::ostream& operator<<(std::ostream& out, const PersistenceDiagram<Data>& pd)
{ return (pd.operator<<(out)); }
// Function: init_diagram_vector()
template<class Diagrams, class Iterator, class Evaluator, class DimensionExtractor>
void init_diagrams(Diagrams& diagrams,
- Iterator bg, Iterator end,
- const Evaluator& evaluator = Evaluator(),
+ Iterator bg, Iterator end,
+ const Evaluator& evaluator = Evaluator(),
const DimensionExtractor& dimension = DimensionExtractor());
template<class Diagrams, class Iterator, class Evaluator, class DimensionExtractor, class Visitor>
void init_diagrams(Diagrams& diagrams,
- Iterator bg, Iterator end,
- const Evaluator& evaluator = Evaluator(),
+ Iterator bg, Iterator end,
+ const Evaluator& evaluator = Evaluator(),
const DimensionExtractor& dimension = DimensionExtractor(),
const Visitor& visitor = Visitor());
-// Class: Linfty
+// Class: Linfty
// Functor that computes L infinity norm between two points
template<class Point1, class Point2>
struct Linfty
{
- RealType operator()(const Point1& p1, const Point2& p2) const { return std::max(std::abs(p1.x() - p2.x()),
+ RealType operator()(const Point1& p1, const Point2& p2) const { return std::max(std::abs(p1.x() - p2.x()),
std::abs(p1.y() - p2.y())); }
-
+
template<class Point>
RealType diagonal(const Point& p) const { return std::abs(p.y() - p.x())/2; }
};
// Function: bottleneck_distance(dgm1, dgm2)
// Computes bottleneck distance between the two diagrams.
-template<class Diagram1,
- class Diagram2,
+template<class Diagram1,
+ class Diagram2,
class Norm>
RealType bottleneck_distance(const Diagram1& dgm1, const Diagram2& dgm2, const Norm& norm = Norm());
-template<class Diagram1,
+template<class Diagram1,
class Diagram2>
RealType bottleneck_distance(const Diagram1& dgm1, const Diagram2& dgm2)
{ return bottleneck_distance(dgm1, dgm2, Linfty<typename Diagram1::Point, typename Diagram2::Point>()); }
--- a/include/topology/persistence-diagram.hpp Mon Jul 25 23:21:29 2011 -0700
+++ b/include/topology/persistence-diagram.hpp Thu May 10 14:24:40 2012 -0700
@@ -62,10 +62,10 @@
RealType y = Infinity;
if (&*(i->pair) != &*i)
y = evaluator(&*(i->pair));
-
+
Point p(x,y);
visitor.point(i, p);
-
+
if (x == y) return boost::optional<Point>();
return p;
@@ -116,8 +116,8 @@
template<class Archive>
void
PDPoint<D>::
-serialize(Archive& ar, version_type )
-{
+serialize(Archive& ar, version_type )
+{
ar & make_nvp("x", x());
ar & make_nvp("y", y());
ar & make_nvp("data", data());
@@ -127,9 +127,9 @@
template<class Archive>
void
PersistenceDiagram<D>::
-serialize(Archive& ar, version_type )
-{
- ar & make_nvp("points", points_);
+serialize(Archive& ar, version_type )
+{
+ ar & make_nvp("points", points_);
}
@@ -148,7 +148,7 @@
Edge(unsigned v1, unsigned v2, RealType d):
Parent(v1, v2), distance(d) {}
- bool operator<(const Edge& other) const { return distance < other.distance; }
+ bool operator<(const Edge& other) const { return distance < other.distance; }
RealType distance;
};
@@ -242,6 +242,6 @@
EV_counting_const_iterator(edges.end()),
edges.begin(),
CardinaliyComparison(max_size, edges.begin()));
-
+
return (*bdistance)->distance;
}