--- a/bindings/python/alphashapes2d.cpp Thu Dec 03 13:21:23 2009 -0800
+++ b/bindings/python/alphashapes2d.cpp Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/bindings/python/alphashapes3d.cpp Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/bindings/python/chain.cpp Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/bindings/python/dionysus.cpp Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/bindings/python/dionysus/__init__.py Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/bindings/python/filtration.cpp Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/bindings/python/filtration.h Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/bindings/python/simplex.cpp Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/bindings/python/simplex.h Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/bindings/python/static-persistence.cpp Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/bindings/python/static-persistence.h Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/bindings/python/utils.h Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/doc/examples/alphashape.rst Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/doc/examples/rips.rst Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/doc/examples/triangle.rst Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/doc/python/alphashapes.rst Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/doc/python/filtration.rst Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/doc/python/rips.rst Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/doc/python/static-persistence.rst Sun Dec 20 08:27:01 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 Thu Dec 03 13:21:23 2009 -0800
+++ b/doc/tutorial.rst Sun Dec 20 08:27:01 2009 -0800
@@ -22,23 +22,24 @@
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`
-for more details, and :ref:`alpha-shape-example` for a full example.
+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, and :ref:`alpha-shape-example` for a
+full example.
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 +79,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.
@@ -190,14 +189,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/CMakeLists.txt Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/CMakeLists.txt Sun Dec 20 08:27:01 2009 -0800
@@ -4,7 +4,7 @@
add_subdirectory (consistency)
add_subdirectory (cohomology)
add_subdirectory (fitness)
-#add_subdirectory (grid)
+add_subdirectory (grid)
add_subdirectory (triangle)
add_subdirectory (poincare)
add_subdirectory (rips)
--- a/examples/alphashapes/alphashapes.py Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/alphashapes/alphashapes.py Sun Dec 20 08:27:01 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/alphashapes/alphashapes2d.cpp Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/alphashapes/alphashapes2d.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -9,7 +9,7 @@
#include <fstream>
-typedef Filtration<AlphaSimplex2DVector> AlphaFiltration;
+typedef Filtration<AlphaSimplex2D> AlphaFiltration;
typedef StaticPersistence<> Persistence;
typedef PersistenceDiagram<> PDgm;
@@ -41,12 +41,12 @@
}
rInfo("Delaunay triangulation computed");
- AlphaSimplex2DVector complex;
- fill_complex(Dt, complex);
- rInfo("Simplices: %i", complex.size());
+ AlphaFiltration af;
+ fill_complex(Dt, af);
+ rInfo("Simplices: %i", af.size());
// Create the alpha-shape filtration
- AlphaFiltration af(complex.begin(), complex.end(), AlphaSimplex2D::AlphaOrder());
+ af.sort(AlphaSimplex2D::AlphaOrder());
rInfo("Filtration initialized");
Persistence p(af);
@@ -55,12 +55,11 @@
p.pair_simplices();
rInfo("Simplices paired");
- std::map<Dimension, PDgm> dgms;
+ Persistence::SimplexMap<AlphaFiltration> m = p.make_simplex_map(af);
+ std::map<Dimension, PDgm> dgms;
init_diagrams(dgms, p.begin(), p.end(),
- evaluate_through_map(OffsetMap<Persistence::OrderIndex, AlphaFiltration::Index>(p.begin(), af.begin()),
- evaluate_through_filtration(af, AlphaSimplex2D::AlphaValueEvaluator())),
- evaluate_through_map(OffsetMap<Persistence::OrderIndex, AlphaFiltration::Index>(p.begin(), af.begin()),
- evaluate_through_filtration(af, AlphaSimplex2D::DimensionExtractor())));
+ evaluate_through_map(m, AlphaSimplex2D::AlphaValueEvaluator()),
+ evaluate_through_map(m, AlphaSimplex2D::DimensionExtractor()));
#if 1
std::cout << 0 << std::endl << dgms[0] << std::endl;
--- a/examples/alphashapes/alphashapes2d.h Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/alphashapes/alphashapes2d.h Sun Dec 20 08:27:01 2009 -0800
@@ -74,7 +74,8 @@
typedef std::vector<AlphaSimplex2D> AlphaSimplex2DVector;
void fill_simplex_set(const Delaunay2D& Dt, AlphaSimplex2D::SimplexSet& simplices);
-void fill_complex(const Delaunay2D& Dt, AlphaSimplex2DVector& alpha_order);
+template<class Filtration>
+void fill_complex(const Delaunay2D& Dt, Filtration& filtration);
std::ostream& operator<<(std::ostream& out, const AlphaSimplex2D& s) { return s.operator<<(out); }
--- a/examples/alphashapes/alphashapes2d.hpp Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/alphashapes/alphashapes2d.hpp Sun Dec 20 08:27:01 2009 -0800
@@ -1,4 +1,5 @@
#include <utilities/log.h>
+#include <boost/foreach.hpp>
AlphaSimplex2D::
AlphaSimplex2D(const Delaunay2D::Vertex& v): alpha_(0), attached_(false)
@@ -104,16 +105,14 @@
rInfo("Vertices inserted");
}
-void fill_complex(const Delaunay2D& Dt, AlphaSimplex2DVector& alpha_order)
+template<class Filtration>
+void fill_complex(const Delaunay2D& Dt, Filtration& filtration)
{
// Compute all simplices with their alpha values and attachment information
+ // TODO: this can be optimized; the new Filtration can act as a SimplexSet
AlphaSimplex2D::SimplexSet simplices;
fill_simplex_set(Dt, simplices);
-
- // Sort simplices by their alpha values
- alpha_order.resize(simplices.size());
- std::copy(simplices.begin(), simplices.end(), alpha_order.begin());
- // std::sort(alpha_order.begin(), alpha_order.end(), AlphaSimplex2D::AlphaOrder());
- std::sort(alpha_order.begin(), alpha_order.end(), AlphaSimplex2D::VertexComparison());
+ BOOST_FOREACH(const AlphaSimplex2D& s, simplices)
+ filtration.push_back(s);
}
--- a/examples/alphashapes/alphashapes3d-cohomology.cpp Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/alphashapes/alphashapes3d-cohomology.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -75,6 +75,7 @@
rInfo("Delaunay triangulation computed");
// Set up the alpha shape filtration
+ typedef std::vector<AlphaSimplex3D> AlphaSimplex3DVector;
AlphaSimplex3DVector complex;
fill_complex(Dt, complex);
rInfo("Simplices: %d", complex.size());
--- a/examples/alphashapes/alphashapes3d.cpp Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/alphashapes/alphashapes3d.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -14,7 +14,7 @@
#include <boost/program_options.hpp>
-typedef Filtration<AlphaSimplex3DVector> AlphaFiltration;
+typedef Filtration<AlphaSimplex3D> AlphaFiltration;
typedef StaticPersistence<> Persistence;
typedef PersistenceDiagram<> PDgm;
@@ -72,12 +72,12 @@
}
rInfo("Delaunay triangulation computed");
- AlphaSimplex3DVector complex;
- fill_complex(Dt, complex);
- rInfo("Simplices: %d", complex.size());
+ AlphaFiltration af;
+ fill_complex(Dt, af);
+ rInfo("Simplices: %d", af.size());
// Create the alpha-shape filtration
- AlphaFiltration af(complex.begin(), complex.end(), AlphaSimplex3D::AlphaOrder());
+ af.sort(AlphaSimplex3D::AlphaOrder());
rInfo("Filtration initialized");
Persistence p(af);
@@ -89,12 +89,11 @@
rInfo("Simplices paired");
persistence_timer.check("Persistence timer");
+ Persistence::SimplexMap<AlphaFiltration> m = p.make_simplex_map(af);
std::map<Dimension, PDgm> dgms;
init_diagrams(dgms, p.begin(), p.end(),
- evaluate_through_map(OffsetMap<Persistence::OrderIndex, AlphaFiltration::Index>(p.begin(), af.begin()),
- evaluate_through_filtration(af, AlphaSimplex3D::AlphaValueEvaluator())),
- evaluate_through_map(OffsetMap<Persistence::OrderIndex, AlphaFiltration::Index>(p.begin(), af.begin()),
- evaluate_through_filtration(af, AlphaSimplex3D::DimensionExtractor())));
+ evaluate_through_map(m, AlphaSimplex3D::AlphaValueEvaluator()),
+ evaluate_through_map(m, AlphaSimplex3D::DimensionExtractor()));
#if 0
std::cout << 0 << std::endl << dgms[0] << std::endl;
std::cout << 1 << std::endl << dgms[1] << std::endl;
--- a/examples/alphashapes/alphashapes3d.h Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/alphashapes/alphashapes3d.h Sun Dec 20 08:27:01 2009 -0800
@@ -77,9 +77,9 @@
bool attached_;
};
-typedef std::vector<AlphaSimplex3D> AlphaSimplex3DVector;
void fill_simplex_set(const Delaunay3D& Dt, AlphaSimplex3D::SimplexSet& simplices);
-void fill_complex(const Delaunay3D& Dt, AlphaSimplex3DVector& alpha_order);
+template<class Filtration>
+void fill_complex(const Delaunay3D& Dt, Filtration& filtration);
std::ostream& operator<<(std::ostream& out, const AlphaSimplex3D& s) { return s.operator<<(out); }
--- a/examples/alphashapes/alphashapes3d.hpp Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/alphashapes/alphashapes3d.hpp Sun Dec 20 08:27:01 2009 -0800
@@ -1,4 +1,5 @@
#include <utilities/log.h>
+#include <boost/foreach.hpp>
AlphaSimplex3D::
AlphaSimplex3D(const Delaunay3D::Vertex& v): alpha_(0), attached_(false)
@@ -159,15 +160,12 @@
rInfo("Vertices inserted");
}
-void fill_complex(const Delaunay3D& Dt, AlphaSimplex3DVector& alpha_order)
+template<class Filtration>
+void fill_complex(const Delaunay3D& Dt, Filtration& filtration)
{
AlphaSimplex3D::SimplexSet simplices;
fill_simplex_set(Dt, simplices);
-
- // Sort simplices by their alpha values
- alpha_order.resize(simplices.size());
- std::copy(simplices.begin(), simplices.end(), alpha_order.begin());
- //std::sort(alpha_order.begin(), alpha_order.end(), AlphaSimplex3D::AlphaOrder());
- std::sort(alpha_order.begin(), alpha_order.end(), AlphaSimplex3D::VertexComparison());
+ BOOST_FOREACH(const AlphaSimplex3D& s, simplices)
+ filtration.push_back(s);
}
--- a/examples/cech-complex/cech-complex.cpp Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/cech-complex/cech-complex.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -15,8 +15,7 @@
typedef std::vector<Point> PointContainer;
typedef unsigned int PointIndex;
typedef Simplex<PointIndex, double> Smplx;
-typedef std::vector<Smplx> SimplexVector;
-typedef Filtration<SimplexVector> CechFiltration;
+typedef Filtration<Smplx> CechFiltration;
typedef StaticPersistence<> Persistence;
//typedef DynamicPersistenceTrails<> Persistence;
typedef PersistenceDiagram<> PDgm;
@@ -30,7 +29,7 @@
return numerator/denominator;
}
-void add_simplices(SimplexVector& sv, int d, const PointContainer& points)
+void add_simplices(CechFiltration& sv, int d, const PointContainer& points)
{
PointIndex indices[d+1];
for (int i = 0; i < d+1; ++i)
@@ -97,21 +96,18 @@
rInfo("Points read: %d", points.size());
// Compute simplices with their Cech values
- SimplexVector sv;
int num_simplices = 0;
for (int i = 0; i <= homology_d + 1; ++i)
num_simplices += choose(points.size(), i+1);
- sv.reserve(num_simplices);
rInfo("Reserved SimplexVector of size: %d", num_simplices);
+ CechFiltration cf;
for (int i = 0; i <= homology_d + 1; ++i)
- add_simplices(sv, i, points);
- rInfo("Size of SimplexVector: %d", sv.size());
-
- std::sort(sv.begin(), sv.end(), Smplx::VertexComparison());
+ add_simplices(cf, i, points);
+ rInfo("Size of SimplexVector: %d", cf.size());
- // Set up the filtration
- CechFiltration cf(sv.begin(), sv.end(), DataDimensionComparison<Smplx>());
+ // Sort the filtration
+ cf.sort(DataDimensionComparison<Smplx>());
rInfo("Filtration initialized");
// Compute persistence
@@ -120,12 +116,11 @@
p.pair_simplices();
rInfo("Simplices paired");
+ Persistence::SimplexMap<CechFiltration> m = p.make_simplex_map(cf);
std::map<Dimension, PDgm> dgms;
init_diagrams(dgms, p.begin(), p.end(),
- evaluate_through_map(OffsetMap<Persistence::OrderIndex, CechFiltration::Index>(p.begin(), cf.begin()),
- evaluate_through_filtration(cf, Smplx::DataEvaluator())),
- evaluate_through_map(OffsetMap<Persistence::OrderIndex, CechFiltration::Index>(p.begin(), cf.begin()),
- evaluate_through_filtration(cf, Smplx::DimensionExtractor())));
+ evaluate_through_map(m, Smplx::DataEvaluator()),
+ evaluate_through_map(m, Smplx::DimensionExtractor()));
for (int i = 0; i <= homology_d; ++i)
{
--- a/examples/fitness/avida-distance.cpp Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/fitness/avida-distance.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -10,7 +10,7 @@
typedef Simplex<AvidaOrganismDetail::IDType, double> Smplx;
typedef std::vector<Smplx> Complex;
-typedef Filtration<Complex> Fltr;
+typedef Filtration<Smplx> Fltr;
typedef StaticPersistence<> Persistence;
int main(int argc, char** argv)
@@ -55,7 +55,6 @@
simplices.back().add(organisms[j].id());
}
}
- std::sort(simplices.begin(), simplices.end(), Smplx::VertexComparison());
rInfo("Average distance: %f", float(avg_distance)/
((organisms.size()*organisms.size() - organisms.size())/2));
@@ -67,14 +66,14 @@
typedef std::vector<RealType> DeathVector;
DeathVector deaths;
Smplx::DataEvaluator eval;
- OffsetMap<Persistence::OrderIndex, Fltr::Index> m(p.begin(), f.begin());
- for (Persistence::OrderIndex i = p.begin(); i != p.end(); ++i)
+ Persistence::SimplexMap<Fltr> m = p.make_simplex_map(f);
+ for (Persistence::iterator i = p.begin(); i != p.end(); ++i)
{
- if (i == i->pair) continue;
+ if (i->unpaired()) continue;
if (i->sign())
{
- const Smplx& s = f.simplex(m[i]);
- const Smplx& t = f.simplex(m[i->pair]);
+ const Smplx& s = m[i];
+ const Smplx& t = m[i->pair];
AssertMsg(s.dimension() == 0, "Expecting only 0-dimensional diagram");
AssertMsg(eval(s) == 0, "Expecting only 0 birth values in 0-D diagram ");
deaths.push_back(eval(t));
--- a/examples/fitness/avida-rips-distance.cpp Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/fitness/avida-rips-distance.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -13,9 +13,8 @@
typedef ExplicitDistances<AvidaPopulationDetail> ExplicitDist;
typedef Rips<ExplicitDist> RipsGen;
typedef RipsGen::Simplex Smplx;
-typedef std::vector<Smplx> Complex;
-typedef Filtration<Complex, unsigned> Fltr;
+typedef Filtration<Smplx> Fltr;
typedef StaticPersistence<> Persistence;
int main(int argc, char** argv)
@@ -51,12 +50,11 @@
*/
rInfo("Starting to generate rips complex");
- Complex c;
- rips.generate(1, rips.max_distance()/2, make_push_back_functor(c));
- std::sort(c.begin(), c.end(), Smplx::VertexComparison());
+ Fltr f;
+ rips.generate(1, rips.max_distance()/2, make_push_back_functor(f));
rInfo("Generated Rips complex, filling filtration");
- Fltr f(c.begin(), c.end(), RipsGen::Comparison(rips.distances()));
+ f.sort(RipsGen::Comparison(rips.distances()));
Persistence p(f);
p.pair_simplices();
@@ -64,14 +62,14 @@
std::cout << "Outputting histogram of death values" << std::endl;
typedef std::vector<RealType> DeathVector;
DeathVector deaths;
- OffsetMap<Persistence::OrderIndex, Fltr::Index> m(p.begin(), f.begin());
- for (Persistence::OrderIndex i = p.begin(); i != p.end(); ++i)
+ Persistence::SimplexMap<Fltr> m = p.make_simplex_map(f);
+ for (Persistence::iterator i = p.begin(); i != p.end(); ++i)
{
- if (i == i->pair) continue;
+ if (i->unpaired()) continue;
if (i->sign())
{
- const Smplx& s = f.simplex(m[i]);
- const Smplx& t = f.simplex(m[i->pair]);
+ const Smplx& s = m[i];
+ const Smplx& t = m[i->pair];
AssertMsg(s.dimension() == 0, "Expecting only 0-dimensional diagram");
AssertMsg(evaluator(s) == 0, "Expecting only 0 birth values in 0-D diagram ");
deaths.push_back(evaluator(t));
--- a/examples/grid/CMakeLists.txt Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/grid/CMakeLists.txt Sun Dec 20 08:27:01 2009 -0800
@@ -1,5 +1,6 @@
set (targets
test-grid2D
+ test-grid2D-vineyard
combustion-vineyard)
if (use_dsrpdb)
@@ -7,9 +8,17 @@
pdbdistance-vineyard)
endif (use_dsrpdb)
-
-add_definitions (${cgal_cxxflags})
-foreach (t ${targets})
- add_executable (${t} ${t}.cpp)
- target_link_libraries (${t} ${libraries} ${cgal_libraries} ${dsrpdb_LIBRARY})
-endforeach (t ${targets})
+if (CGAL_FOUND)
+ # add_definitions (${CGAL_CXX_FLAGS_INIT})
+ add_definitions (-frounding-math)
+ include_directories (${CGAL_INCLUDE_DIRS})
+
+ foreach (t ${targets})
+ add_executable (${t} ${t}.cpp)
+ target_link_libraries (${t} ${libraries} ${CGAL_LIBRARY} ${dsrpdb_LIBRARY})
+ endforeach (t ${targets})
+
+ # Extra special for program_options
+ add_executable (pl-vineyard pl-vineyard.cpp)
+ target_link_libraries (pl-vineyard ${libraries} ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+endif (CGAL_FOUND)
--- a/examples/grid/combustion-vineyard.cpp Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/grid/combustion-vineyard.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -10,72 +10,79 @@
#include <algorithm>
#include "grid2D.h"
-#include "grid2Dvineyard.h"
+#include "lsvineyard.h"
const int xsize = 600;
const int ysize = 600;
-const int var = 0; // which variable to use out of nc of them in each record in the file
+const int var = 0; // which variable to use out of nc of them in each record in the file
template<typename T>
void read(std::ifstream& ifs, T& var)
{
- ifs.read(reinterpret_cast<char*>(&var), sizeof(T));
+ ifs.read(reinterpret_cast<char*>(&var), sizeof(T));
}
int main(int argc, char** argv)
{
- if (argc < 3)
- {
- std::cout << "Usage: combustion-vineyard FRAME1 FRAME2" << std::endl;
- exit(0);
- }
-
- int size0, nc0;
- int size1, nc1;
+ if (argc < 3)
+ {
+ std::cout << "Usage: combustion-vineyard FRAME1 FRAME2" << std::endl;
+ exit(0);
+ }
+
+ int size0, nc0;
+ int size1, nc1;
- std::cout << "Reading: " << argv[1] << std::endl;
- std::ifstream ifs0(argv[1], std::ios::binary);
- std::cout << "Reading: " << argv[2] << std::endl;
- std::ifstream ifs1(argv[2], std::ios::binary);
+ std::cout << "Reading: " << argv[1] << std::endl;
+ std::ifstream ifs0(argv[1], std::ios::binary);
+ std::cout << "Reading: " << argv[2] << std::endl;
+ std::ifstream ifs1(argv[2], std::ios::binary);
- if (!ifs0 || !ifs1)
- {
- std::cout << "Could not open the frames" << std::endl;
- exit(0);
- }
+ if (!ifs0 || !ifs1)
+ {
+ std::cout << "Could not open the frames" << std::endl;
+ exit(0);
+ }
- read(ifs0, size0); read(ifs0, nc0);
- read(ifs1, size1); read(ifs1, nc1);
+ read(ifs0, size0); read(ifs0, nc0);
+ read(ifs1, size1); read(ifs1, nc1);
- assert(size0 == size1); assert(nc0 == nc1);
- assert(size0 == xsize*ysize);
-
- Grid2D g0(xsize, ysize), g1(xsize, ysize);
-
- for (int y = 0; y < ysize; ++y)
- for (int x = 0; x < xsize; ++x)
- for (int d = 0; d < nc0; ++d)
- {
- float val0, val1;
- read(ifs0, val0);
- read(ifs1, val1);
- if (d == var)
- {
- g0(x,y) = val0;
- g1(x,y) = val1;
- }
- }
- std::cout << "Grids read" << std::endl;
-
- // Generate filtration and compute pairing
- Grid2DVineyard v(&g0);
- std::cout << "Filtration generated, size: " << v.filtration()->size() << std::endl;
- v.compute_pairing();
- std::cout << "Pairing computed" << std::endl;
-
- // Compute vineyard
- v.compute_vineyard(&g1, true);
- std::cout << "Vineyard computed" << std::endl;
+ assert(size0 == size1); assert(nc0 == nc1);
+ assert(size0 == xsize*ysize);
+
+ Grid2D g0(xsize, ysize), g1(xsize, ysize);
+
+ for (int y = 0; y < ysize; ++y)
+ for (int x = 0; x < xsize; ++x)
+ for (int d = 0; d < nc0; ++d)
+ {
+ float val0, val1;
+ read(ifs0, val0);
+ read(ifs1, val1);
+ if (d == var)
+ {
+ g0(x,y) = val0;
+ g1(x,y) = val1;
+ }
+ }
+ std::cout << "Grids read" << std::endl;
+
+ // Generate the complex, initialize the vineyard (which also computes the pairing)
+ typedef LSVineyard<Grid2D::CoordinateIndex, Grid2D> Grid2DVineyard;
- v.vineyard()->save_edges("combustion");
+ Grid2DVineyard::LSFiltration simplices;
+ g0.complex_generator(make_push_back_functor(simplices));
+ Grid2DVineyard::VertexComparison vcmp(g0);
+ Grid2DVineyard::SimplexComparison scmp(vcmp);
+ simplices.sort(scmp);
+
+ Grid2DVineyard v(g0.begin(), g0.end(), simplices, g0);
+ std::cout << "Filtration generated, size: " << v.filtration().size() << std::endl;
+ std::cout << "Pairing computed" << std::endl;
+
+ // Compute vineyard
+ v.compute_vineyard(g1, true);
+ std::cout << "Vineyard computed" << std::endl;
+
+ v.vineyard().save_edges("combustion");
}
--- a/examples/grid/grid2D.h Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/grid/grid2D.h Sun Dec 20 08:27:01 2009 -0800
@@ -18,67 +18,79 @@
#include <boost/serialization/base_object.hpp>
#include <boost/serialization/split_member.hpp>
#include <boost/serialization/nvp.hpp>
+#include <boost/serialization/export.hpp>
+
+#include <boost/iterator/counting_iterator.hpp>
#include "utilities/types.h"
-
-#include <boost/serialization/export.hpp>
+#include "topology/simplex.h"
/**
* Grid2D stores a grid
*/
class Grid2D
{
- public:
- typedef RealType ValueType;
- typedef unsigned int CoordinateIndex;
-
- typedef std::vector<ValueType> ValueVector;
+ public:
+ typedef RealType ValueType;
+ typedef ValueType result_type;
+
+ typedef unsigned int CoordinateIndex;
+ typedef boost::counting_iterator<CoordinateIndex> iterator;
+
+ typedef Simplex<CoordinateIndex> Smplx;
+ typedef std::vector<ValueType> ValueVector;
- public:
- Grid2D(CoordinateIndex xx = 1, CoordinateIndex yy = 1);
+ public:
+ Grid2D(CoordinateIndex xx = 1, CoordinateIndex yy = 1);
- /// Sets the grid dimensions to (xx,yy)
- void change_dimensions(CoordinateIndex xx, CoordinateIndex yy);
+ /// Sets the grid dimensions to (xx,yy)
+ void change_dimensions(CoordinateIndex xx, CoordinateIndex yy);
- ValueType& operator()(CoordinateIndex i, CoordinateIndex j) { return data[i*x + j]; }
- const ValueType& operator()(CoordinateIndex i, CoordinateIndex j) const { return data[i*x + j]; }
- ValueType& operator()(CoordinateIndex i) { return data[i]; }
- const ValueType& operator()(CoordinateIndex i) const { return data[i]; }
+ ValueType& operator()(CoordinateIndex i, CoordinateIndex j) { return data[i*x + j]; }
+ const ValueType& operator()(CoordinateIndex i, CoordinateIndex j) const { return data[i*x + j]; }
+ ValueType& operator()(CoordinateIndex i) { return data[i]; }
+ const ValueType& operator()(CoordinateIndex i) const { return data[i]; }
- CoordinateIndex xsize() const { return x; }
- CoordinateIndex ysize() const { return y; }
- CoordinateIndex size() const { return x*y; }
-
- /* Given a sequential index of an element return its coordinates */
- CoordinateIndex xpos(CoordinateIndex i) const { return i / x; }
- CoordinateIndex ypos(CoordinateIndex i) const { return i % x; }
- CoordinateIndex seq(CoordinateIndex i, CoordinateIndex j) const;
+ CoordinateIndex xsize() const { return x; }
+ CoordinateIndex ysize() const { return y; }
+ CoordinateIndex size() const { return x*y; }
- std::ostream& operator<<(std::ostream& out) const;
-
- static const CoordinateIndex INVALID_INDEX = -1;
+ iterator begin() const { return iterator(0); }
+ iterator end() const { return iterator(size()); }
+
+ /* Given a sequential index of an element return its coordinates */
+ CoordinateIndex xpos(CoordinateIndex i) const { return i / x; }
+ CoordinateIndex ypos(CoordinateIndex i) const { return i % x; }
+ CoordinateIndex seq(CoordinateIndex i, CoordinateIndex j) const;
- private:
- CoordinateIndex x,y;
- ValueVector data;
+ template<class Functor>
+ void complex_generator(const Functor& f);
+
+ std::ostream& operator<<(std::ostream& out) const;
+
+ static const CoordinateIndex INVALID_INDEX = -1;
+
+ private:
+ CoordinateIndex x,y;
+ ValueVector data;
#if 0
- private:
- // Serialization
- friend class boost::serialization::access;
+ private:
+ // Serialization
+ friend class boost::serialization::access;
- template<class Archive> void save(Archive& ar, version_type ) const;
- template<class Archive> void load(Archive& ar, version_type );
+ template<class Archive> void save(Archive& ar, version_type ) const;
+ template<class Archive> void load(Archive& ar, version_type );
- BOOST_SERIALIZATION_SPLIT_MEMBER()
+ BOOST_SERIALIZATION_SPLIT_MEMBER()
#endif
};
//BOOST_CLASS_EXPORT(Grid2D)
-
+
-std::ostream& operator<<(std::ostream& out, const Grid2D& grid) { return grid.operator<<(out); }
+std::ostream& operator<<(std::ostream& out, const Grid2D& grid) { return grid.operator<<(out); }
-
+
#include "grid2D.hpp"
#endif // __GRID2D_H__
--- a/examples/grid/grid2D.hpp Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/grid/grid2D.hpp Sun Dec 20 08:27:01 2009 -0800
@@ -5,40 +5,40 @@
Grid2D::
Grid2D(CoordinateIndex xx, CoordinateIndex yy):
- x(xx), y(yy), data(x*y)
+ x(xx), y(yy), data(x*y)
{}
-void
+void
Grid2D::
change_dimensions(CoordinateIndex xx, CoordinateIndex yy)
{
- x = xx; y = yy;
- data.resize(x*y);
+ x = xx; y = yy;
+ data.resize(x*y);
}
Grid2D::CoordinateIndex
Grid2D::
seq(CoordinateIndex i, CoordinateIndex j) const
{
- // Do not forget to check if less than 0, if Index is made signed --- dangerous
- if (i >= x || j >= y)
- return INVALID_INDEX;
+ // Do not forget to check if less than 0, if Index is made signed --- dangerous
+ if (i >= x || j >= y)
+ return INVALID_INDEX;
- return i*x + j;
+ return i*x + j;
}
-std::ostream&
+std::ostream&
Grid2D::
operator<<(std::ostream& out) const
{
- for (Grid2D::CoordinateIndex i = 0; i < xsize(); ++i)
- {
- for (Grid2D::CoordinateIndex j = 0; j < ysize(); ++j)
- std::cout << operator()(i, j) << ' ';
- std::cout << std::endl;
- }
- return out;
-}
+ for (Grid2D::CoordinateIndex i = 0; i < xsize(); ++i)
+ {
+ for (Grid2D::CoordinateIndex j = 0; j < ysize(); ++j)
+ std::cout << operator()(i, j) << ' ';
+ std::cout << std::endl;
+ }
+ return out;
+}
#if 0
using boost::serialization::make_nvp;
@@ -48,18 +48,74 @@
Grid2D::
save(Archive& ar, version_type ) const
{
- ar << BOOST_SERIALIZATION_NVP(x);
- ar << BOOST_SERIALIZATION_NVP(y);
- ar << make_nvp("data", data);
+ ar << BOOST_SERIALIZATION_NVP(x);
+ ar << BOOST_SERIALIZATION_NVP(y);
+ ar << make_nvp("data", data);
}
-template<class Archive>
+template<class Archive>
void
Grid2D::
load(Archive& ar, version_type )
{
- ar >> make_nvp("x", x);
- ar >> make_nvp("y", y);
- ar >> make_nvp("data", data);
+ ar >> make_nvp("x", x);
+ ar >> make_nvp("y", y);
+ ar >> make_nvp("data", data);
}
#endif
+
+template<class Functor>
+void
+Grid2D::
+complex_generator(const Functor& f)
+{
+ for (CoordinateIndex x = 0; x < xsize() - 1; ++x)
+ for (CoordinateIndex y = 0; y < ysize() - 1; ++y)
+ {
+ CoordinateIndex v(seq(x,y));
+ CoordinateIndex vh(seq(x+1,y));
+ CoordinateIndex vv(seq(x,y+1));
+ CoordinateIndex vd(seq(x+1,y+1));
+
+ Smplx sh; sh.add(v);
+ f(sh);
+ sh.add(vh); f(sh); // Horizontal edge
+ sh.add(vd); f(sh); // "Horizontal" triangle
+
+ Smplx sv; sv.add(v);
+ sv.add(vv); f(sv); // Vertical edge
+ sv.add(vd); f(sv); // "Vertical" triangle
+
+ Smplx sd; sd.add(v);
+ sd.add(vd); f(sd); // Diagonal edge
+
+ if (y == ysize() - 2)
+ {
+ Smplx s; s.add(vv);
+ s.add(vd); f(s); // Top edge
+ }
+ if (x == xsize() - 2)
+ {
+ Smplx s; s.add(vh);
+ s.add(vd); f(s); // Right edge
+ }
+ }
+
+ // Last row
+ for (CoordinateIndex x = 0; x < xsize(); ++x)
+ {
+ std::cout << x << " " << ysize() - 1 << " " << seq(x, ysize() - 1) << std::endl;
+ CoordinateIndex v(seq(x, ysize() - 1));
+ Smplx s; s.add(v);
+ f(s);
+ }
+
+ // Last column
+ for (CoordinateIndex y = 0; y < ysize() - 1; ++y)
+ {
+ std::cout << xsize() - 1 << " " << y << " " << seq(xsize() - 1, y) << std::endl;
+ CoordinateIndex v(seq(xsize() - 1, y));
+ Smplx s; s.add(v);
+ f(s);
+ }
+}
--- a/examples/grid/grid2Dvineyard.h Thu Dec 03 13:21:23 2009 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,190 +0,0 @@
-/*
- * Author: Dmitriy Morozov
- * Department of Computer Science, Duke University, 2005 -- 2008
- */
-
-#ifndef __GRID2DVINEYARD_H__
-#define __GRID2DVINEYARD_H__
-
-#include "grid2D.h"
-#include "topology/lowerstarfiltration.h"
-
-#include <CGAL/Kinetic/Inexact_simulation_traits.h>
-#include <CGAL/Kinetic/Sort.h>
-#include <CGAL/Kinetic/Sort_visitor_base.h>
-
-#include <vector>
-
-
-class Grid2DVineyard
-{
- public:
- typedef Grid2DVineyard Self;
-
- typedef Grid2D::CoordinateIndex CoordinateIndex;
- typedef Grid2D::ValueType ValueType;
-
- class KineticVertexType;
- typedef std::vector<KineticVertexType> VertexVector;
- typedef VertexVector::iterator VertexIndex;
-
- typedef LowerStarFiltration<VertexIndex> LSFiltration;
-
- class StaticEvaluator;
- class KineticEvaluator;
- class VertexComparison;
-
- typedef LSFiltration::Index Index;
- typedef LSFiltration::Simplex Simplex;
- typedef LSFiltration::VertexOrderIndex VertexOrderIndex;
- typedef LSFiltration::VertexType<CoordinateIndex> LSVertexType;
-
- typedef LSFiltration::Vineyard Vineyard;
- typedef Vineyard::Evaluator Evaluator;
-
- class SortVisitor;
- typedef CGAL::Kinetic::Inexact_simulation_traits Traits;
- typedef CGAL::Kinetic::Sort<Traits, SortVisitor> Sort;
- typedef Traits::Simulator Simulator;
- typedef Traits::Active_points_1_table ActivePointsTable;
- typedef ActivePointsTable::Key Key;
- typedef std::map<Key, VertexOrderIndex> KeyOrderMap;
-
- typedef std::vector<Grid2D*> GridStackContainer;
-
- public:
- Grid2DVineyard(Grid2D* g);
- ~Grid2DVineyard();
-
- void compute_pairing();
- void compute_vineyard(Grid2D* grid, bool explicit_events = false);
-
- Grid2D* grid() const { return grid_stack_.back(); }
- Grid2D* grid(int i) const { return grid_stack_[i]; }
- int num_grids() const { return grid_stack_.size(); }
- const LSFiltration* filtration() const { return filtration_; }
- const Vineyard* vineyard() const { return vineyard_; }
-
- public:
- // For Kinetic Sort
- void swap(Key a, Key b);
-
- protected:
- // Do something cleverer
- virtual bool neighbors(VertexIndex v1, VertexIndex v2) const { return true; }
-
- private:
- void add_simplices();
- void change_evaluator(Evaluator* eval);
-
- private:
- GridStackContainer grid_stack_;
- VertexVector vertices_;
- LSFiltration* filtration_;
- Vineyard* vineyard_;
- Evaluator* evaluator_;
-
- KeyOrderMap kinetic_map_;
-
-#if 0
- private:
- // Serialization
- friend class boost::serialization::access;
-
- Grid2DVineyard() {}
-
- template<class Archive>
- void serialize(Archive& ar, version_type )
- {
- ar & BOOST_SERIALIZATION_NVP(grid_stack_);
- ar & BOOST_SERIALIZATION_NVP(vertices_);
- ar & BOOST_SERIALIZATION_NVP(filtration_);
- };
-#endif
-};
-
-//BOOST_CLASS_EXPORT(Grid2DVineyard)
-
-class Grid2DVineyard::KineticVertexType: public LSVertexType
-{
- public:
- typedef LSVertexType Parent;
-
- Key kinetic_key() const { return key_; }
- void set_kinetic_key(Key k) { key_ = k; }
-
- private:
- Key key_;
-};
-
-std::ostream& operator<<(std::ostream& out, const Grid2DVineyard::VertexIndex& vi) { return out << vi->index(); }
-
-class Grid2DVineyard::VertexComparison
-{
- public:
- VertexComparison(const Grid2D* g): grid(g) {}
- bool operator()(VertexIndex i, VertexIndex j) const { return (*grid)(i->index()) <
- (*grid)(j->index()); }
-
- private:
- const Grid2D* grid;
-
-#if 0
- private:
- // Serialization
- friend class boost::serialization::access;
-
- VertexComparison() {}
-
- template<class Archive>
- void serialize(Archive& ar, version_type ) { ar & BOOST_SERIALIZATION_NVP(grid); }
-#endif
-};
-
-class Grid2DVineyard::StaticEvaluator: public Evaluator
-{
- public:
- StaticEvaluator(Grid2D* grid, RealType time):
- time_(time), grid_(grid) {}
-
- virtual RealType time() { return time_; }
- virtual RealType value(const Simplex& s) { return (*grid_)(s.get_attachment()->index()); }
-
- private:
- RealType time_;
- Grid2D* grid_;
-};
-
-class Grid2DVineyard::KineticEvaluator: public Evaluator
-{
- public:
- KineticEvaluator(Simulator::Handle sp, ActivePointsTable::Handle apt, RealType time_offset):
- sp_(sp), apt_(apt), time_offset_(time_offset) {}
-
- virtual RealType time() { return time_offset_ + CGAL::to_double(get_time()); }
- virtual RealType value(const Simplex& s) { return CGAL::to_double(apt_->at(s.get_attachment()->kinetic_key()).x()(get_time())); }
-
- private:
- Simulator::Time get_time() { return sp_->current_time(); }
-
- Simulator::Handle sp_;
- ActivePointsTable::Handle apt_;
- RealType time_offset_;
-};
-
-
-class Grid2DVineyard::SortVisitor: public CGAL::Kinetic::Sort_visitor_base
-{
- public:
- SortVisitor(Grid2DVineyard* gv): gv_(gv) {}
-
- template<class Vertex_handle>
- void before_swap(Vertex_handle a, Vertex_handle b) const { gv_->swap(*a,*b); }
-
- private:
- Grid2DVineyard* gv_;
-};
-
-#include "grid2Dvineyard.hpp"
-
-#endif // __GRID2DVINEYARD_H__
--- a/examples/grid/grid2Dvineyard.hpp Thu Dec 03 13:21:23 2009 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,147 +0,0 @@
-#include <utilities/log.h>
-
-/* Implementation */
-
-Grid2DVineyard::
-Grid2DVineyard(Grid2D* g): vertices_(g->size())
-{
- grid_stack_.push_back(g);
- for (CoordinateIndex i = 0; i < g->size(); ++i)
- vertices_[i].set_index(i);
-
- evaluator_ = new StaticEvaluator(grid(), 0);
- vineyard_ = new Vineyard(evaluator_);
-
- filtration_ = new LSFiltration(vertices_.begin(), vertices_.end(), VertexComparison(grid()), vineyard_);
- add_simplices();
-}
-
-Grid2DVineyard::
-~Grid2DVineyard()
-{
- delete filtration_;
- delete vineyard_;
- delete evaluator_;
-}
-
-void
-Grid2DVineyard::
-compute_pairing()
-{
- filtration_->fill_simplex_index_map();
- filtration_->pair_simplices(filtration_->begin(), filtration_->end());
- vineyard_->start_vines(filtration_->begin(), filtration_->end());
-}
-
-void
-Grid2DVineyard::
-compute_vineyard(Grid2D* g, bool explicit_events)
-{
- AssertMsg(filtration_->is_paired(), "Simplices must be paired for a vineyard to be computed");
-
- typedef Traits::Kinetic_kernel::Point_1 Point;
- typedef Traits::Kinetic_kernel::Function_kernel::Construct_function CF;
- typedef Traits::Kinetic_kernel::Motion_function F;
-
- Traits tr(0,1);
- Simulator::Handle sp = tr.simulator_handle();
- ActivePointsTable::Handle apt = tr.active_points_1_table_handle();
- Sort sort(tr, SortVisitor(this));
-
- // Setup the (linear) trajectories
- std::cout << "Setting up trajectories" << std::endl;
- CF cf;
- kinetic_map_.clear();
- for (VertexIndex cur = vertices_.begin(); cur != vertices_.end(); ++cur)
- {
- ValueType val0 = (*grid())(cur->index());
- ValueType val1 = (*g)(cur->index());
- F x = cf(F::NT(val0), F::NT(val1 - val0)); // x = val0 + (val1 - val0)*t
- Point p(x);
- cur->set_kinetic_key(apt->insert(p));
- kinetic_map_[cur->kinetic_key()] = cur->get_order();
- if (cur->index() % 10000 == 0)
- std::cout << "Added trajectory: " << cur->index() << " " << val0 << " " << val1 << std::endl;
- }
-
- // Process all the events (compute the vineyard in the process)
- change_evaluator(new KineticEvaluator(sp, apt, num_grids() - 1));
- if (explicit_events)
- {
- while (sp->next_event_time() < 1)
- {
- std::cout << "Next event time: " << sp->next_event_time() << std::endl;
- sp->set_current_event_number(sp->current_event_number() + 1);
- std::cout << "Processed event" << std::endl;
- }
- } else
- sp->set_current_time(1.0);
- std::cout << "Processed " << sp->current_event_number() << " events" << std::endl;
-
- // Add the grid to the stack
- grid_stack_.push_back(g);
- change_evaluator(new StaticEvaluator(grid(), num_grids() - 1));
- vineyard_->record_diagram(filtration_->begin(), filtration_->end());
-}
-
-void
-Grid2DVineyard::
-swap(Key a, Key b)
-{
- VertexOrderIndex ao = kinetic_map_[a], bo = kinetic_map_[b];
- AssertMsg(filtration_->get_vertex_cmp()(ao, bo), "In swap(a,b), a must precede b");
- filtration_->transpose_vertices(ao);
- AssertMsg(filtration_->get_vertex_cmp()(bo, ao), "In swap(a,b), b must precede a after the transposition");
-}
-
-void
-Grid2DVineyard::
-add_simplices()
-{
- // Takes advantage of LowerStarFiltration's smart append (which allows faces
- // to be inserted after cofaces, since everything is rearranged in the
- // proper lower star order anyway). Also note that vertices were added by
- // LowerStarFiltration's constructor
- for (CoordinateIndex x = 0; x < grid()->xsize() - 1; ++x)
- for (CoordinateIndex y = 0; y < grid()->ysize() - 1; ++y)
- {
- VertexIndex v(&vertices_[grid()->seq(x,y)]);
- VertexIndex vh(&vertices_[grid()->seq(x+1,y)]);
- VertexIndex vv(&vertices_[grid()->seq(x,y+1)]);
- VertexIndex vd(&vertices_[grid()->seq(x+1,y+1)]);
-
- Simplex sh(2, v);
- sh.add(vh); filtration_->append(sh); // Horizontal edge
- sh.add(vd); filtration_->append(sh); // "Horizontal" triangle
-
- Simplex sv(2, v);
- sv.add(vv); filtration_->append(sv); // Vertical edge
- sv.add(vd); filtration_->append(sv); // "Vertical" triangle
-
- Simplex sd(2, v);
- sd.add(vd); filtration_->append(sd); // Diagonal edge
-
- if (y == grid()->ysize() - 2)
- {
- Simplex s(1, vv);
- s.add(vd); filtration_->append(s); // Top edge
- }
- if (x == grid()->xsize() - 2)
- {
- Simplex s(1, vh);
- s.add(vd); filtration_->append(s); // Right edge
- }
- }
-}
-
-void
-Grid2DVineyard::
-change_evaluator(Evaluator* eval)
-{
- AssertMsg(evaluator_ != 0, "change_evaluator() assumes that existing evaluator is not null");
-
- delete evaluator_;
- evaluator_ = eval;
- vineyard_->set_evaluator(evaluator_);
-}
-
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/grid/lsvineyard.h Sun Dec 20 08:27:01 2009 -0800
@@ -0,0 +1,245 @@
+/**
+ * Author: Dmitriy Morozov
+ * Department of Computer Science, Duke University, 2005 -- 2009
+ */
+
+#ifndef __LSVINEYARD_H__
+#define __LSVINEYARD_H__
+
+#include <iostream>
+
+#include "topology/simplex.h"
+#include "topology/dynamic-persistence.h"
+#include "topology/lowerstarfiltration.h"
+#include "topology/vineyard.h"
+
+#include <utilities/indirect.h>
+
+#include <CGAL/Kinetic/Inexact_simulation_traits.h>
+#include <CGAL/Kinetic/Sort.h>
+#include <CGAL/Kinetic/Sort_visitor_base.h>
+
+#include <list>
+
+
+template<class Vertex_, class VertexEvaluator_, class Simplex_ = Simplex<Vertex_>, class Filtration_ = Filtration<Simplex_> >
+class LSVineyard
+{
+ public:
+ typedef LSVineyard Self;
+
+ typedef Vertex_ Vertex;
+ typedef VertexEvaluator_ VertexEvaluator;
+ typedef typename VertexEvaluator::result_type VertexValue;
+
+ typedef Simplex_ Simplex;
+ typedef Filtration_ LSFiltration;
+ typedef typename LSFiltration::Index LSFIndex;
+
+ class KineticVertexType;
+ class KineticVertexComparison;
+ typedef std::list<KineticVertexType> VertexContainer;
+ typedef typename VertexContainer::iterator VertexIndex;
+
+ struct AttachmentData: public VineData
+ {
+ void set_attachment(VertexIndex v) { attachment = v; }
+ VertexIndex attachment;
+ };
+ typedef DynamicPersistenceTrails<AttachmentData> Persistence;
+ typedef typename Persistence::OrderIndex Index;
+ typedef typename Persistence::iterator iterator;
+
+ typedef typename Persistence::template SimplexMap<LSFiltration>
+ PFMap;
+
+ class Evaluator;
+ class StaticEvaluator;
+ class KineticEvaluator;
+ class DimensionFromIterator;
+
+ typedef ThroughEvaluatorComparison<VertexEvaluator> VertexComparison;
+ typedef MaxVertexComparison<Simplex, VertexComparison> SimplexComparison;
+
+ class TranspositionVisitor;
+ friend class TranspositionVisitor;
+
+ typedef Vineyard<Index, iterator, Evaluator> Vnrd;
+
+ class SortVisitor;
+ typedef CGAL::Kinetic::Inexact_simulation_traits Traits;
+ typedef CGAL::Kinetic::Sort<Traits, SortVisitor> Sort;
+ typedef Traits::Simulator Simulator;
+ typedef Traits::Active_points_1_table ActivePointsTable;
+ typedef ActivePointsTable::Key Key;
+ typedef std::map<Key, VertexIndex> KeyVertexMap;
+
+ public:
+ template<class VertexIterator>
+ LSVineyard(VertexIterator begin, VertexIterator end,
+ LSFiltration& filtration,
+ const VertexEvaluator& veval = VertexEvaluator());
+ ~LSVineyard();
+
+ void compute_vineyard(const VertexEvaluator& veval, bool explicit_events = false);
+ bool transpose_vertices(VertexIndex vi);
+
+ const LSFiltration& filtration() const { return filtration_; }
+ const Vnrd& vineyard() const { return vineyard_; }
+ const Persistence& persistence() const { return persistence_; }
+ const VertexComparison& vertex_comparison() const { return vcmp_; }
+ const VertexEvaluator& vertex_evaluator() const { return veval_; }
+ const SimplexComparison& simplex_comparison() const { return scmp_; }
+
+ VertexValue vertex_value(const Vertex& v) const { return veval_(v); }
+ VertexValue simplex_value(const Simplex& s) const { return vertex_value(*std::max_element(s.vertices().begin(), s.vertices().end(), vcmp_)); }
+ const Simplex& pfmap(iterator i) const { return pfmap_[i]; }
+ const Simplex& pfmap(Index i) const { return pfmap_[i]; }
+
+ Index index(iterator i) const { return persistence_.index(i); }
+
+ public:
+ // For Kinetic Sort
+ void swap(Key a, Key b);
+
+ private:
+ void change_evaluator(Evaluator* eval);
+ void set_attachment(iterator i, VertexIndex vi) { persistence_.modifier()(i, boost::bind(&AttachmentData::set_attachment, bl::_1, vi)); }
+ void transpose_filtration(iterator i) { filtration_.transpose(filtration_.begin() + (i - persistence_.begin())); }
+
+ private:
+ VertexContainer vertices_;
+
+ VertexEvaluator veval_;
+ VertexComparison vcmp_;
+ SimplexComparison scmp_;
+
+ LSFiltration& filtration_;
+ Persistence persistence_;
+ PFMap pfmap_;
+
+ Vnrd vineyard_;
+ Evaluator* evaluator_;
+ unsigned time_count_;
+
+ KeyVertexMap kinetic_map_;
+
+#if 0
+ private:
+ // Serialization
+ friend class boost::serialization::access;
+
+ LSVineyard() {}
+
+ template<class Archive>
+ void serialize(Archive& ar, version_type )
+ {
+ ar & BOOST_SERIALIZATION_NVP(grid_stack_);
+ ar & BOOST_SERIALIZATION_NVP(vertices_);
+ ar & BOOST_SERIALIZATION_NVP(filtration_);
+ };
+#endif
+};
+
+//BOOST_CLASS_EXPORT(LSVineyard)
+
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::KineticVertexType
+{
+ public:
+ KineticVertexType(const Vertex& v):
+ vertex_(v) {}
+
+ Key kinetic_key() const { return key_; }
+ void set_kinetic_key(Key k) { key_ = k; }
+
+ Vertex vertex() const { return vertex_; }
+ void set_vertex(Vertex v) { vertex_ = v; }
+
+ iterator simplex_index() const { return simplex_index_; }
+ void set_simplex_index(iterator i) { simplex_index_ = i; }
+
+ private:
+ Key key_;
+ Vertex vertex_;
+ iterator simplex_index_;
+};
+
+template<class V, class VE, class S, class C>
+std::ostream&
+operator<<(std::ostream& out, const typename LSVineyard<V,VE,S,C>::VertexIndex& vi)
+{ return out << vi->vertex(); }
+
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::KineticVertexComparison: public std::binary_function<const KineticVertexType&, const KineticVertexType&, bool>
+{
+ public:
+ KineticVertexComparison(const VertexComparison& vcmp):
+ vcmp_(vcmp) {}
+
+ bool operator()(const KineticVertexType& v1, const KineticVertexType& v2) const
+ { return vcmp_(v1.vertex(), v2.vertex()); }
+
+ private:
+ VertexComparison vcmp_;
+};
+
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::TranspositionVisitor: public Persistence::TranspositionVisitor
+{
+ public:
+ typedef typename Persistence::TranspositionVisitor Parent;
+ typedef typename LSVineyard<V,VE,S,C>::iterator iterator;
+ typedef typename LSVineyard<V,VE,S,C>::Index Index;
+
+ TranspositionVisitor(LSVineyard& v): lsvineyard_(v) {}
+
+ void transpose(iterator i) { lsvineyard_.transpose_filtration(i); }
+ void switched(iterator i, SwitchType type) { lsvineyard_.vineyard_.switched(index(i), index(boost::next(i))); }
+
+ private:
+ Index index(iterator i) { return lsvineyard_.index(i); }
+
+ LSVineyard& lsvineyard_;
+};
+
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::Evaluator: public std::unary_function<Index, RealType>
+{
+ public:
+ virtual RealType time() const =0;
+ virtual RealType operator()(Index i) const =0;
+ virtual Dimension dimension(Index i) const =0;
+ virtual RealType operator()(iterator i) const { return operator()(&*i); }
+ virtual Dimension dimension(iterator i) const { return dimension(&*i); }
+};
+
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::SortVisitor: public CGAL::Kinetic::Sort_visitor_base
+{
+ public:
+ SortVisitor(LSVineyard& v):
+ vineyard_(v) {}
+
+ template<class Vertex_handle>
+ void pre_swap(Vertex_handle a, Vertex_handle b) const { vineyard_.swap(*a,*b); }
+
+ private:
+ LSVineyard& vineyard_;
+};
+
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::DimensionFromIterator: std::unary_function<iterator, Dimension>
+{
+ public:
+ DimensionFromIterator(const PFMap& pfmap): pfmap_(pfmap) {}
+
+ Dimension operator()(iterator i) const { return pfmap_[i].dimension(); }
+
+ private:
+ const PFMap& pfmap_;
+};
+
+#include "lsvineyard.hpp"
+
+#endif // __LSVINEYARD_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/grid/lsvineyard.hpp Sun Dec 20 08:27:01 2009 -0800
@@ -0,0 +1,254 @@
+#include <utilities/log.h>
+
+#ifdef LOGGING
+static rlog::RLogChannel* rlLSVineyard = DEF_CHANNEL("lsvineyard/info", rlog::Log_Debug);
+static rlog::RLogChannel* rlLSVineyardDebug = DEF_CHANNEL("lsvineyard/debug", rlog::Log_Debug);
+#endif // LOGGING
+
+#ifdef COUNTERS
+static Counter* cVertexTransposition = GetCounter("lsfiltration/transposition"); // counts number of vertex transpositions
+static Counter* cAttachment = GetCounter("lsfiltration/attachment"); // counts the number of attachment changes
+#endif
+
+
+template<class V, class VE, class S, class F>
+template<class VertexIterator>
+LSVineyard<V,VE,S,F>::
+LSVineyard(VertexIterator begin, VertexIterator end,
+ LSFiltration& filtration,
+ const VertexEvaluator& veval):
+ vertices_(begin, end), filtration_(filtration),
+ persistence_(filtration_),
+ veval_(veval), vcmp_(veval_), scmp_(vcmp_),
+ pfmap_(persistence_.make_simplex_map(filtration_)),
+ time_count_(0)
+{
+ // TODO: LSVineyard really should sort the filtration_ itself
+ // filtration_.sort(scmp_);
+ // persistence_.init(filtration_);
+
+ rLog(rlLSVineyardDebug, "Initializing LSVineyard");
+ persistence_.pair_simplices();
+ rLog(rlLSVineyardDebug, "Simplices paired");
+
+ vertices_.sort(KineticVertexComparison(vcmp_)); // sort vertices w.r.t. vcmp_
+#if LOGGING
+ rLog(rlLSVineyardDebug, "Vertex order:");
+ for (VertexIndex cur = vertices_.begin(); cur != vertices_.end(); ++cur)
+ rLog(rlLSVineyardDebug, " %d", cur->vertex());
+#endif
+
+ // Initialize simplex_index() and attachment
+ VertexIndex vi = boost::prior(vertices_.begin());
+ for (iterator i = persistence_.begin(); i != persistence_.end(); ++i)
+ {
+ const Simplex& s = pfmap_[i];
+ if (s.dimension() == 0)
+ {
+ ++vi;
+ AssertMsg(s.vertices().front() == vi->vertex(), "In constructor, simplices and vertices must match.");
+ vi->set_simplex_index(i);
+ }
+ set_attachment(i, vi);
+ rLog(rlLSVineyardDebug, "%s attached to %d", tostring(pfmap(i)).c_str(), vi->vertex());
+ }
+
+ evaluator_ = new StaticEvaluator(*this, time_count_);
+ vineyard_.set_evaluator(evaluator_);
+ vineyard_.start_vines(persistence_.begin(), persistence_.end());
+}
+
+template<class V, class VE, class S, class F>
+LSVineyard<V,VE,S,F>::
+~LSVineyard()
+{
+ delete evaluator_;
+}
+
+template<class V, class VE, class S, class F_>
+void
+LSVineyard<V,VE,S,F_>::
+compute_vineyard(const VertexEvaluator& veval, bool explicit_events)
+{
+ typedef Traits::Kinetic_kernel::Point_1 Point;
+ typedef Traits::Kinetic_kernel::Function_kernel::Construct_function CF;
+ typedef Traits::Kinetic_kernel::Motion_function F;
+
+ Traits tr(0,1);
+ Simulator::Handle sp = tr.simulator_handle();
+ ActivePointsTable::Handle apt = tr.active_points_1_table_handle();
+ Sort sort(tr, SortVisitor(*this));
+
+ // Setup the (linear) trajectories
+ rLog(rlLSVineyard, "Setting up trajectories");
+ CF cf;
+ kinetic_map_.clear();
+ for (VertexIndex cur = vertices_.begin(); cur != vertices_.end(); ++cur)
+ {
+ VertexValue val0 = veval_(cur->vertex());
+ VertexValue val1 = veval(cur->vertex());
+ rLog(rlLSVineyardDebug, "Vertex %d: %f -> %f", cur->vertex(), val0, val1);
+ F x = cf(F::NT(val0), F::NT(val1 - val0)); // x = val0 + (val1 - val0)*t = (1-t)*val0 + t*val1
+ Point p(x);
+ cur->set_kinetic_key(apt->insert(p));
+ kinetic_map_[cur->kinetic_key()] = cur;
+ rLog(rlLSVineyardDebug, "Scheduling: %d %s", cur->vertex(), tostring(x).c_str());
+ }
+
+ // Process all the events (compute the vineyard in the process)
+ change_evaluator(new KineticEvaluator(*this, sp, apt, time_count_));
+ if (explicit_events)
+ {
+ while (sp->next_event_time() < 1)
+ {
+ rLog(rlLSVineyardDebug, "Next event time: %f", sp->next_event_time());
+ sp->set_current_event_number(sp->current_event_number() + 1);
+ rLog(rlLSVineyardDebug, "Processed event");
+ }
+ } else
+ sp->set_current_time(1.0);
+ rLog(rlLSVineyard, "Processed %d events", sp->current_event_number());
+
+ veval_ = veval;
+ change_evaluator(new StaticEvaluator(*this, ++time_count_));
+ vineyard_.record_diagram(persistence_.begin(), persistence_.end());
+}
+
+template<class V, class VE, class S, class F>
+void
+LSVineyard<V,VE,S,F>::
+swap(Key a, Key b)
+{
+ rLog(rlLSVineyardDebug, "Entered swap");
+ VertexIndex ao = kinetic_map_[a], bo = kinetic_map_[b];
+ AssertMsg(vcmp_(ao->vertex(), bo->vertex()), "In swap(a,b), a must precede b");
+ transpose_vertices(ao);
+ // AssertMsg(vcmp_(bo->vertex(), ao->vertex()), "In swap(a,b), b must precede a after the transposition");
+}
+
+template<class V, class VE, class S, class F>
+void
+LSVineyard<V,VE,S,F>::
+change_evaluator(Evaluator* eval)
+{
+ AssertMsg(evaluator_ != 0, "change_evaluator() assumes that existing evaluator is not null");
+
+ delete evaluator_;
+ evaluator_ = eval;
+ vineyard_.set_evaluator(evaluator_);
+}
+
+#include <boost/function.hpp>
+#include <boost/bind.hpp>
+#include <boost/lambda/lambda.hpp>
+namespace bl = boost::lambda;
+
+template<class V, class VE, class S, class F>
+bool
+LSVineyard<V,VE,S,F>::
+transpose_vertices(VertexIndex vi)
+{
+ Count(cVertexTransposition);
+ rLog(rlLSVineyard, "Transposing vertices (%d, %d)", vi->vertex(), boost::next(vi)->vertex());
+
+ DimensionFromIterator dim(pfmap_);
+ TranspositionVisitor visitor(*this);
+
+ iterator i = vi->simplex_index();
+ iterator i_prev = boost::prior(i);
+ iterator i_next = boost::next(vi)->simplex_index();
+ iterator i_next_prev = boost::prior(i_next); // transpositions are done in terms of the first index in the pair
+ iterator j = boost::next(i_next);
+
+ VertexIndex vi_next = boost::next(vi);
+ const Vertex& v = vi->vertex();
+
+ bool result = false; // has a switch in pairing occurred
+
+ // First, move the vertex --- this can be sped up if we devise special "vertex transpose" operation
+ rLog(rlLSVineyardDebug, "Starting to move the vertex");
+ while (i_next_prev != i_prev)
+ {
+ rLog(rlLSVineyardDebug, " Transposing %s %s", tostring(pfmap(i_next_prev)).c_str(),
+ tostring(pfmap(boost::next(i_next_prev))).c_str());
+ result |= persistence_.transpose(i_next_prev, dim, visitor);
+ i_next_prev = boost::prior(i_next);
+ }
+ rLog(rlLSVineyardDebug, "Done moving the vertex");
+
+ // Second, move the simplices attached to it
+ rLog(rlLSVineyardDebug, "Moving attached simplices");
+ while (j != persistence_.end() && j->attachment == vi_next)
+ {
+ rLog(rlLSVineyardDebug, " Considering %s", tostring(pfmap(j)).c_str());
+ if (pfmap(j).contains(v)) // j becomes attached to v and does not move
+ {
+ Count(cAttachment);
+ rLog(rlLSVineyardDebug, " Attachment changed for ", tostring(pfmap(j)).c_str());
+ set_attachment(j, vi);
+ ++j;
+ continue;
+ }
+
+ iterator j_prev = j; ++j;
+ while ((--j_prev)->attachment != vi_next) // i.e., until we have reached vi_next (and the simplices that follow it) again
+ {
+ rLog(rlLSVineyardDebug, " Moving: %s, %s",
+ tostring(pfmap(j_prev)).c_str(),
+ tostring(pfmap(boost::next(j_prev))).c_str());
+ AssertMsg(j_prev->attachment == vi, "Simplex preceding the one being moved must be attached to v");
+ result |= persistence_.transpose(j_prev, dim, visitor);
+ --j_prev;
+ }
+ }
+ rLog(rlLSVineyard, "Done moving attached simplices");
+ vertices_.splice(vi, vertices_, vi_next); // swap vi and vi_next
+
+ return result;
+}
+
+
+/* Evaluators */
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::StaticEvaluator: public Evaluator
+{
+ public:
+ StaticEvaluator(const LSVineyard& v, RealType time):
+ time_(time), vineyard_(v) {}
+
+ virtual RealType time() const { return time_; }
+ virtual RealType operator()(Index i) const { return vineyard_.simplex_value(vineyard_.pfmap(i)); }
+ virtual Dimension dimension(Index i) const { return vineyard_.pfmap(i).dimension(); }
+
+ private:
+ RealType time_;
+ const LSVineyard& vineyard_;
+};
+
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::KineticEvaluator: public Evaluator
+{
+ public:
+ KineticEvaluator(const LSVineyard& v, Simulator::Handle sp, ActivePointsTable::Handle apt, RealType time_offset):
+ vineyard_(v), sp_(sp), apt_(apt), time_offset_(time_offset) {}
+
+ virtual RealType time() const { return time_offset_ + CGAL::to_double(get_time()); }
+ virtual RealType operator()(Index i) const
+ {
+ rLog(rlLSVineyard, "%s (attached to %d): %s(%f) = %f", tostring(vineyard_.pfmap(i)).c_str(),
+ i->attachment->vertex(),
+ tostring(apt_->at(i->attachment->kinetic_key()).x()).c_str(),
+ get_time(),
+ CGAL::to_double(apt_->at(i->attachment->kinetic_key()).x()(get_time())));
+ return CGAL::to_double(apt_->at(i->attachment->kinetic_key()).x()(get_time()));
+ }
+ virtual Dimension dimension(Index i) const { return vineyard_.pfmap(i).dimension(); }
+
+ private:
+ Simulator::Time get_time() const { return sp_->current_time(); }
+
+ const LSVineyard& vineyard_;
+ Simulator::Handle sp_;
+ ActivePointsTable::Handle apt_;
+ RealType time_offset_;
+};
--- a/examples/grid/pdbdistance-vineyard.cpp Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/grid/pdbdistance-vineyard.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -2,7 +2,7 @@
#include "utilities/log.h"
#include "pdbdistance.h"
-#include "grid2Dvineyard.h"
+#include "lsvineyard.h"
#include <fstream>
#include <string>
@@ -51,10 +51,20 @@
// Compute initial filtration
int f = 0; int sf = 0;
std::ifstream in(frame_filename(infilename, f, sf++).c_str());
- Grid2DVineyard v(new PDBDistanceGrid(in, cas_only));
+ PDBDistanceGrid ginit(in, cas_only);
in.close();
- std::cout << "Filtration generated, size: " << v.filtration()->size() << std::endl;
- v.compute_pairing();
+
+ typedef LSVineyard<Grid2D::CoordinateIndex, Grid2D> Grid2DVineyard;
+
+ Grid2DVineyard::LSFiltration simplices;
+ ginit.complex_generator(make_push_back_functor(simplices));
+ Grid2DVineyard::VertexComparison vcmp(ginit);
+ Grid2DVineyard::SimplexComparison scmp(vcmp);
+ simplices.sort(scmp);
+ std::cout << "Complex generated, size: " << simplices.size() << std::endl;
+
+ Grid2DVineyard v(ginit.begin(), ginit.end(), simplices, ginit);
+ std::cout << "Filtration generated, size: " << v.filtration().size() << std::endl;
std::cout << "Pairing computed" << std::endl;
// Process frames computing the vineyard
@@ -63,13 +73,13 @@
std::string fn = frame_filename(infilename, f, sf++);
std::cout << "Processing " << fn << std::endl;
in.open(fn.c_str());
- v.compute_vineyard(new PDBDistanceGrid(in, cas_only));
+ v.compute_vineyard(PDBDistanceGrid(in, cas_only));
in.close();
if (sf == lastsubframe) { sf = 0; ++f; }
}
std::cout << "Vineyard computed" << std::endl;
- v.vineyard()->save_edges(outfilename);
+ v.vineyard().save_edges(outfilename);
#if 0
std::ofstream ofs(outfilename.c_str(), std::ios::binary);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/grid/pl-vineyard.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -0,0 +1,169 @@
+#include <utilities/log.h>
+
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+
+#include "lsvineyard.h"
+
+#include <boost/program_options.hpp>
+#include <boost/iterator/counting_iterator.hpp>
+#include <boost/function.hpp>
+#include <boost/bind.hpp>
+#include <boost/lambda/lambda.hpp>
+namespace bl = boost::lambda;
+
+
+typedef float VertexValue;
+typedef unsigned Vertex;
+typedef std::vector<VertexValue> VertexVector;
+struct SubscriptFunctor: public std::unary_function<Vertex, VertexValue>
+{
+ SubscriptFunctor(const VertexVector& v): vec(&v) {}
+ float operator()(Vertex i) const { return (*vec)[i]; }
+ SubscriptFunctor& operator=(const SubscriptFunctor& other) { vec = other.vec; return *this; }
+ const VertexVector* vec;
+};
+typedef SubscriptFunctor VertexEvaluator;
+typedef std::vector<VertexVector> VertexVectorVector;
+typedef LSVineyard<Vertex, VertexEvaluator> PLVineyard;
+typedef PLVineyard::Simplex Smplx; // gotta start using namespaces
+
+void program_options(int argc, char* argv[], std::string& complex_fn,
+ std::string& values_fn,
+ std::string& output_prefix,
+ bool& skip_infinite_vines,
+ bool& save_vines,
+ bool& explicit_events);
+void read_simplices(const std::string& complex_fn, PLVineyard::LSFiltration& simplices);
+void read_vertices(const std::string& vertex_fn, VertexVectorVector& vertices);
+
+
+int main(int argc, char** argv)
+{
+#ifdef LOGGING
+ rlog::RLogInit(argc, argv);
+#endif
+
+ std::string complex_fn, values_fn, output_prefix;
+ bool skip_infinite_vines = false, explicit_events = false, save_vines = false;
+ program_options(argc, argv, complex_fn, values_fn, output_prefix, skip_infinite_vines, save_vines, explicit_events);
+
+
+ // Read in the complex
+ PLVineyard::LSFiltration simplices;
+ read_simplices(complex_fn, simplices);
+ std::cout << "Complex read, size: " << simplices.size() << std::endl;
+
+ // Read in vertex values
+ VertexVectorVector vertices;
+ read_vertices(values_fn, vertices);
+
+ // Setup the vineyard
+ VertexEvaluator veval(vertices[0]);
+ PLVineyard::VertexComparison vcmp(veval);
+ PLVineyard::SimplexComparison scmp(vcmp);
+ simplices.sort(scmp);
+ PLVineyard v(boost::counting_iterator<Vertex>(0),
+ boost::counting_iterator<Vertex>(vertices[0].size()),
+ simplices, veval);
+ std::cout << "Pairing computed" << std::endl;
+
+ // Compute vineyard
+ for (size_t i = 1; i < vertices.size(); ++i)
+ {
+ veval = VertexEvaluator(vertices[i]);
+ v.compute_vineyard(veval, explicit_events);
+ std::cout << "Processed frame: " << i << std::endl;
+ }
+ std::cout << "Vineyard computed" << std::endl;
+
+ if (save_vines)
+ v.vineyard().save_vines(output_prefix, skip_infinite_vines);
+ else
+ v.vineyard().save_edges(output_prefix, skip_infinite_vines);
+}
+
+
+void read_simplices(const std::string& complex_fn, PLVineyard::LSFiltration& simplices)
+{
+ std::ifstream in(complex_fn.c_str());
+ std::string line;
+ while (std::getline(in, line))
+ {
+ std::istringstream strin(line);
+ simplices.push_back(Smplx(std::istream_iterator<Vertex>(strin), std::istream_iterator<Vertex>()));
+ }
+ std::cout << "Simplices read:" << std::endl;
+ std::copy(simplices.begin(), simplices.end(), std::ostream_iterator<Smplx>(std::cout, "\n"));
+}
+
+void read_vertices(const std::string& vertex_fn, VertexVectorVector& vertices)
+{
+ std::ifstream in(vertex_fn.c_str());
+ std::string line;
+ while (std::getline(in, line))
+ {
+ std::istringstream strin(line);
+ vertices.push_back(VertexVector(std::istream_iterator<VertexValue>(strin), std::istream_iterator<VertexValue>()));
+ }
+ std::cout << "Vertex values read:" << std::endl;
+ for (size_t i = 0; i < vertices.size(); ++i)
+ {
+ std::copy(vertices[i].begin(), vertices[i].end(), std::ostream_iterator<VertexValue>(std::cout, " "));
+ std::cout << std::endl;
+ }
+}
+
+void program_options(int argc, char* argv[], std::string& complex_fn,
+ std::string& values_fn,
+ std::string& output_prefix,
+ bool& skip_infinite_vines,
+ bool& save_vines,
+ bool& explicit_events)
+{
+ namespace po = boost::program_options;
+
+ po::options_description hidden("Hidden options");
+ hidden.add_options()
+ ("complex-file", po::value<std::string>(&complex_fn), "file listing the simplices of the complex")
+ ("values-file", po::value<std::string>(&values_fn), "file listing the values at the vertices")
+ ("output-prefix", po::value<std::string>(&output_prefix), "output prefix");
+
+ po::options_description visible("Allowed options", 100);
+ visible.add_options()
+ ("skip-infinite,s", po::bool_switch(&skip_infinite_vines), "skip infinite vines in the output")
+ ("explicit-events,e", po::bool_switch(&explicit_events), "process kinetic sort events one by one")
+ ("save-vines,v", po::bool_switch(&save_vines), "save vines instead of edges")
+ ("help,h", "produce help message");
+#if LOGGING
+ std::vector<std::string> log_channels;
+ visible.add_options()
+ ("log,l", po::value< std::vector<std::string> >(&log_channels), "log channels to turn on (info, debug, etc)");
+#endif
+
+ po::positional_options_description pos;
+ pos.add("complex-file", 1);
+ pos.add("values-file", 1);
+ pos.add("output-prefix", 1);
+
+ po::options_description all; all.add(visible).add(hidden);
+
+ po::variables_map vm;
+ po::store(po::command_line_parser(argc, argv).
+ options(all).positional(pos).run(), vm);
+ po::notify(vm);
+
+#if LOGGING
+ for (std::vector<std::string>::const_iterator cur = log_channels.begin(); cur != log_channels.end(); ++cur)
+ stderrLog.subscribeTo( RLOG_CHANNEL(cur->c_str()) );
+#endif
+
+ if (vm.count("help") || !vm.count("complex-file") || !vm.count("values-file") || !vm.count("output-prefix"))
+ {
+ std::cout << "Usage: " << argv[0] << " [options] complex-file values-file output-prefix" << std::endl;
+ std::cout << visible << std::endl;
+ std::abort();
+ }
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/grid/test-grid2D-vineyard.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -0,0 +1,54 @@
+#include <utilities/log.h>
+
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+
+#include "grid2D.h"
+#include "lsvineyard.h"
+
+int main(int argc, char** argv)
+{
+#ifdef LOGGING
+ rlog::RLogInit(argc, argv);
+ // stdoutLog.subscribeTo(RLOG_CHANNEL("topology"));
+ // stdoutLog.subscribeTo(RLOG_CHANNEL("topology/persistence/transpositions"));
+ // stdoutLog.subscribeTo(RLOG_CHANNEL("topology/vineyard"));
+ stdoutLog.subscribeTo(RLOG_CHANNEL("lsvineyard"));
+#endif
+
+ Grid2D g0(2, 2), g1(2, 2);
+ g0(0,0) = 1; g0(0,1) = 2; g0(1,0) = 3; g0(1,1) = 0;
+ g1(0,0) = 4; g1(0,1) = 2; g1(1,0) = 3; g1(1,1) = 5;
+
+ // Generate the complex, initialize the vineyard (which also computes the pairing)
+ typedef LSVineyard<Grid2D::CoordinateIndex, Grid2D> Grid2DVineyard;
+
+ Grid2DVineyard::LSFiltration simplices;
+ g0.complex_generator(make_push_back_functor(simplices));
+ Grid2DVineyard::VertexComparison vcmp(g0);
+ Grid2DVineyard::SimplexComparison scmp(vcmp);
+ simplices.sort(scmp);
+ std::cout << "Complex generated, size: " << simplices.size() << std::endl;
+ std::copy(simplices.begin(), simplices.end(), std::ostream_iterator<Grid2D::Smplx>(std::cout, "\n"));
+
+ Grid2DVineyard v(g0.begin(), g0.end(), simplices, g0);
+ std::cout << "Filtration generated, size: " << v.filtration().size() << std::endl;
+ std::cout << "Pairing computed" << std::endl;
+
+ // Simplex order before
+ std::cout << "Simplex order:" << std::endl;
+ for (Grid2DVineyard::LSFiltration::Index cur = v.filtration().begin(); cur != v.filtration().end(); ++cur)
+ std::cout << " " << v.filtration().simplex(cur) << std::endl;
+
+ // Compute vineyard
+ v.compute_vineyard(g1, true);
+ std::cout << "Vineyard computed" << std::endl;
+
+ // Simplex order after
+ std::cout << "Simplex order:" << std::endl;
+ for (Grid2DVineyard::LSFiltration::Index cur = v.filtration().begin(); cur != v.filtration().end(); ++cur)
+ std::cout << " " << v.filtration().simplex(cur) << std::endl;
+
+ v.vineyard().save_edges("test-vineyard");
+}
--- a/examples/lsfiltration.py Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/lsfiltration.py Sun Dec 20 08:27:01 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/poincare/poincare.cpp Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/poincare/poincare.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -10,8 +10,7 @@
#include <boost/program_options.hpp>
typedef Simplex<unsigned, unsigned> Smplx;
-typedef std::vector<Smplx> Complex;
-typedef Filtration<Complex> Fltr;
+typedef Filtration<Smplx> Fltr;
typedef StaticPersistence<> Persistence;
typedef PersistenceDiagram<> PDgm;
@@ -48,8 +47,7 @@
}
- Complex c;
-
+ Fltr f;
std::ifstream in(infilename.c_str());
unsigned int i = 0;
std::string s;
@@ -66,21 +64,19 @@
linestream >> vertex;
}
std::cout << simplex << std::endl;
- c.push_back(simplex);
+ f.push_back(simplex);
std::getline(in, s);
}
- std::sort(c.begin(), c.end(), Smplx::VertexComparison());
- Fltr f(c.begin(), c.end(), Smplx::DataComparison());
+ f.sort(Smplx::DataComparison());
Persistence pers(f);
pers.pair_simplices();
+ Persistence::SimplexMap<Fltr> m = pers.make_simplex_map(f);
std::map<Dimension, PDgm> dgms;
init_diagrams(dgms, pers.begin(), pers.end(),
- evaluate_through_map(OffsetMap<Persistence::OrderIndex, Fltr::Index>(pers.begin(), f.begin()),
- evaluate_through_filtration(f, Smplx::DataEvaluator())),
- evaluate_through_map(OffsetMap<Persistence::OrderIndex, Fltr::Index>(pers.begin(), f.begin()),
- evaluate_through_filtration(f, Smplx::DimensionExtractor())));
+ evaluate_through_map(m, Smplx::DataEvaluator()),
+ evaluate_through_map(m, Smplx::DimensionExtractor()));
std::cout << 0 << std::endl << dgms[0] << std::endl;
std::cout << 1 << std::endl << dgms[1] << std::endl;
--- a/examples/rips/rips-pairwise.cpp Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/rips/rips-pairwise.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -21,10 +21,9 @@
typedef Rips<PairDistances> Generator;
typedef Generator::Simplex Smplx;
-typedef std::vector<Smplx> SimplexVector;
-typedef Filtration<SimplexVector, unsigned> Fltr;
-typedef StaticPersistence<> Persistence;
-// typedef DynamicPersistenceChains<> Persistence;
+typedef Filtration<Smplx> Fltr;
+// typedef StaticPersistence<> Persistence;
+typedef DynamicPersistenceChains<> Persistence;
typedef PersistenceDiagram<> PDgm;
void program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, std::string& diagram_name);
@@ -45,44 +44,45 @@
PairDistances distances(points);
Generator rips(distances);
Generator::Evaluator size(distances);
- SimplexVector complex;
+ Fltr f;
// Generate 2-skeleton of the Rips complex for epsilon = 50
- rips.generate(skeleton, max_distance, make_push_back_functor(complex));
- std::sort(complex.begin(), complex.end(), Smplx::VertexComparison()); // unnecessary
- std::cout << "Generated complex of size: " << complex.size() << std::endl;
+ rips.generate(skeleton, max_distance, make_push_back_functor(f));
+ std::cout << "# Generated complex of size: " << f.size() << std::endl;
// Generate filtration with respect to distance and compute its persistence
- Fltr f(complex.begin(), complex.end(), Generator::Comparison(distances));
+ f.sort(Generator::Comparison(distances));
Timer persistence_timer; persistence_timer.start();
Persistence p(f);
p.pair_simplices();
persistence_timer.stop();
+#if 1
// Output cycles
- for (Persistence::OrderIndex cur = p.begin(); cur != p.end(); ++cur)
+ Persistence::SimplexMap<Fltr> m = p.make_simplex_map(f);
+ for (Persistence::iterator cur = p.begin(); cur != p.end(); ++cur)
{
- Persistence::Cycle& cycle = cur->cycle;
+ const Persistence::Cycle& cycle = cur->cycle;
if (!cur->sign()) // only negative simplices have non-empty cycles
{
Persistence::OrderIndex birth = cur->pair; // the cycle that cur killed was born when we added birth (another simplex)
- const Smplx& b = f.simplex(f.begin() + (birth - p.begin())); // eventually this will be encapsulated behind an interface
- const Smplx& d = f.simplex(f.begin() + (cur - p.begin()));
+ const Smplx& b = m[birth];
+ const Smplx& d = m[cur];
// if (b.dimension() != 1) continue;
// std::cout << "Pair: (" << size(b) << ", " << size(d) << ")" << std::endl;
if (b.dimension() >= skeleton) continue;
diagram_out << b.dimension() << " " << size(b) << " " << size(d) << std::endl;
- } else if (cur->pair == cur) // positive could be unpaired
+ } else if (cur->unpaired()) // positive could be unpaired
{
- const Smplx& b = f.simplex(f.begin() + (cur - p.begin()));
+ const Smplx& b = m[cur];
// if (b.dimension() != 1) continue;
// std::cout << "Unpaired birth: " << size(b) << std::endl;
- // cycle = cur->chain;
+ // cycle = cur->chain; // TODO
if (b.dimension() >= skeleton) continue;
diagram_out << b.dimension() << " " << size(b) << " inf" << std::endl;
}
@@ -91,13 +91,14 @@
// for (Persistence::Cycle::const_iterator si = cycle.begin();
// si != cycle.end(); ++si)
// {
- // const Smplx& s = f.simplex(f.begin() + (*si - p.begin()));
+ // const Smplx& s = m[*si];
// //std::cout << s.dimension() << std::endl;
// const Smplx::VertexContainer& vertices = s.vertices(); // std::vector<Vertex> where Vertex = Distances::IndexType
// AssertMsg(vertices.size() == s.dimension() + 1, "dimension of a simplex is one less than the number of its vertices");
// std::cout << vertices[0] << " " << vertices[1] << std::endl;
// }
}
+#endif
persistence_timer.check("# Persistence timer");
}
--- a/examples/rips/rips-pairwise.py Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/rips/rips-pairwise.py Sun Dec 20 08:27:01 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-weighted.cpp Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/rips/rips-weighted.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -27,8 +27,7 @@
typedef WeightedRips<PairDistances> Generator;
typedef Generator::Simplex Smplx;
-typedef std::vector<Smplx> SimplexVector;
-typedef Filtration<SimplexVector, unsigned> Fltr;
+typedef Filtration<Smplx> Fltr;
typedef StaticPersistence<> Persistence;
//typedef DynamicPersistenceChains<> Persistence;
typedef PersistenceDiagram<> PDgm;
@@ -56,39 +55,38 @@
Generator rips(distances);
Generator::Evaluator size(distances);
Generator::Comparison cmp (distances);
- SimplexVector complex;
+ Fltr complex;
// Generate skeleton of the weighted Rips complex for epsilon = 50
rips.generate(skeleton, max_distance, make_push_back_functor(complex));
-
- std::sort(complex.begin(), complex.end(), Smplx::VertexComparison()); // unnecessary
std::cout << "# Generated complex of size: " << complex.size() << std::endl;
// Generate filtration with respect to distance and compute its persistence
- Fltr f(complex.begin(), complex.end(), cmp);
+ complex.sort(cmp);
Timer persistence_timer; persistence_timer.start();
- Persistence p(f);
+ Persistence p(complex);
p.pair_simplices();
persistence_timer.stop();
// Output cycles
- for (Persistence::OrderIndex cur = p.begin(); cur != p.end(); ++cur)
+ Persistence::SimplexMap<Fltr> m = p.make_simplex_map(complex);
+ for (Persistence::iterator cur = p.begin(); cur != p.end(); ++cur)
{
- Persistence::Cycle& cycle = cur->cycle;
+ const Persistence::Cycle& cycle = cur->cycle;
if (!cur->sign()) // only negative simplices have non-empty cycles
{
Persistence::OrderIndex birth = cur->pair; // the cycle that cur killed was born when we added birth (another simplex)
- const Smplx& b = f.simplex(f.begin() + (birth - p.begin())); // eventually this will be encapsulated behind an interface
- const Smplx& d = f.simplex(f.begin() + (cur - p.begin()));
+ const Smplx& b = m[birth];
+ const Smplx& d = m[cur];
if (b.dimension() >= skeleton) continue;
std::cout << b.dimension() << " " << size(b) << " " << size(d) << std::endl;
- } else if (cur->pair == cur) // positive could be unpaired
+ } else if (cur->unpaired()) // positive could be unpaired
{
- const Smplx& b = f.simplex(f.begin() + (cur - p.begin()));
+ const Smplx& b = m[cur];
if (b.dimension() >= skeleton) continue;
std::cout << b.dimension() << " " << size(b) << " inf" << std::endl;
--- a/examples/rips/rips.cpp Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/rips/rips.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -25,10 +25,9 @@
//typedef Rips<ExplicitDistances<Distances> > Generator;
typedef Rips<Distances> Generator;
typedef Generator::Simplex Smplx;
-typedef std::vector<Smplx> SimplexVector;
-typedef Filtration<SimplexVector, unsigned> Fltr;
-//typedef StaticPersistence<> Persistence;
-typedef DynamicPersistenceChains<> Persistence;
+typedef Filtration<Smplx> Fltr;
+typedef StaticPersistence<> Persistence;
+// typedef DynamicPersistenceChains<> Persistence;
typedef PersistenceDiagram<> PDgm;
@@ -50,26 +49,25 @@
Generator rips(distances);
Generator::Evaluator size(distances);
- SimplexVector complex;
+ Fltr f;
// Generate 2-skeleton of the Rips complex for epsilon = 50
- rips.generate(2, 10, make_push_back_functor(complex));
- std::sort(complex.begin(), complex.end(), Smplx::VertexComparison()); // unnecessary
- rInfo("Generated complex of size: %d", complex.size());
+ rips.generate(2, 10, make_push_back_functor(f));
+ rInfo("Generated complex of size: %d", f.size());
// Generate filtration with respect to distance and compute its persistence
- Fltr f(complex.begin(), complex.end(), Generator::Comparison(distances));
+ f.sort(Generator::Comparison(distances));
Persistence p(f);
p.pair_simplices();
rInfo("Simplices paired");
+ Persistence::SimplexMap<Fltr> m = p.make_simplex_map(f);
+
// Record the persistence intervals in the persistence diagrams
std::map<Dimension, PDgm> dgms;
init_diagrams(dgms, p.begin(), p.end(),
- evaluate_through_map(make_offset_map(p.begin(), f.begin()),
- evaluate_through_filtration(f, size)),
- evaluate_through_map(make_offset_map(p.begin(), f.begin()),
- evaluate_through_filtration(f, Smplx::DimensionExtractor())));
+ evaluate_through_map(m, size),
+ evaluate_through_map(m, Smplx::DimensionExtractor()));
// Serialize the diagrams to a file
std::ofstream ofs("rips-diagrams");
@@ -77,33 +75,32 @@
oa << dgms;
// Output cycles
- for (Persistence::OrderIndex cur = p.begin(); cur != p.end(); ++cur)
+ for (Persistence::iterator cur = p.begin(); cur != p.end(); ++cur)
{
- Persistence::Cycle& cycle = cur->cycle;
+ const Persistence::Cycle& cycle = cur->cycle;
if (!cur->sign()) // only negative simplices have non-empty cycles
{
Persistence::OrderIndex birth = cur->pair; // the cycle that cur killed was born when we added birth (another simplex)
- const Smplx& b = f.simplex(f.begin() + (birth - p.begin())); // eventually this will be encapsulated behind an interface
- const Smplx& d = f.simplex(f.begin() + (cur - p.begin()));
+ const Smplx& b = m[birth];
+ const Smplx& d = m[cur];
if (b.dimension() != 1) continue;
std::cout << "Pair: (" << size(b) << ", " << size(d) << ")" << std::endl;
- } else if (cur->pair == cur) // positive could be unpaired
+ } else if (cur->unpaired()) // positive could be unpaired
{
- const Smplx& b = f.simplex(f.begin() + (cur - p.begin()));
+ const Smplx& b = m[cur];
if (b.dimension() != 1) continue;
std::cout << "Unpaired birth: " << size(b) << std::endl;
- cycle = cur->chain;
+ // cycle = cur->chain;
}
// Iterate over the cycle
- for (Persistence::Cycle::const_iterator si = cycle.begin();
- si != cycle.end(); ++si)
+ for (Persistence::Cycle::const_iterator si = cycle.begin(); si != cycle.end(); ++si)
{
- const Smplx& s = f.simplex(f.begin() + (*si - p.begin()));
+ const Smplx& s = m[*si];
//std::cout << s.dimension() << std::endl;
const Smplx::VertexContainer& vertices = s.vertices(); // std::vector<Vertex> where Vertex = Distances::IndexType
AssertMsg(vertices.size() == s.dimension() + 1, "dimension of a simplex is one less than the number of its vertices");
--- a/examples/rips/rips.py Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/rips/rips.py Sun Dec 20 08:27:01 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-zigzag.py Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/triangle/triangle-zigzag.py Sun Dec 20 08:27:01 2009 -0800
@@ -27,7 +27,7 @@
b += 1
# Remove all the simplices
-for s in sorted(complex.keys(), reverse = True):
+for s in sorted(complex.keys(), data_cmp, reverse = True):
print "%d: Removing %s" % (b, s)
d = zz.remove(complex[s], b)
del complex[s]
--- a/examples/triangle/triangle.cpp Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/triangle/triangle.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -2,7 +2,7 @@
#include "topology/simplex.h"
#include "topology/filtration.h"
-//#include "topology/static-persistence.h"
+#include "topology/static-persistence.h"
#include "topology/dynamic-persistence.h"
#include "topology/persistence-diagram.h"
#include <utilities/indirect.h>
@@ -19,15 +19,34 @@
#include <boost/serialization/vector.hpp>
#endif
-typedef unsigned Vertex;
-typedef Simplex<Vertex, double> Smplx;
-typedef std::vector<Smplx> Complex;
-typedef Filtration<Complex, unsigned> Fltr;
-//typedef StaticPersistence<> Persistence;
-typedef DynamicPersistenceTrails<> Persistence;
-typedef PersistenceDiagram<> PDgm;
+typedef unsigned Vertex;
+typedef Simplex<Vertex, double> Smplx;
+typedef Filtration<Smplx> Fltr;
+// typedef StaticPersistence<> Persistence;
+typedef DynamicPersistenceTrails<> Persistence;
+typedef PersistenceDiagram<> PDgm;
+typedef OffsetBeginMap<Persistence, Fltr,
+ Persistence::iterator,
+ Fltr::Index> PersistenceFiltrationMap;
+typedef OffsetBeginMap<Fltr, Persistence,
+ Fltr::Index,
+ Persistence::iterator> FiltrationPersistenceMap;
-void fillTriangleSimplices(Complex& c)
+// Transposes elements of the filtration together with the
+struct FiltrationTranspositionVisitor: public Persistence::TranspositionVisitor
+{
+ typedef Persistence::iterator iterator;
+
+ FiltrationTranspositionVisitor(const Persistence& p,
+ Fltr& f):
+ p_(p), f_(f) {}
+ void transpose(iterator i) { f_.transpose(f_.begin() + (i - p_.begin())); }
+
+ const Persistence& p_;
+ Fltr& f_;
+};
+
+void fillTriangleSimplices(Fltr& c)
{
typedef std::vector<Vertex> VertexVector;
VertexVector vertices(4);
@@ -43,8 +62,6 @@
c.push_back(Smplx(bg + 1, bg + 3, 2.9)); // BC
c.push_back(Smplx(bg + 2, end, 3.5)); // CA
c.push_back(Smplx(bg, bg + 3, 5)); // ABC
-
- std::sort(c.begin(), c.end(), Smplx::VertexComparison());
}
int main(int argc, char** argv)
@@ -57,28 +74,28 @@
//stdoutLog.subscribeTo(RLOG_CHANNEL("topology/vineyard"));
#endif
- Complex c;
- fillTriangleSimplices(c);
+ Fltr f;
+ fillTriangleSimplices(f);
std::cout << "Simplices filled" << std::endl;
- for (Complex::const_iterator cur = c.begin(); cur != c.end(); ++cur)
+ for (Fltr::Index cur = f.begin(); cur != f.end(); ++cur)
std::cout << " " << *cur << std::endl;
-#if 1 // testing serialization of Complex (really Simplex)
+#if 1 // testing serialization of the Filtration (really Simplex)
{
std::ofstream ofs("complex");
boost::archive::text_oarchive oa(ofs);
- oa << c;
- c.clear();
+ oa << f;
+ f.clear();
}
{
std::ifstream ifs("complex");
boost::archive::text_iarchive ia(ifs);
- ia >> c;
+ ia >> f;
}
#endif
- Fltr f(c.begin(), c.end(), Smplx::DataComparison());
+ f.sort(Smplx::DataComparison());
std::cout << "Filtration initialized" << std::endl;
std::cout << f << std::endl;
@@ -88,36 +105,36 @@
p.pair_simplices();
std::cout << "Simplices paired" << std::endl;
+ Persistence::SimplexMap<Fltr> m = p.make_simplex_map(f);
std::map<Dimension, PDgm> dgms;
init_diagrams(dgms, p.begin(), p.end(),
- evaluate_through_map(make_offset_map(p.begin(), f.begin()),
- evaluate_through_filtration(f, Smplx::DataEvaluator())),
- evaluate_through_map(make_offset_map(p.begin(), f.begin()),
- evaluate_through_filtration(f, Smplx::DimensionExtractor())));
+ evaluate_through_map(m, Smplx::DataEvaluator()),
+ evaluate_through_map(m, Smplx::DimensionExtractor()));
std::cout << 0 << std::endl << dgms[0] << std::endl;
std::cout << 1 << std::endl << dgms[1] << std::endl;
+ PersistenceFiltrationMap pfmap(p, f);
+ DimensionFunctor<PersistenceFiltrationMap, Fltr> dim(pfmap, f);
+
// Transpositions
- p.transpose(p.begin()); // transposition case 1.2 special
-
-#if 0
+ FiltrationPersistenceMap fpmap(f, p);
+ FiltrationTranspositionVisitor visitor(p, f);
Smplx A; A.add(0);
std::cout << A << std::endl;
- std::cout << *tf.get_index(A) << std::endl;
- std::cout << "Transposing A: " << tf.transpose(tf.get_index(A)) << std::endl;
- std::cout << tf;
+ std::cout << "Transposing A: " << p.transpose(fpmap[f.find(A)], dim, visitor) << std::endl; // 1.2 unpaired
Smplx BC; BC.add(1); BC.add(2);
Smplx AB; AB.add(0); AB.add(1);
std::cout << BC << std::endl;
- std::cout << *tf.get_index(BC) << std::endl;
- tf.transpose(tf.get_index(BC));
- std::cout << tf;
+ std::cout << p.transpose(fpmap[f.find(BC)], dim, visitor) << std::endl; // 3.1
+ // p.transpose(fpmap[f.find(BC)], dim, visitor);
std::cout << AB << std::endl;
- std::cout << *tf.get_index(AB) << std::endl;
- tf.transpose(tf.get_index(AB));
- std::cout << tf;
-#endif
+ std::cout << p.transpose(fpmap[f.find(AB)], dim, visitor) << std::endl; // 2.1
+ // p.transpose(fpmap[f.find(AB)], dim, visitor);
+
+ std::cout << p.transpose(p.begin(), dim, visitor) << std::endl; // transposition case 1.2 special
+ std::cout << p.transpose(boost::next(p.begin()), dim, visitor) << std::endl;
+ std::cout << p.transpose(boost::next(p.begin(),3), dim, visitor) << std::endl;
}
--- a/examples/triangle/triangle.py Thu Dec 03 13:21:23 2009 -0800
+++ b/examples/triangle/triangle.py Sun Dec 20 08:27:01 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()])
--- a/include/topology/chain.h Thu Dec 03 13:21:23 2009 -0800
+++ b/include/topology/chain.h Sun Dec 20 08:27:01 2009 -0800
@@ -79,6 +79,7 @@
/// @{
void remove(iterator i) { ChainRepresentation::erase(i); Size::operator--(); }
void remove(const_reference x) { remove(std::find(begin(), end(), x)); }
+ bool remove_if_contains(const_reference x);
template<class ConsistencyComparison>
void append(const_reference x, const ConsistencyComparison& cmp);
@@ -95,7 +96,7 @@
const_reference top(const OrderComparison& cmp) const; ///< First element in cmp order
boost::optional<const_iterator> contains(const_reference x) const; ///< tests whether chain contains x
- boost::optional<iterator> contains(reference x); ///< tests whether chain contains x
+ boost::optional<iterator> contains(const_reference x); ///< tests whether chain contains x
/// @}
/// \name Debugging
--- a/include/topology/chain.hpp Thu Dec 03 13:21:23 2009 -0800
+++ b/include/topology/chain.hpp Sun Dec 20 08:27:01 2009 -0800
@@ -101,19 +101,33 @@
contains(const_reference x) const
{
const_iterator res = std::find(begin(), end(), x);
- return make_optional(res != end(), res);
+ return boost::make_optional(res != end(), res);
}
template<class C>
boost::optional<typename ChainWrapper<C>::iterator>
ChainWrapper<C>::
-contains(reference x)
+contains(const_reference x)
{
iterator res = std::find(begin(), end(), x);
return boost::make_optional(res != end(), res);
}
template<class C>
+bool
+ChainWrapper<C>::
+remove_if_contains(const_reference x)
+{
+ boost::optional<iterator> i = contains(x);
+ if (i)
+ {
+ remove(*i);
+ return true;
+ } else
+ return false;
+}
+
+template<class C>
template<class OutputMap>
std::string
ChainWrapper<C>::
--- a/include/topology/complex-traits.h Thu Dec 03 13:21:23 2009 -0800
+++ b/include/topology/complex-traits.h Sun Dec 20 08:27:01 2009 -0800
@@ -43,6 +43,9 @@
static Index begin(const Complex& complex) { return complex.begin(); }
static Index end(const Complex& complex) { return complex.end(); }
+
+ static const Simplex&
+ simplex(Index i) { return *i; }
};
@@ -66,6 +69,9 @@
static Index begin(const Complex& complex) { return complex.begin(); }
static Index end(const Complex& complex) { return complex.end(); }
+
+ static const Simplex&
+ simplex(Index i) { return *i; }
};
#endif // __COMPLEX_TRAITS_H__
--- a/include/topology/dynamic-persistence.h Thu Dec 03 13:21:23 2009 -0800
+++ b/include/topology/dynamic-persistence.h Sun Dec 20 08:27:01 2009 -0800
@@ -4,42 +4,39 @@
#include "static-persistence.h"
#include <utilities/types.h>
+#include <boost/bind.hpp>
+#include <boost/lambda/lambda.hpp>
+namespace bl = boost::lambda;
+
#ifdef COUNTERS
static Counter* cTrailLength = GetCounter("persistence/pair/traillength"); // the size of matrix U in RU decomposition
static Counter* cChainLength = GetCounter("persistence/pair/chainlength"); // the size of matrix V in R=DV decomposition
#endif // COUNTERS
-template<class Data_, class ChainTraits_, class ContainerTraits_, class ConsistencyIndex_>
-struct TrailData: public PairCycleData<Data_, ChainTraits_, ContainerTraits_,
- TrailData<Data_, ChainTraits_, ContainerTraits_, ConsistencyIndex_> >
+template<class Data_, class ChainTraits_>
+struct TrailData: public PairCycleData<Data_, ChainTraits_, TrailData<Data_, ChainTraits_> >
{
typedef Data_ Data;
- typedef ConsistencyIndex_ ConsistencyIndex;
-
- typedef PairCycleData<Data_, ChainTraits_, ContainerTraits_, TrailData> Parent;
- typedef TrailData<Data_, ChainTraits_, ContainerTraits_, ConsistencyIndex_> Self;
- // typedef typename ContainerTraits_::template rebind<Self>::other ContainerTraits;
- // typedef typename ContainerTraits::Index Index;
- // typedef typename ChainTraits_::template rebind<Index>::other ChainTraits;
- // typedef typename ChainTraits::Chain Chain;
+ typedef PairCycleData<Data_, ChainTraits_, TrailData> Parent;
+ typedef TrailData<Data_, ChainTraits_> Self;
+
typedef typename Parent::Index Index;
+ typedef typename Parent::Cycle Cycle;
typedef typename Parent::Chain Chain;
typedef Chain Trail;
-
- template<class Comparison>
- struct ConsistencyComparison: public std::binary_function<const Index&, const Index&, bool>
- {
- ConsistencyComparison(const Comparison& cmp = Comparison()):
- cmp_(cmp) {}
+
+ // Modifiers
+ template<class Cmp>
+ void trail_append(Index i, const Cmp& cmp) { trail.append(i, cmp); }
+ template<class Cmp>
+ void trail_add(const Trail& t, const Cmp& cmp) { trail.add(t, cmp); }
- bool operator()(const Index& a, const Index& b) const { return cmp_(a->consistency, b->consistency); }
+ template<class Cmp>
+ void cycle_add(const Cycle& z, const Cmp& cmp) { cycle.add(z, cmp); }
- Comparison cmp_;
- };
-
- Trail trail;
- ConsistencyIndex consistency;
+ using Parent::cycle;
+ Trail trail;
};
/**
@@ -59,28 +56,31 @@
// position, or one could provide consistency that is references into the complex
template<class Data_ = Empty<>,
class ChainTraits_ = VectorChains<>,
- class Comparison_ = GreaterComparison<>,
- class ContainerTraits_ = VectorContainer<>,
- class ConsistencyIndex_ = size_t,
- class ConsistencyComparison_ = std::less<ConsistencyIndex_>,
- class Element_ = TrailData<Data_, ChainTraits_, ContainerTraits_, ConsistencyIndex_> >
+ class ContainerTraits_ = OrderConsistencyContainer<>,
+ class Element_ = TrailData<Data_, ChainTraits_>,
+ class Comparison_ = ElementComparison<typename ContainerTraits_::template rebind<Element_>::other::Container,
+ std::greater<typename ContainerTraits_::template rebind<Element_>::other::Container::iterator> >,
+ class ConsistencyComparison_ = ElementComparison<typename ContainerTraits_::template rebind<Element_>::other::ConsistentContainer,
+ std::greater<typename ContainerTraits_::template rebind<Element_>::other::ConsistentContainer::iterator> >
+ >
class DynamicPersistenceTrails:
- public StaticPersistence<Data_, ChainTraits_, Comparison_, ContainerTraits_, Element_>
+ public StaticPersistence<Data_, ChainTraits_, ContainerTraits_, Element_, Comparison_>
{
public:
typedef Data_ Data;
typedef Element_ Element;
- typedef StaticPersistence<Data_, ChainTraits_, Comparison_,
- ContainerTraits_, Element_> Parent;
+ typedef StaticPersistence<Data_, ChainTraits_,
+ ContainerTraits_, Element_, Comparison_> Parent;
typedef typename Parent::ContainerTraits Traits;
typedef typename Parent::Order Order;
typedef typename Parent::OrderComparison OrderComparison;
typedef typename Parent::OrderIndex OrderIndex;
- typedef ConsistencyIndex_ ConsistencyIndex;
- typedef ThreeOutcomeCompare<
- typename Element::
- template ConsistencyComparison<ConsistencyComparison_> > ConsistencyComparison;
+ typedef ConsistencyComparison_ ConsistencyComparison;
+ typedef typename Parent::iterator iterator;
+
+ typedef typename Element::Trail Trail;
+ typedef typename Element::Cycle Cycle;
/**
* Constructor: DynamicPersistenceTrails()
@@ -89,22 +89,23 @@
* Template parameters:
* Filtration - filtration of the complex whose persistence we are computing
*/
- template<class Filtration> DynamicPersistenceTrails(const Filtration& f,
- const OrderComparison& ocmp = OrderComparison(),
- const ConsistencyComparison& ccmp = ConsistencyComparison());
+ template<class Filtration> DynamicPersistenceTrails(const Filtration& f);
void pair_simplices();
// Function: transpose(i)
// Tranpose i and the next element.
// Returns: true iff the pairing switched.
- bool transpose(OrderIndex i) { return transpose(i, TranspositionVisitor()); }
+ template<class DimensionFunctor, class Visitor>
+ bool transpose(iterator i, const DimensionFunctor& dimension, Visitor visitor = Visitor());
- template<class Visitor>
- bool transpose(OrderIndex i, const Visitor& visitor = Visitor());
+ template<class DimensionFunctor>
+ bool transpose(iterator i, const DimensionFunctor& dimension) { return transpose(i,dimension,TranspositionVisitor()); }
using Parent::begin;
using Parent::end;
+ using Parent::iterator_to;
+ using Parent::index;
using Parent::size;
// Struct: TranspositionVisitor
@@ -116,66 +117,72 @@
// This function is called before transposition is processed
// (at the very beginning of <transpose(i, visitor)>). It is meant to update any structures
// that may need to be updated, but perhaps it has other uses as well.
- void transpose(OrderIndex i) const {}
+ void transpose(iterator i) const {}
// Function: switched(i, type)
// This function is called after the transposition if the switch in pairing has occured.
// `i` is the index of the preceding simplex after the transposition.
// `type` indicates the <SwitchType>.
- void switched(OrderIndex i, SwitchType type) const {}
+ void switched(iterator i, SwitchType type) const {}
};
protected:
using Parent::order;
+ using Parent::set_pair;
+ using Parent::swap_cycle;
+
+ bool trail_remove_if_contains
+ (iterator i, OrderIndex j) { TrailRemover rm(j); order().modify(i, rm); return rm.result; }
+ void cycle_add(iterator i, const Cycle& z) { order().modify(i, boost::bind(&Element::template cycle_add<ConsistencyComparison>, bl::_1, boost::ref(z), ccmp_)); } // i->cycle_add(z, ccmp_)
+ void trail_add(iterator i, const Trail& t) { order().modify(i, boost::bind(&Element::template trail_add<ConsistencyComparison>, bl::_1, boost::ref(t), ccmp_)); } // i->trail_add(t, ccmp_)
private:
- void swap(OrderIndex i, OrderIndex j);
- void pairing_switch(OrderIndex i, OrderIndex j);
+ void swap(iterator i, iterator j);
+ void pairing_switch(iterator i, iterator j);
struct PairingTrailsVisitor: public Parent::PairVisitor
{
- // TODO: this is specialized for std::vector
- PairingTrailsVisitor(OrderIndex bg, ConsistencyComparison ccmp, unsigned size):
- Parent::PairVisitor(size), bg_(bg), ccmp_(ccmp) {}
+ PairingTrailsVisitor(Order& order, ConsistencyComparison ccmp, unsigned size):
+ Parent::PairVisitor(size), order_(order), ccmp_(ccmp) {}
- void init(OrderIndex i) const { i->consistency = i - bg_; i->trail.append(i, ccmp_); Count(cTrailLength); }
- void update(OrderIndex j, OrderIndex i) const { i->pair->trail.append(j, ccmp_); Count(cTrailLength); }
- void finished(OrderIndex i) const { Parent::PairVisitor::finished(i); }
+ void init(iterator i) const { order_.modify(i, boost::bind(&Element::template trail_append<ConsistencyComparison>, bl::_1, &*i, ccmp_)); Count(cTrailLength); } // i->trail_append(&*i, ccmp)
+ void update(iterator j, iterator i) const { order_.modify(order_.iterator_to(*(i->pair)), boost::bind(&Element::template trail_append<ConsistencyComparison>, bl::_1, &*j, ccmp_)); Count(cTrailLength); } // i->pair->trail_append(&*j, ccmp)
+ void finished(iterator i) const { Parent::PairVisitor::finished(i); }
- OrderIndex bg_;
+ Order& order_;
ConsistencyComparison ccmp_;
};
+ struct TrailRemover;
+
ConsistencyComparison ccmp_;
};
-template<class Data_, class ChainTraits_, class ContainerTraits_, class ConsistencyIndex_>
-struct ChainData: public PairCycleData<Data_, ChainTraits_, ContainerTraits_,
- ChainData<Data_, ChainTraits_, ContainerTraits_, ConsistencyIndex_> >
+/* Chains */
+template<class Data_, class ChainTraits_>
+struct ChainData: public PairCycleData<Data_, ChainTraits_, ChainData<Data_, ChainTraits_> >
{
typedef Data_ Data;
- typedef ConsistencyIndex_ ConsistencyIndex;
- typedef PairCycleData<Data_, ChainTraits_, ContainerTraits_, ChainData> Parent;
- typedef ChainData<Data_, ChainTraits_, ContainerTraits_, ConsistencyIndex_> Self;
+ typedef PairCycleData<Data_, ChainTraits_, ChainData> Parent;
+ typedef ChainData<Data_, ChainTraits_> Self;
typedef typename Parent::Index Index;
+ typedef typename Parent::Cycle Cycle;
typedef typename Parent::Chain Chain;
typedef Chain Trail;
- template<class Comparison>
- struct ConsistencyComparison: public std::binary_function<const Index&, const Index&, bool>
- {
- ConsistencyComparison(const Comparison& cmp = Comparison()):
- cmp_(cmp) {}
+ // Modifiers
+ template<class Cmp>
+ void chain_append(Index i, const Cmp& cmp) { chain.append(i, cmp); }
+ template<class Cmp>
+ void chain_add(const Chain& c, const Cmp& cmp) { chain.add(c, cmp); }
- bool operator()(const Index& a, const Index& b) const { return cmp_(a->consistency, b->consistency); }
+ template<class Cmp>
+ void cycle_add(const Cycle& z, const Cmp& cmp) { cycle.add(z, cmp); }
- Comparison cmp_;
- };
-
- Chain chain;
- ConsistencyIndex consistency;
+ using Parent::cycle;
+ Chain chain;
};
/**
@@ -194,28 +201,31 @@
*/
template<class Data_ = Empty<>,
class ChainTraits_ = VectorChains<>,
- class Comparison_ = GreaterComparison<>,
- class ContainerTraits_ = VectorContainer<>,
- class ConsistencyIndex_ = size_t,
- class ConsistencyComparison_ = std::less<ConsistencyIndex_>,
- class Element_ = ChainData<Data_, ChainTraits_, ContainerTraits_, ConsistencyIndex_> >
+ class ContainerTraits_ = OrderConsistencyContainer<>,
+ class Element_ = ChainData<Data_, ChainTraits_>,
+ class Comparison_ = ElementComparison<typename ContainerTraits_::template rebind<Element_>::other::Container,
+ std::greater<typename ContainerTraits_::template rebind<Element_>::other::Container::iterator> >,
+ class ConsistencyComparison_ = ElementComparison<typename ContainerTraits_::template rebind<Element_>::other::ConsistentContainer,
+ std::greater<typename ContainerTraits_::template rebind<Element_>::other::ConsistentContainer::iterator> >
+ >
class DynamicPersistenceChains:
- public StaticPersistence<Data_, ChainTraits_, Comparison_, ContainerTraits_, Element_>
+ public StaticPersistence<Data_, ChainTraits_, ContainerTraits_, Element_, Comparison_>
{
public:
typedef Data_ Data;
typedef Element_ Element;
- typedef StaticPersistence<Data_, ChainTraits_, Comparison_,
- ContainerTraits_, Element_> Parent;
+ typedef StaticPersistence<Data_, ChainTraits_,
+ ContainerTraits_, Element_, Comparison_> Parent;
typedef typename Parent::ContainerTraits Traits;
typedef typename Parent::Order Order;
typedef typename Parent::OrderComparison OrderComparison;
typedef typename Parent::OrderIndex OrderIndex;
- typedef ConsistencyIndex_ ConsistencyIndex;
- typedef ThreeOutcomeCompare<
- typename Element::
- template ConsistencyComparison<ConsistencyComparison_> > ConsistencyComparison;
+ typedef ConsistencyComparison_ ConsistencyComparison;
+ typedef typename Parent::iterator iterator;
+
+ typedef typename Element::Chain Chain;
+ typedef typename Element::Cycle Cycle;
/**
* Constructor: DynamicPersistenceChains()
@@ -224,9 +234,7 @@
* Template parameters:
* Filtration - filtration of the complex whose persistence we are computing
*/
- template<class Filtration> DynamicPersistenceChains(const Filtration& f,
- const OrderComparison& ocmp = OrderComparison(),
- const ConsistencyComparison& ccmp = ConsistencyComparison());
+ template<class Filtration> DynamicPersistenceChains(const Filtration& f);
void pair_simplices();
@@ -238,10 +246,12 @@
// TODO: the main missing piece to be dynamic
//template<class Visitor>
- //bool transpose(OrderIndex i, const Visitor& visitor = Visitor());
-
+ //bool transpose(OrderIndex i, Visitor& visitor = Visitor());
+
using Parent::begin;
using Parent::end;
+ using Parent::iterator_to;
+ using Parent::index;
using Parent::size;
// Struct: TranspositionVisitor
@@ -253,17 +263,22 @@
// This function is called before transposition is processed
// (at the very beginning of <transpose(i, visitor)>). It is meant to update any structures
// that may need to be updated, but perhaps it has other uses as well.
- void transpose(OrderIndex i) const {}
+ void transpose(iterator i) const {}
// Function: switched(i, type)
// This function is called after the transposition if the switch in pairing has occured.
// `i` is the index of the preceding simplex after the transposition.
// `type` indicates the <SwitchType>.
- void switched(OrderIndex i, SwitchType type) const {}
+ void switched(iterator i, SwitchType type) const {}
};
protected:
using Parent::order;
+ using Parent::set_pair;
+ using Parent::swap_cycle;
+
+ void cycle_add(iterator i, const Cycle& z) { order().modify(i, boost::bind(&Element::template cycle_add<ConsistencyComparison>, bl::_1, boost::ref(z), ccmp_)); } // i->cycle_add(z, ccmp_)
+ void chain_add(iterator i, const Chain& c) { order().modify(i, boost::bind(&Element::template chain_add<ConsistencyComparison>, bl::_1, boost::ref(c), ccmp_)); } // i->chain_add(c, ccmp_)
private:
void swap(OrderIndex i, OrderIndex j);
@@ -271,15 +286,14 @@
struct PairingChainsVisitor: public Parent::PairVisitor
{
- // TODO: this is specialized for std::vector
- PairingChainsVisitor(OrderIndex bg, ConsistencyComparison ccmp, unsigned size):
- Parent::PairVisitor(size), bg_(bg), ccmp_(ccmp) {}
+ PairingChainsVisitor(Order& order, ConsistencyComparison ccmp, unsigned size):
+ Parent::PairVisitor(size), order_(order), ccmp_(ccmp) {}
- void init(OrderIndex i) const { i->consistency = i - bg_; i->chain.append(i, ccmp_); }
- void update(OrderIndex j, OrderIndex i) const { j->chain.add(i->pair->chain, ccmp_); }
- void finished(OrderIndex i) const { Parent::PairVisitor::finished(i); CountBy(cChainLength, i->chain.size()); }
+ void init(iterator i) const { order_.modify(i, boost::bind(&Element::template chain_append<ConsistencyComparison>, bl::_1, &*i, ccmp_)); } // i->chain_append(&*i, ccmp)
+ void update(iterator j, iterator i) const { order_.modify(j, boost::bind(&Element::template chain_add<ConsistencyComparison>, bl::_1, i->pair->chain, ccmp_)); } // j->chain.add(i->pair->chain, ccmp_)
+ void finished(iterator i) const { Parent::PairVisitor::finished(i); CountBy(cChainLength, i->chain.size()); }
- OrderIndex bg_;
+ Order& order_;
ConsistencyComparison ccmp_;
};
--- a/include/topology/dynamic-persistence.hpp Thu Dec 03 13:21:23 2009 -0800
+++ b/include/topology/dynamic-persistence.hpp Sun Dec 20 08:27:01 2009 -0800
@@ -22,48 +22,45 @@
/* Trails */
-template<class D, class CT, class Cmp, class OT, class CI, class CC, class E>
+template<class D, class CT, class OT, class E, class Cmp, class CCmp>
template<class Filtration>
-DynamicPersistenceTrails<D,CT,Cmp,OT,CI,CC,E>::
-DynamicPersistenceTrails(const Filtration& f, const OrderComparison& ocmp, const ConsistencyComparison& ccmp):
- Parent(f, ocmp), ccmp_(ccmp)
+DynamicPersistenceTrails<D,CT,OT,E,Cmp,CCmp>::
+DynamicPersistenceTrails(const Filtration& f):
+ Parent(f), ccmp_(order().get<consistency>())
{}
-template<class D, class CT, class Cmp, class OT, class CI, class CC, class E>
+template<class D, class CT, class OT, class E, class Cmp, class CCmp>
void
-DynamicPersistenceTrails<D,CT,Cmp,OT,CI,CC,E>::
+DynamicPersistenceTrails<D,CT,OT,E,Cmp,CCmp>::
pair_simplices()
{
- Parent::pair_simplices(begin(), end(), true, PairingTrailsVisitor(begin(), ccmp_, size()));
+ Parent::pair_simplices(begin(), end(), true, PairingTrailsVisitor(order(), ccmp_, size()));
}
-template<class D, class CT, class Cmp, class OT, class CI, class CC, class E>
-template<class Visitor>
+template<class D, class CT, class OT, class E, class Cmp, class CCmp>
+template<class DimensionFunctor, class Visitor>
bool
-DynamicPersistenceTrails<D,CT,Cmp,OT,CI,CC,E>::
-transpose(OrderIndex i, const Visitor& visitor)
+DynamicPersistenceTrails<D,CT,OT,E,Cmp,CCmp>::
+transpose(iterator i, const DimensionFunctor& dimension, Visitor visitor)
{
#if LOGGING
typename Traits::OutputMap outmap(order());
#endif
Count(cTransposition);
- typedef OrderIndex Index;
typedef typename Element::Trail::iterator TrailIterator;
visitor.transpose(i);
- Index i_prev = i++;
+ iterator i_prev = i++;
-#if 0 // Persistence no longer has the notion of dimension
- if (i_prev->dimension() != i->dimension())
+ if (dimension(i_prev) != dimension(i))
{
swap(i_prev, i);
rLog(rlTranspositions, "Different dimension");
Count(cTranspositionDiffDim);
return false;
}
-#endif
bool si = i_prev->sign(), sii = i->sign();
if (si && sii)
@@ -71,15 +68,15 @@
rLog(rlTranspositions, "Trail prev: %s", i_prev->trail.tostring(outmap).c_str());
// Case 1
- boost::optional<TrailIterator> i_in_i_prev = i_prev->trail.contains(i);
- if (i_in_i_prev)
- {
+ if (trail_remove_if_contains(i_prev, index(i)))
rLog(rlTranspositions, "Case 1, U[i,i+1] = 1");
- i_prev->trail.remove(*i_in_i_prev);
- }
- Index k = i_prev->pair;
- Index l = i->pair;
+ iterator k = iterator_to(i_prev->pair);
+ iterator l = iterator_to(i->pair);
+
+ // rLog(rlTranspositions, "(i_prev, k), (i, l): (%s, %s), (%s, %s)",
+ // outmap(i_prev).c_str(), outmap(k).c_str(),
+ // outmap(i).c_str(), outmap(l).c_str());
// Explicit treatment of unpaired simplex
if (l == i)
@@ -91,7 +88,7 @@
return false;
} else if (k == i_prev)
{
- if (!(l->cycle.contains(i_prev)))
+ if (!(l->cycle.contains(index(i_prev))))
{
// Case 1.2
swap(i_prev, i);
@@ -104,7 +101,7 @@
// Case 1.2 --- special version (plain swap, but pairing switches)
swap(i_prev, i);
pairing_switch(i_prev, i);
- visitor.switched(i_prev, Case12);
+ visitor.switched(i, Case12);
rLog(rlTranspositions, "Case 1.2 --- unpaired (pairing switch)");
rLog(rlTranspositions, outmap(i_prev).c_str());
Count(cTranspositionCase12s);
@@ -113,7 +110,7 @@
}
rLog(rlTranspositions, "l cycle: %s", l->cycle.tostring(outmap).c_str());
- if (!(l->cycle.contains(i_prev)))
+ if (!(l->cycle.contains(index(i_prev))))
{
// Case 1.2
swap(i_prev, i);
@@ -123,12 +120,12 @@
} else
{
// Case 1.1
- if (not2(ccmp_)(k,l))
+ if (std::not2(ccmp_)(index(k),index(l)))
{
// Case 1.1.1
swap(i_prev, i);
- l->cycle.add(k->cycle, ccmp_); // Add column k to l
- k->trail.add(l->trail, ccmp_); // Add row l to k
+ cycle_add(l, k->cycle); // Add column k to l
+ trail_add(k, l->trail); // Add row l to k
rLog(rlTranspositions, "Case 1.1.1");
Count(cTranspositionCase111);
return false;
@@ -136,10 +133,10 @@
{
// Case 1.1.2
swap(i_prev, i);
- k->cycle.add(l->cycle, ccmp_); // Add column l to k
- l->trail.add(k->trail, ccmp_); // Add row k to l
+ cycle_add(k, l->cycle); // Add column l to k
+ trail_add(l, k->trail); // Add row k to l
pairing_switch(i_prev, i);
- visitor.switched(i_prev, Case112);
+ visitor.switched(i, Case112);
rLog(rlTranspositions, "Case 1.1.2");
Count(cTranspositionCase112);
return true;
@@ -148,7 +145,7 @@
} else if (!si && !sii)
{
// Case 2
- if (!(i_prev->trail.contains(i)))
+ if (!(i_prev->trail.contains(index(i))))
{
// Case 2.2
swap(i_prev, i);
@@ -158,18 +155,18 @@
} else
{
// Case 2.1
- Index low_i = i_prev->pair;
- Index low_ii = i->pair;
- i_prev->trail.add(i->trail, ccmp_); // Add row i to i_prev
- i->cycle.add(i_prev->cycle, ccmp_); // Add column i_prev to i
+ iterator low_i = iterator_to(i_prev->pair);
+ iterator low_ii = iterator_to(i->pair);
+ trail_add(i_prev, i->trail); // Add row i to i_prev
+ cycle_add(i, i_prev->cycle); // Add column i_prev to i
swap(i_prev, i);
- if (not2(ccmp_)(low_ii, low_i))
+ if (std::not2(ccmp_)(index(low_ii), index(low_i)))
{
// Case 2.1.2
- i_prev->cycle.add(i->cycle, ccmp_); // Add column i to i_prev (after transposition)
- i->trail.add(i_prev->trail, ccmp_); // Add row i to i_prev
+ cycle_add(i_prev, i->cycle); // Add column i to i_prev (after transposition)
+ trail_add(i, i_prev->trail); // Add row i to i_prev
pairing_switch(i_prev, i);
- visitor.switched(i_prev, Case212);
+ visitor.switched(i, Case212);
rLog(rlTranspositions, "Case 2.1.2");
Count(cTranspositionCase212);
return true;
@@ -183,7 +180,7 @@
} else if (!si && sii)
{
// Case 3
- if (!(i_prev->trail.contains(i)))
+ if (!(i_prev->trail.contains(index(i))))
{
// Case 3.2
swap(i_prev, i);
@@ -193,13 +190,13 @@
} else
{
// Case 3.1
- i_prev->trail.add(i->trail, ccmp_); // Add row i to i_prev
- i->cycle.add(i_prev->cycle, ccmp_); // Add column i_prev to i
+ trail_add(i_prev, i->trail); // Add row i to i_prev
+ cycle_add(i, i_prev->cycle); // Add column i_prev to i
swap(i_prev, i);
- i_prev->cycle.add(i->cycle, ccmp_); // Add column i_prev to i (after transposition)
- i->trail.add(i_prev->trail, ccmp_); // Add row i to i_prev
+ cycle_add(i_prev, i->cycle); // Add column i_prev to i (after transposition)
+ trail_add(i, i_prev->trail); // Add row i to i_prev
pairing_switch(i_prev, i);
- visitor.switched(i_prev, Case31);
+ visitor.switched(i, Case31);
rLog(rlTranspositions, "Case 3.1");
Count(cTranspositionCase31);
return true;
@@ -207,12 +204,8 @@
} else if (si && !sii)
{
// Case 4
- boost::optional<TrailIterator> i_in_i_prev = i_prev->trail.contains(i);
- if (i_in_i_prev)
- {
+ if (trail_remove_if_contains(i_prev, index(i)))
rLog(rlTranspositions, "Case 4, U[i,i+1] = 1");
- i_prev->trail.remove(*i_in_i_prev);
- }
swap(i_prev, i);
rLog(rlTranspositions, "Case 4");
Count(cTranspositionCase4);
@@ -222,58 +215,69 @@
return false; // to avoid compiler complaints; we should never reach this point
}
-template<class D, class CT, class Cmp, class OT, class CI, class CC, class E>
+template<class D, class CT, class OT, class E, class Cmp, class CCmp>
void
-DynamicPersistenceTrails<D,CT,Cmp,OT,CI,CC,E>::
-swap(OrderIndex i, OrderIndex j)
+DynamicPersistenceTrails<D,CT,OT,E,Cmp,CCmp>::
+swap(iterator i, iterator j)
{
- std::swap<Data>(*i, *j);
- std::swap(i->pair, j->pair);
-
- std::swap(i->cycle, j->cycle); // TODO: double-check that the STL container specializations actually get invoked
- std::swap(i->trail, j->trail);
+ order().relocate(i,j);
}
-template<class D, class CT, class Cmp, class OT, class CI, class CC, class E>
+template<class D, class CT, class OT, class E, class Cmp, class CCmp>
void
-DynamicPersistenceTrails<D,CT,Cmp,OT,CI,CC,E>::
-pairing_switch(OrderIndex i, OrderIndex j)
+DynamicPersistenceTrails<D,CT,OT,E,Cmp,CCmp>::
+pairing_switch(iterator i, iterator j)
{
OrderIndex i_pair = i->pair;
OrderIndex j_pair = j->pair;
- if (i_pair == i)
- j->pair = j;
+ // rLog(rlTranspositions, " (%x %x %x) (%x %x %x)", &*i, i->pair, i->pair->pair, &*j, j->pair, j->pair->pair);
+ if (i_pair == index(i))
+ set_pair(j, j);
else
{
- j->pair = i_pair;
- i_pair->pair = j;
+ set_pair(j, i_pair);
+ set_pair(i_pair, j);
}
- if (j_pair == j)
- i->pair = i;
+ if (j_pair == index(j))
+ set_pair(i, i);
else
{
- i->pair = j_pair;
- j_pair->pair = i;
+ set_pair(i, j_pair);
+ set_pair(j_pair, i);
}
+ // rLog(rlTranspositions, " (%x %x %x) (%x %x %x)", &*i, i->pair, i->pair->pair, &*j, j->pair, j->pair->pair);
}
+// Helper classes
+template<class D, class CT, class OT, class E, class Cmp, class CCmp>
+struct DynamicPersistenceTrails<D,CT,OT,E,Cmp,CCmp>::TrailRemover:
+ public std::unary_function<Element&, void>
+{
+ TrailRemover(OrderIndex i):
+ i_(i) {}
+
+ void operator()(Element& e) { result = e.trail.remove_if_contains(i_); }
+
+ OrderIndex i_;
+ bool result;
+};
+
/* Chains */
-template<class D, class CT, class Cmp, class OT, class CI, class CC, class E>
+template<class D, class CT, class OT, class E, class Cmp, class CCmp>
template<class Filtration>
-DynamicPersistenceChains<D,CT,Cmp,OT,CI,CC,E>::
-DynamicPersistenceChains(const Filtration& f, const OrderComparison& ocmp, const ConsistencyComparison& ccmp):
- Parent(f, ocmp), ccmp_(ccmp)
+DynamicPersistenceChains<D,CT,OT,E,Cmp,CCmp>::
+DynamicPersistenceChains(const Filtration& f):
+ Parent(f), ccmp_(order().get<consistency>())
{}
-template<class D, class CT, class Cmp, class OT, class CI, class CC, class E>
+template<class D, class CT, class OT, class E, class Cmp, class CCmp>
void
-DynamicPersistenceChains<D,CT,Cmp,OT,CI,CC,E>::
+DynamicPersistenceChains<D,CT,OT,E,Cmp,CCmp>::
pair_simplices()
{
- Parent::pair_simplices(begin(), end(), true, PairingChainsVisitor(begin(), ccmp_, size()));
+ Parent::pair_simplices(begin(), end(), true, PairingChainsVisitor(order(), ccmp_, size()));
}
-
--- a/include/topology/filtration.h Thu Dec 03 13:21:23 2009 -0800
+++ b/include/topology/filtration.h Sun Dec 20 08:27:01 2009 -0800
@@ -3,9 +3,25 @@
#include <vector>
#include <iostream>
+
#include "complex-traits.h"
+
#include "utilities/indirect.h"
#include "utilities/property-maps.h"
+#include "utilities/types.h"
+
+#include <boost/multi_index_container.hpp>
+#include <boost/multi_index/ordered_index.hpp>
+#include <boost/multi_index/identity.hpp>
+#include <boost/multi_index/random_access_index.hpp>
+
+#include <boost/serialization/access.hpp>
+#include <boost/serialization/nvp.hpp>
+#include <boost/serialization/serialization.hpp>
+
+
+namespace b = boost;
+namespace bmi = boost::multi_index;
// Class: Filtration
@@ -13,61 +29,69 @@
// Filtration keeps track of the ordering of the simplices in a complex.
// The most significant function it provides is <boundary()> which converts
// the boundary of a simplex at a given index into a list of indices.
-//
-// TODO: this is really specialized for an std::vector<> Complex; eventually generalize
-// TODO: should we derive from Order?
-template<class Complex_,
- class Index_ = size_t,
- class ComplexTraits_ = ComplexTraits<Complex_> >
+template<class Simplex_,
+ class SimplexOrderIndex_ = bmi::ordered_unique<bmi::identity<Simplex_>,
+ typename Simplex_::VertexComparison> >
class Filtration
{
+ private:
+ struct order {}; // tag
+
public:
// Typedefs: Template parameters
- typedef Index_ IntermediateIndex;
- typedef Complex_ Complex;
- typedef ComplexTraits_ ComplexTraits;
+ typedef Simplex_ Simplex;
+ typedef SimplexOrderIndex_ SimplexOrderIndex;
- // Typedefs: Complex
- typedef typename ComplexTraits::Index ComplexIndex;
- typedef typename ComplexTraits::Simplex Simplex;
- typedef typename ComplexTraits::SimplexIndexMap SimplexIndexMap;
- typedef std::vector<IntermediateIndex> IndexBoundary;
+ typedef b::multi_index_container<Simplex,
+ bmi::indexed_by<SimplexOrderIndex,
+ bmi::random_access<bmi::tag<order> >
+ >
+ > Container;
+ typedef typename Container::value_type value_type;
- // Typedefs: Order
- typedef std::vector<ComplexIndex> Order;
+ // Typedefs: Complex and Order views
+ typedef typename Container::template nth_index<0>::type Complex;
+ typedef typename Container::template nth_index<1>::type Order;
typedef typename Order::const_iterator Index;
- typedef std::vector<IntermediateIndex> ReverseOrder;
- typedef typename ReverseOrder::const_iterator ReverseOrderIndex;
+
+ Filtration() {}
// Constructor: Filtration(bg, end, cmp)
- template<class Comparison>
- Filtration(ComplexIndex bg, ComplexIndex end, const Comparison& cmp = Comparison());
-
- const Simplex& simplex(Index i) const { return **i; }
+ template<class ComplexIndex, class Comparison>
+ Filtration(ComplexIndex bg, ComplexIndex end, const Comparison& cmp = Comparison()):
+ container_(bg, end) { sort(cmp); }
- // Function: boundary(i, bdry, map)
- // Computes boundary of a given index `i` in terms of other indices
- template<class Cycle, class Map>
- void boundary(const Index& i, Cycle& bdry, const Map& map) const;
+ // Lookup
+ const Simplex& simplex(Index i) const { return *i; }
+ Index find(const Simplex& s) const { return bmi::project<order>(container_, container_.find(s)); }
+
+ // Modifiers
+ template<class Comparison>
+ void sort(const Comparison& cmp = Comparison()) { container_.get<order>().sort(cmp); }
+ void push_back(const Simplex& s) { container_.get<order>().push_back(s); }
+ void transpose(Index i) { container_.get<order>().relocate(i, i+1); }
+ void clear() { container_.get<order>().clear(); }
- Index begin() const { return order_.begin(); }
- Index end() const { return order_.end(); }
- size_t size() const { return order_.size(); }
+ Index begin() const { return container_.get<order>().begin(); }
+ Index end() const { return container_.get<order>().end(); }
+ size_t size() const { return container_.size(); }
- std::ostream& operator<<(std::ostream& out) const;
+ std::ostream& operator<<(std::ostream& out) const { std::copy(begin(), end(), std::ostream_iterator<Simplex>(out, "\n")); return out; }
private:
- Order order_;
- ReverseOrder reverse_order_;
- OffsetMap<ComplexIndex,
- ReverseOrderIndex>
- complex_order_map_;
- SimplexIndexMap simplex_index_map_;
+ Container container_;
+
+ private:
+ // Serialization
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int)
+ { ar & boost::serialization::make_nvp("order", container_); }
};
-template<class C, class I, class CT>
+template<class S, class SOI>
std::ostream&
-operator<<(std::ostream& out, const Filtration<C,I,CT>& f) { return f.operator<<(out); }
+operator<<(std::ostream& out, const Filtration<S,SOI>& f) { return f.operator<<(out); }
template<class Functor_, class Filtration_>
@@ -97,6 +121,29 @@
evaluate_through_filtration(const Filtration& filtration, const Functor& functor)
{ return ThroughFiltration<Functor, Filtration>(filtration, functor); }
+
+template<class Map, class Filtration>
+class DimensionFunctor
+{
+ public:
+ DimensionFunctor(const Map& map, const Filtration& filtration):
+ map_(map), filtration_(filtration)
+ {}
+
+ template<class key_type>
+ Dimension operator()(key_type i) const { return filtration_.simplex(map_[i]).dimension(); }
+
+ private:
+ const Map& map_;
+ const Filtration& filtration_;
+};
+
+template<class Map, class Filtration>
+DimensionFunctor<Map, Filtration>
+make_dimension_functor(const Map& map, const Filtration& filtration)
+{ return DimensionFunctor<Map, Filtration>(map, filtration); }
+
+
#include "filtration.hpp"
#endif // __FILTRATION_H__
--- a/include/topology/filtration.hpp Thu Dec 03 13:21:23 2009 -0800
+++ b/include/topology/filtration.hpp Sun Dec 20 08:27:01 2009 -0800
@@ -1,46 +1,6 @@
-#include "utilities/containers.h"
-#include <boost/iterator/counting_iterator.hpp>
-
-template<class C, class I, class CT>
-template<class Comparison>
-Filtration<C,I,CT>::
-Filtration(ComplexIndex bg, ComplexIndex end, const Comparison& cmp):
- order_(boost::counting_iterator<ComplexIndex>(bg),
- boost::counting_iterator<ComplexIndex>(end)),
- reverse_order_(order_.size()),
- complex_order_map_(bg, reverse_order_.begin()),
- simplex_index_map_(bg, end)
-{
- std::sort(order_.begin(), order_.end(), make_indirect_comparison(cmp));
- for (Index obg = order_.begin(), cur = obg; cur != order_.end(); ++cur)
- reverse_order_[*cur - bg] = cur - obg;
-}
+#include <utilities/log.h>
-template<class C, class I, class CT>
-template<class Cycle, class Map>
-void
-Filtration<C,I,CT>::
-boundary(const Index& i, Cycle& bdry, const Map& map) const
-{
- AssertMsg(bdry.empty(), "We are initializing the boundary from scratch");
- ContainerTraits<Cycle>::reserve(bdry, (*i)->boundary_end() - (*i)->boundary_begin());
- typename Map::template rebind_from<IntermediateIndex>::other order_bdry_map(0, map.to());
-
- for (typename Simplex::BoundaryIterator cur = (*i)->boundary_begin(); cur != (*i)->boundary_end(); ++cur)
- {
- //std::cout << *cur << std::endl;
- //std::cout << simplex_index_map_[*cur] - complex_order_map_.from() << std::endl;
- bdry.push_back(order_bdry_map[*complex_order_map_[simplex_index_map_[*cur]]]);
- //std::cout << bdry.back() - order_bdry_map.to() << std::endl;
- }
-}
-
-template<class C, class I, class CT>
-std::ostream&
-Filtration<C,I,CT>::
-operator<<(std::ostream& out) const
-{
- for (Index i = begin(); i != end(); ++i)
- out << **i << std::endl;
- return out;
-}
+#ifdef LOGGING
+static rlog::RLogChannel* rlFiltration = DEF_CHANNEL("topology/filtration/info", rlog::Log_Debug);
+static rlog::RLogChannel* rlFiltrationDebug = DEF_CHANNEL("topology/filtration/debug", rlog::Log_Debug);
+#endif // LOGGING
--- a/include/topology/lowerstarfiltration.h Thu Dec 03 13:21:23 2009 -0800
+++ b/include/topology/lowerstarfiltration.h Sun Dec 20 08:27:01 2009 -0800
@@ -13,21 +13,30 @@
/**
* Struct: MaxVertexComparison
*
- * Functor that determines which simplex has a higher vertex with respect to VertexComparison_
+ * Functor that determines which simplex has a higher vertex with respect to VertexComparison_, breaking ties by dimension
*/
template<class Simplex_, class VertexComparison_>
struct MaxVertexComparison
{
typedef VertexComparison_ VertexComparison;
typedef Simplex_ Simplex;
+ typedef typename Simplex::Vertex Vertex;
MaxVertexComparison(const VertexComparison& vcmp):
vcmp_(vcmp) {}
bool operator()(const Simplex& s1, const Simplex& s2) const
{
- return std::max_element(s1.vertices().begin(), s1.vertices().end(), vcmp) <
- std::max_element(s2.vertices().begin(), s2.vertices().end(), vcmp);
+ const Vertex& max1 = *std::max_element(s1.vertices().begin(), s1.vertices().end(), vcmp_);
+ const Vertex& max2 = *std::max_element(s2.vertices().begin(), s2.vertices().end(), vcmp_);
+
+ bool less = vcmp_(max1, max2),
+ more = vcmp_(max2, max1);
+
+ if (!less && !more) // equal
+ return s1.dimension() < s2.dimension();
+
+ return less;
}
VertexComparison vcmp_;
@@ -37,17 +46,18 @@
/**
* Map from i-th vertex to its index in the filtration.
*/
-template<class Index_>
+template<class Index_, class Filtration_>
class VertexSimplexMap
{
public:
typedef Index_ Index;
- typedef std::vector<FiltrationIndex> VertexVector;
+ typedef Filtration_ Filtration;
+ typedef std::vector<Index> VertexVector;
- VertexSimplexMap(Index begin, Index end, const Map& m)
+ VertexSimplexMap(Index begin, Index end, const Filtration& f)
{
- for (FiltrationIndex cur = begin; cur != end; ++cur)
- if (m[cur].dimension() == 0)
+ for (Index cur = begin; cur != end; ++cur)
+ if (f.simplex(cur).dimension() == 0)
vertices_.push_back(cur);
}
@@ -55,6 +65,5 @@
VertexVector vertices_;
};
-// TODO: transpose_vertices(Index, Filtration, Persistence, Visitor);
#endif // __LOWERSTARFILTRATION_H__
--- a/include/topology/order.h Thu Dec 03 13:21:23 2009 -0800
+++ b/include/topology/order.h Sun Dec 20 08:27:01 2009 -0800
@@ -3,48 +3,111 @@
#include "utilities/types.h"
#include "utilities/indirect.h"
+#include "utilities/property-maps.h"
#include <vector>
+#include <list>
+
+#include <boost/multi_index_container.hpp>
+#include <boost/multi_index/random_access_index.hpp>
+namespace bmi = boost::multi_index;
//#include <iostream>
#include <sstream>
#include <string>
-template<class Element_ = Empty<> >
-struct VectorContainer
+/* Tags */
+// TODO: pollutes global namespace; find a place to localize
+struct order {};
+struct consistency {};
+
+template<class Container>
+class OffsetOutputMap
{
- typedef Element_ Element;
- typedef std::vector<Element> Container;
-
- typedef typename Container::iterator Index;
+ public:
+ typedef const typename Container::value_type* const_element_pointer;
+ typedef typename Container::iterator iterator;
+
+ OffsetOutputMap(const Container& order):
+ om(order.begin(),0), order_(order) {}
- class OutputMap
- {
- public:
- OutputMap(const Container& order):
- bg_(order.begin()) {}
+ // returns a string with (i - bg_)
+ std::string operator()(iterator i) const
+ {
+ std::stringstream s;
+ s << om[i];
+ return s.str();
+ }
+
+ std::string operator()(const_element_pointer p) const
+ {
+ return (*this)(order_.iterator_to(*p));
+ }
- // returns a string with (i - bg_)
- std::string operator()(Index i) const
- {
- std::stringstream s; s << (i - bg_);
- return s.str();
- }
+ private:
+ OffsetMap<typename Container::iterator, unsigned> om;
+ const Container& order_;
+};
- private:
- typename Container::const_iterator bg_;
- };
+template<class Element_ = Empty<> >
+struct OrderContainer
+{
+ typedef Element_ Element;
+ typedef boost::multi_index_container<Element,
+ bmi::indexed_by<
+ bmi::random_access<bmi::tag<order> > /* order index */
+ >
+ > Container;
+ typedef typename Container::template index<order>::type OrderedContainer;
+
+ typedef OffsetOutputMap<Container> OutputMap;
template<class U> struct rebind
- { typedef VectorContainer<U> other; };
+ { typedef OrderContainer<U> other; };
+};
+
+
+template<class Container_, class Comparison_>
+struct ElementComparison: public std::binary_function<const typename Container_::value_type*,
+ const typename Container_::value_type*,
+ bool>
+{
+ typedef Container_ Container;
+ typedef Comparison_ Comparison;
+ typedef typename Container::value_type Element;
+
+ ElementComparison(const Container& container,
+ const Comparison& cmp = Comparison()):
+ container_(container), cmp_(cmp) {}
+
+ bool operator()(const Element* a, const Element* b) const { return cmp_(container_.iterator_to(*a), container_.iterator_to(*b)); }
+
+ const Container& container_;
+ const Comparison& cmp_;
};
-template<class Index = int>
-struct GreaterComparison: public std::greater<Index>
+
+
+template<class Element_ = Empty<> >
+struct OrderConsistencyContainer
{
+ typedef Element_ Element;
+ typedef boost::multi_index_container<Element,
+ bmi::indexed_by<
+ bmi::random_access<bmi::tag<order> >, /* current index */
+ bmi::random_access<bmi::tag<consistency> > /* original index */
+ >
+ > Container;
+
+ typedef typename Container::template index<order>::type OrderedContainer;
+ typedef typename Container::template index<consistency>::type ConsistentContainer;
+
+ typedef OffsetOutputMap<Container> OutputMap;
+
template<class U> struct rebind
- { typedef GreaterComparison<U> other; };
+ { typedef OrderConsistencyContainer<U> other; };
};
+
#endif // __ORDER_H__
--- a/include/topology/persistence-diagram.hpp Thu Dec 03 13:21:23 2009 -0800
+++ b/include/topology/persistence-diagram.hpp Sun Dec 20 08:27:01 2009 -0800
@@ -58,10 +58,10 @@
boost::optional<Point>
make_point(Iterator i, const Evaluator& evaluator, const Visitor& visitor)
{
- RealType x = evaluator(i);
+ RealType x = evaluator(&*i);
RealType y = Infinity;
- if (i->pair != i)
- y = evaluator(i->pair);
+ if (&*(i->pair) != &*i)
+ y = evaluator(&*(i->pair));
Point p(x,y);
visitor.point(i, p);
@@ -98,7 +98,7 @@
{
boost::optional<typename PDiagram::Point> p = make_point<typename PDiagram::Point>(cur, evaluator, visitor);
if (p)
- diagrams[dimension(cur)].push_back(*p);
+ diagrams[dimension(&*cur)].push_back(*p);
}
}
--- a/include/topology/simplex.h Thu Dec 03 13:21:23 2009 -0800
+++ b/include/topology/simplex.h Sun Dec 20 08:27:01 2009 -0800
@@ -95,7 +95,7 @@
/// \name Vertex manipulation
/// @{
- //bool contains(const Vertex& v) const;
+ bool contains(const Vertex& v) const;
bool contains(const Self& s) const;
void add(const Vertex& v);
template<class Iterator>
--- a/include/topology/simplex.hpp Thu Dec 03 13:21:23 2009 -0800
+++ b/include/topology/simplex.hpp Sun Dec 20 08:27:01 2009 -0800
@@ -57,7 +57,6 @@
return BoundaryIterator(vertices().end(), vertices());
}
-#if 0
template<class V, class T>
bool
Simplex<V,T>::
@@ -67,7 +66,6 @@
typename VertexContainer::const_iterator location = std::lower_bound(vertices().begin(), vertices().end(), v);
return ((location != vertices().end()) && (*location == v));
}
-#endif
template<class V, class T>
bool
--- a/include/topology/static-persistence.h Thu Dec 03 13:21:23 2009 -0800
+++ b/include/topology/static-persistence.h Sun Dec 20 08:27:01 2009 -0800
@@ -5,20 +5,23 @@
#include "cycles.h"
#include "filtration.h"
+#include <boost/ref.hpp>
+#include <boost/lambda/lambda.hpp>
+namespace bl = boost::lambda;
+
#include <utilities/types.h>
#include <boost/progress.hpp>
-template<class Data_, class ChainTraits_, class ContainerTraits_, class Element_ = use_default>
+// Element_ should derive from PairCycleData
+template<class Data_, class ChainTraits_, class Element_ = use_default>
struct PairCycleData: public Data_
{
typedef Data_ Data;
typedef typename if_default<Element_, PairCycleData>::type Element;
- typedef PairCycleData<Data, ChainTraits_, ContainerTraits_, Element> Self;
- typedef typename ContainerTraits_::template rebind<Element>::other ContainerTraits;
- typedef typename ContainerTraits::Index Index;
+ typedef const Element* Index;
typedef typename ChainTraits_::template rebind<Index>::other ChainTraits;
typedef typename ChainTraits::Chain Chain;
typedef Chain Cycle;
@@ -28,6 +31,10 @@
{}
bool sign() const { return cycle.empty(); }
+ bool unpaired() const { return pair == this; }
+
+ void swap_cycle(Cycle& z) { cycle.swap(z); }
+ void set_pair(Index i) { pair = i; }
Index pair;
Cycle cycle;
@@ -46,9 +53,10 @@
*/
template<class Data_ = Empty<>,
class ChainTraits_ = VectorChains<>,
- class Comparison_ = GreaterComparison<>,
- class ContainerTraits_ = VectorContainer<>,
- class Element_ = PairCycleData<Data_, ChainTraits_, ContainerTraits_> >
+ class ContainerTraits_ = OrderContainer<>,
+ class Element_ = PairCycleData<Data_, ChainTraits_>,
+ class Comparison_ = ElementComparison<typename ContainerTraits_::template rebind<Element_>::other::Container,
+ std::greater<typename ContainerTraits_::template rebind<Element_>::other::Container::iterator> > >
class StaticPersistence
{
public:
@@ -61,15 +69,16 @@
template rebind<Element>::other ContainerTraits;
typedef typename ContainerTraits::Container Container;
typedef Container Order;
- typedef typename ContainerTraits::Index OrderIndex;
- typedef typename ContainerTraits::Element OrderElement;
+ typedef typename Element::Index OrderIndex;
+ typedef Element OrderElement;
+ typedef typename Order::iterator iterator;
+
typedef typename ChainTraits_::
template rebind<OrderIndex>::other ChainTraits;
typedef typename ChainTraits::Chain Chain;
typedef Chain Cycle;
- typedef typename Comparison_::
- template rebind<OrderIndex>::other OrderComparison;
+ typedef Comparison_ OrderComparison;
/* Constructor: StaticPersistence()
* TODO: write a description
@@ -77,8 +86,7 @@
* Template parameters:
* Filtration - filtration of the complex whose persistence we are computing
*/
- template<class Filtration> StaticPersistence(const Filtration& f,
- const OrderComparison& ocmp = OrderComparison());
+ template<class Filtration> StaticPersistence(const Filtration& f);
// Function: pair_simplices()
// Compute persistence of the filtration
@@ -88,16 +96,35 @@
// begin() - returns OrderIndex of the first element
// end() - returns OrderIndex of one past the last element
// size() - number of elements in the StaticPersistence
- OrderIndex begin() { return order_.begin(); }
- OrderIndex end() { return order_.end(); }
+ iterator begin() const { return order_.begin(); }
+ iterator end() const { return order_.end(); }
+ iterator iterator_to(OrderIndex i) const { return order_.iterator_to(*i); }
+ OrderIndex index(iterator i) const { return &*i; }
size_t size() const { return order_.size(); }
const OrderComparison& order_comparison() const { return ocmp_; }
+ // A map to extract simplices
+ template<class Filtration> class SimplexMap;
+ template<class Filtration>
+ SimplexMap<Filtration> make_simplex_map(const Filtration& filtration) const { return SimplexMap<Filtration>(*this, filtration); }
+
+ class OrderModifier
+ {
+ public:
+ OrderModifier(Order& order): order_(order) {}
+ template<class Functor>
+ void operator()(iterator i, const Functor& f) { order_.modify(i, f); }
+
+ private:
+ Order& order_;
+ };
+ OrderModifier modifier() { return OrderModifier(order()); }
+
protected:
// Function: pair_simplices(bg, end)
// Compute persistence of the simplices in filtration between bg and end
template<class Visitor>
- void pair_simplices(OrderIndex bg, OrderIndex end, bool store_negative = false, const Visitor& visitor = Visitor());
+ void pair_simplices(iterator bg, iterator end, bool store_negative = false, const Visitor& visitor = Visitor());
// Struct: PairVisitor
// Acts as an archetype and if necessary a base class for visitors passed to <pair_simplices(bg, end, visitor)>.
@@ -107,23 +134,29 @@
// Function: init(i)
// Called after OrderElement pointed to by `i` has been initialized
// (its cycle is set to be its boundary, and pair is set to self, i.e. `i`)
- void init(OrderIndex i) const {}
+ void init(iterator i) const {}
// Function: update(j, i)
// Called after the cycle of `i` has been added to the cycle of `j`,
// this allows the derived class to perform the necessary updates
// (e.g., add `i`'s chain to `j`'s chain)
- void update(OrderIndex j, OrderIndex i) const {}
+ void update(iterator j, iterator i) const {}
// Function: finished(j)
// Called after the processing of `j` is finished.
- void finished(OrderIndex j) const { ++show_progress; }
-
+ void finished(iterator j) const { ++show_progress; }
mutable boost::progress_display
show_progress;
};
const Order& order() const { return order_; }
+ Order& order() { return order_; }
+
+ void set_pair(iterator i, iterator j) { set_pair(i, &*j); }
+ void set_pair(iterator i, OrderIndex j) { order_.modify(i, boost::bind(&OrderElement::set_pair, bl::_1, j)); } // i->set_pair(j)
+ void set_pair(OrderIndex i, iterator j) { set_pair(iterator_to(i), &*j); }
+ void set_pair(OrderIndex i, OrderIndex j) { set_pair(iterator_to(i), j); }
+ void swap_cycle(iterator i, Cycle& z) { order_.modify(i, boost::bind(&OrderElement::swap_cycle, bl::_1, boost::ref(z))); } // i->swap_cycle(z)
private:
Order order_;
--- a/include/topology/static-persistence.hpp Thu Dec 03 13:21:23 2009 -0800
+++ b/include/topology/static-persistence.hpp Sun Dec 20 08:27:01 2009 -0800
@@ -1,8 +1,10 @@
#include <utilities/log.h>
#include <utilities/containers.h>
+#include <utilities/property-maps.h>
+
#include <boost/utility/enable_if.hpp>
#include <boost/utility.hpp>
-#include <utilities/property-maps.h>
+#include <boost/foreach.hpp>
#include <boost/foreach.hpp>
@@ -16,31 +18,34 @@
static Counter* cPersistencePairCycleLength = GetCounter("persistence/pair/cyclelength");
#endif // COUNTERS
-template<class D, class CT, class Cmp, class OT, class E>
+template<class D, class CT, class OT, class E, class Cmp>
template<class Filtration>
-StaticPersistence<D, CT, Cmp, OT, E>::
-StaticPersistence(const Filtration& filtration,
- const OrderComparison& ocmp):
- order_(filtration.size()),
- ocmp_(ocmp)
+StaticPersistence<D, CT, OT, E, Cmp>::
+StaticPersistence(const Filtration& filtration):
+ ocmp_(order_)
{
- OrderIndex ocur = begin();
- OffsetMap<typename Filtration::IntermediateIndex, OrderIndex> om(0, ocur); // TODO: this is customized for std::vector Order
- for (typename Filtration::Index cur = filtration.begin(); cur != filtration.end(); ++cur, ++ocur)
+ order_.assign(filtration.size(), OrderElement());
+ rLog(rlPersistence, "Initializing persistence");
+
+ OffsetMap<typename Filtration::Index, iterator> om(filtration.begin(), begin());
+ for (typename Filtration::Index cur = filtration.begin(); cur != filtration.end(); ++cur)
{
- OrderElement& oe = *ocur;
- oe.pair = ocur;
+ Cycle z;
+ BOOST_FOREACH(const typename Filtration::Simplex& s, std::make_pair(cur->boundary_begin(), cur->boundary_end()))
+ z.push_back(index(om[filtration.find(s)]));
+ z.sort(ocmp_);
- filtration.boundary(cur, oe.cycle, om);
- oe.cycle.sort(ocmp_);
+ iterator ocur = om[cur];
+ swap_cycle(ocur, z);
+ set_pair(ocur, ocur);
}
}
-template<class D, class CT, class Cmp, class OT, class E>
+template<class D, class CT, class OT, class E, class Cmp>
template<class Visitor>
void
-StaticPersistence<D, CT, Cmp, OT, E>::
-pair_simplices(OrderIndex bg, OrderIndex end, bool store_negative, const Visitor& visitor)
+StaticPersistence<D, CT, OT, E, Cmp>::
+pair_simplices(iterator bg, iterator end, bool store_negative, const Visitor& visitor)
{
#if LOGGING
typename ContainerTraits::OutputMap outmap(order_);
@@ -48,13 +53,13 @@
// FIXME: need sane output for logging
rLog(rlPersistence, "Entered: pair_simplices");
- for (OrderIndex j = bg; j != end; ++j)
+ for (iterator j = bg; j != end; ++j)
{
visitor.init(j);
rLog(rlPersistence, "Simplex %s", outmap(j).c_str());
- OrderElement& oe = *j;
- typename OrderElement::Cycle& z = oe.cycle;
+ Cycle z;
+ swap_cycle(j, z);
rLog(rlPersistence, " has boundary: %s", z.tostring(outmap).c_str());
// Sparsify the cycle by removing the negative elements
@@ -68,7 +73,7 @@
}
// --------------------------
- CountNum(cPersistencePairBoundaries, oe.cycle.size());
+ CountNum(cPersistencePairBoundaries, z.size());
Count(cPersistencePair);
while(!z.empty())
@@ -76,33 +81,59 @@
OrderIndex i = z.top(ocmp_); // take the youngest element with respect to the OrderComparison
rLog(rlPersistence, " %s: %s", outmap(i).c_str(), outmap(i->pair).c_str());
// TODO: is this even a meaningful assert?
- AssertMsg(!ocmp_(i, j),
+ AssertMsg(!ocmp_(i, index(j)),
"Simplices in the cycle must precede current simplex: (%s in cycle of %s)",
outmap(i).c_str(), outmap(j).c_str());
// i is not paired, so we pair j with i
- if (i->pair == i)
+ if (iterator_to(i->pair) == iterator_to(i))
{
rLog(rlPersistence, " Pairing %s and %s with cycle %s",
outmap(i).c_str(), outmap(j).c_str(),
z.tostring(outmap).c_str());
- i->pair = j;
- j->pair = i;
- CountNum(cPersistencePairCycleLength, z.size());
- CountBy(cPersistencePairCycleLength, z.size());
+
+ set_pair(i, j);
+ swap_cycle(j, z);
+ set_pair(j, i);
+
+ CountNum(cPersistencePairCycleLength, j->cycle.size());
+ CountBy (cPersistencePairCycleLength, j->cycle.size());
break;
}
// update element
z.add(i->pair->cycle, ocmp_);
- visitor.update(j, i);
+ visitor.update(j, iterator_to(i));
rLog(rlPersistence, " new cycle: %s", z.tostring(outmap).c_str());
}
+ // if z was empty, so is (already) j->cycle, so nothing to do
visitor.finished(j);
rLog(rlPersistence, "Finished with %s: %s",
- outmap(j).c_str(), outmap(oe.pair).c_str());
+ outmap(j).c_str(), outmap(j->pair).c_str());
}
}
+template<class D, class CT, class OT, class E, class Cmp>
+template<class Filtration>
+class StaticPersistence<D, CT, OT, E, Cmp>::SimplexMap
+{
+ public:
+ typedef OrderIndex key_type;
+ // typedef iterator key_type;
+ typedef const typename Filtration::Simplex& value_type;
+
+
+ SimplexMap(const StaticPersistence& persistence,
+ const Filtration& filtration):
+ persistence_(persistence),
+ filtration_(filtration) {}
+
+ value_type operator[](OrderIndex k) const { return (*this)[persistence_.iterator_to(k)]; }
+ value_type operator[](iterator i) const { return filtration_.simplex(filtration_.begin() + (i - persistence_.begin())); }
+
+ private:
+ const StaticPersistence& persistence_;
+ const Filtration& filtration_;
+};
/* Private */
--- a/include/topology/vineyard.h Thu Dec 03 13:21:23 2009 -0800
+++ b/include/topology/vineyard.h Sun Dec 20 08:27:01 2009 -0800
@@ -1,6 +1,6 @@
/*
* Author: Dmitriy Morozov
- * Department of Computer Science, Duke University, 2005 -- 2006
+ * Department of Computer Science, Duke University, 2005 -- 2009
*/
#ifndef __VINEYARD_H__
@@ -14,144 +14,142 @@
#include <boost/serialization/vector.hpp>
#include <boost/serialization/nvp.hpp>
#include <boost/serialization/list.hpp>
-
+
+#include <boost/iterator/iterator_traits.hpp>
-template<class Smplx> class Knee;
-template<class Smplx> class Vine;
+
+class Knee;
+class Vine;
/**
* Vineyard class. Keeps track of vines and knees. switched() is the key function called
- * by filtration when pairing switches after a Filtration::transpose().
+ * when pairing switches.
*
* \ingroup topology
*/
-template<class FltrSmplx>
+template<class Index_, class Iterator_, class Evaluator_>
class Vineyard
{
- public:
- typedef FltrSmplx FiltrationSimplex;
- typedef typename FiltrationSimplex::Simplex Simplex;
- typedef Vine<Simplex> VineS;
- typedef Knee<Simplex> KneeS;
- typedef std::list<VineS> VineList;
- typedef std::list<VineList> VineListList;
- typedef std::vector<typename VineListList::iterator> VineListVector;
- typedef typename FiltrationSimplex::Cycle Cycle;
+ public:
+ typedef Index_ Index;
+ typedef Iterator_ Iterator;
+ typedef Evaluator_ Evaluator;
- typedef typename FiltrationSimplex::Index Index;
- typedef typename FiltrationSimplex::EvaluatorS Evaluator;
-
- public:
- Vineyard(Evaluator* eval = 0):
- evaluator(eval) {}
+ typedef std::list<Vine> VineList;
+ typedef std::list<VineList> VineListList;
+ typedef std::vector<VineListList::iterator> VineListVector;
+
+ public:
+ Vineyard(Evaluator* eval = 0):
+ evaluator(eval) {}
- void start_vines(Index bg, Index end);
- void switched(Index i, Index j);
- void record_knee(Index i);
- void record_diagram(Index bg, Index end);
+ void start_vines(Iterator bg, Iterator end);
+ void switched(Index i, Index j);
+ template<class Iter>
+ void record_knee(Iter i);
+ void record_diagram(Iterator bg, Iterator end);
- void set_evaluator(Evaluator* eval) { evaluator = eval; }
-
- void save_edges(const std::string& filename) const;
+ void set_evaluator(Evaluator* eval) { evaluator = eval; }
- protected:
- typename KneeS::SimplexList resolve_cycle(Index i) const;
+ void save_edges(const std::string& filename, bool skip_infinite = false) const;
+ void save_vines(const std::string& filename, bool skip_infinite = false) const;
- private:
- void start_vine(Index i);
+ private:
+ template<class Iter>
+ void start_vine(Iter i);
- private:
- VineListList vines; // stores vine lists
- VineListVector vines_vector; // stores pointers (iterators) to vine lists
- Evaluator* evaluator;
+ private:
+ VineListList vines; // stores vine lists
+ VineListVector vines_vector; // stores pointers (iterators) to vine lists
+ Evaluator* evaluator;
};
/**
- * Knee class stores the knee in R^3 as well as the cycle that is associated with the Simplex starting from the Knee.
+ * Knee class stores the knee in R^3.
*
* \ingroup topology
*/
-template<class S>
class Knee
{
- public:
- typedef S Simplex;
- typedef std::list<Simplex> SimplexList;
-
- RealType birth;
- RealType death;
- RealType time;
- SimplexList cycle;
-
- // Default parameters for serialization
- Knee(RealType b = 0, RealType d = 0, RealType t = 0):
- birth(b), death(d), time(t)
- {}
- Knee(const Knee& other):
- birth(other.birth), death(other.death), time(other.time)
- {}
+ public:
+ RealType birth;
+ RealType death;
+ RealType time;
+
+ // Default parameters for serialization
+ Knee(RealType b = 0, RealType d = 0, RealType t = 0):
+ birth(b), death(d), time(t)
+ {}
+ Knee(const Knee& other):
+ birth(other.birth), death(other.death), time(other.time)
+ {}
- bool is_diagonal() const { return birth == death; }
- bool is_infinite() const { return (death == Infinity) || (birth == Infinity); }
- void set_cycle(const SimplexList& lst) { cycle = lst; }
+ bool is_diagonal() const { return birth == death; }
+ bool is_infinite() const { return (death == Infinity) || (birth == Infinity); }
- std::ostream& operator<<(std::ostream& out) const { return out << "(" << birth << ", "
- << death << ", "
- << time << ")"; }
-
- private:
- friend class boost::serialization::access;
+ std::ostream& operator<<(std::ostream& out) const { return out << "(" << birth << ", "
+ << death << ", "
+ << time << ")"; }
+
+ private:
+ friend class boost::serialization::access;
- template<class Archive>
- void serialize(Archive& ar, version_type );
+ template<class Archive>
+ void serialize(Archive& ar, version_type );
};
-template<class S>
-std::ostream& operator<<(std::ostream& out, const Knee<S>& k) { return k.operator<<(out); }
+std::ostream& operator<<(std::ostream& out, const Knee& k) { return k.operator<<(out); }
/**
* Vine is a list of Knees
*/
-template<class S>
-class Vine: public std::list<Knee<S> >
-{
- public:
- typedef S Simplex;
- typedef Knee<Simplex> KneeS;
- typedef std::list<KneeS> VineRepresentation;
- typedef typename VineRepresentation::const_iterator const_knee_iterator;
-
- Vine() {}
- Vine(const Vine& other):
- VineRepresentation(other) {}
- Vine(const VineRepresentation& other):
- VineRepresentation(other) {}
- Vine(const KneeS& k) { add(k); }
-
- void add(RealType b, RealType d, RealType t) { push_back(KneeS(b,d,t)); }
- void add(const KneeS& k) { push_back(k); }
+class Vine: public std::list<Knee>
+{
+ public:
+ typedef std::list<Knee> VineRepresentation;
+ typedef VineRepresentation::const_iterator const_knee_iterator;
+
+ Vine() {}
+ Vine(const Vine& other):
+ VineRepresentation(other) {}
+ Vine(const VineRepresentation& other):
+ VineRepresentation(other) {}
+ Vine(const Knee& k) { add(k); }
+
+ void add(RealType b, RealType d, RealType t) { push_back(Knee(b,d,t)); }
+ void add(const Knee& k) { push_back(k); }
- std::ostream& operator<<(std::ostream& out) const { for (const_knee_iterator cur = begin(); cur != end(); ++cur) out << *cur; return out; }
+ std::ostream& operator<<(std::ostream& out) const { std::copy(begin(), end(), std::ostream_iterator<Knee>(out, " ")); return out; }
- using VineRepresentation::begin;
- using VineRepresentation::end;
- using VineRepresentation::front;
- using VineRepresentation::back;
- using VineRepresentation::size;
- using VineRepresentation::empty;
+ using VineRepresentation::begin;
+ using VineRepresentation::end;
+ using VineRepresentation::front;
+ using VineRepresentation::back;
+ using VineRepresentation::size;
+ using VineRepresentation::empty;
- protected:
- using VineRepresentation::push_back;
+ protected:
+ using VineRepresentation::push_back;
- private:
- friend class boost::serialization::access;
+ private:
+ friend class boost::serialization::access;
- template<class Archive>
- void serialize(Archive& ar, version_type );
+ template<class Archive>
+ void serialize(Archive& ar, version_type );
};
-template<class S>
-std::ostream& operator<<(std::ostream& out, const Vine<S>& v) { return v.operator<<(out); }
+std::ostream& operator<<(std::ostream& out, const Vine& v) { return v.operator<<(out); }
+
+
+class VineData
+{
+ public:
+ void set_vine(Vine* vine) const { vine_ = vine; }
+ Vine* vine() const { return vine_; }
+
+ private:
+ mutable Vine* vine_; // cheap trick to work around MultiIndex's constness
+};
#include "vineyard.hpp"
--- a/include/topology/vineyard.hpp Thu Dec 03 13:21:23 2009 -0800
+++ b/include/topology/vineyard.hpp Sun Dec 20 08:27:01 2009 -0800
@@ -6,160 +6,185 @@
#include "utilities/log.h"
#ifdef LOGGING
-static rlog::RLogChannel* rlVineyard = DEF_CHANNEL("topology/vineyard", rlog::Log_Debug);
+static rlog::RLogChannel* rlVineyard = DEF_CHANNEL("topology/vineyard", rlog::Log_Debug);
#endif // LOGGING
-template<class FS>
+template<class I, class It, class E>
void
-Vineyard<FS>::
-start_vines(Index bg, Index end)
+Vineyard<I,It,E>::
+start_vines(Iterator bg, Iterator end)
{
- AssertMsg(evaluator != 0, "Cannot start vines with a null evaluator");
- for (Index cur = bg; cur != end; ++cur)
- {
- if (!cur->sign()) continue;
- Dimension dim = cur->dimension();
-
- if (dim >= vines.size())
- {
- AssertMsg(dim == vines.size(), "New dimension has to be contiguous");
- vines.push_back(VineList());
+ AssertMsg(evaluator != 0, "Cannot start vines with a null evaluator");
+ for (Iterator cur = bg; cur != end; ++cur)
+ {
+ if (!cur->sign()) continue;
+ Dimension dim = evaluator->dimension(cur);
+
+ if (dim >= vines.size())
+ {
+ AssertMsg(dim == vines.size(), "New dimension has to be contiguous");
+ vines.push_back(VineList());
vines_vector.push_back(boost::prior(vines.end()));
- }
+ }
- start_vine(cur);
- record_knee(cur);
- }
+ start_vine(cur);
+ record_knee(cur);
+ }
}
-template<class FS>
+template<class I, class It, class E>
void
-Vineyard<FS>::
+Vineyard<I,It,E>::
switched(Index i, Index j)
{
- VineS* i_vine = i->vine();
- VineS* j_vine = j->vine();
- i->set_vine(j_vine);
- j->set_vine(i_vine);
+ rLog(rlVineyard, "Switching vines");
+
+ Vine* i_vine = i->vine();
+ Vine* j_vine = j->vine();
+ i->set_vine(j_vine);
+ j->set_vine(i_vine);
+
+ // rLog(rlVineyard, " %x %x %x %x", i->pair, i->pair->pair, j->pair, j->pair->pair);
+ // rLog(rlVineyard, " %x %x %x %x", i_vine, i->pair->vine(), j_vine, j->pair->vine());
- // Since the pairing has already been updated, the following assertions should be true
- AssertMsg(i->vine() == i->pair()->vine(), "i's vine must match the vine of its pair");
- AssertMsg(j->vine() == j->pair()->vine(), "j's vine must match the vine of its pair");
+ // Since the pairing has already been updated, the following assertions should be true
+ AssertMsg(i->vine() == i->pair->vine(), "i's vine must match the vine of its pair");
+ AssertMsg(j->vine() == j->pair->vine(), "j's vine must match the vine of its pair");
- if (!i->sign()) i = i->pair();
- if (!j->sign()) j = j->pair();
- record_knee(i);
- record_knee(j);
+ if (!i->sign()) i = i->pair;
+ if (!j->sign()) j = j->pair;
+
+ // std::cout << "i sign: " << i->sign() << std::endl;
+ // std::cout << "j sign: " << j->sign() << std::endl;
+
+ record_knee(i);
+ record_knee(j);
}
-template<class FS>
+template<class I, class It, class E>
+template<class Iter>
void
-Vineyard<FS>::
-start_vine(Index i)
+Vineyard<I,It,E>::
+start_vine(Iter i)
{
- rLog(rlVineyard, "Starting new vine");
- AssertMsg(i->sign(), "Can only start vines for positive simplices");
-
- Dimension dim = i->dimension();
- vines_vector[dim]->push_back(VineS());
- i->set_vine(&vines_vector[dim]->back());
- i->pair()->set_vine(i->vine());
+ rLog(rlVineyard, "Starting new vine");
+ AssertMsg(i->sign(), "Can only start vines for positive simplices");
+
+ Dimension dim = evaluator->dimension(i);
+ vines_vector[dim]->push_back(Vine());
+ i->set_vine(&vines_vector[dim]->back());
+ i->pair->set_vine(i->vine());
}
-
+
/// Records the current diagram in the vineyard
-template<class FS>
+template<class I, class It, class E>
void
-Vineyard<FS>::
-record_diagram(Index bg, Index end)
+Vineyard<I,It,E>::
+record_diagram(Iterator bg, Iterator end)
{
- rLog(rlVineyard, "Entered: record_diagram()");
- AssertMsg(evaluator != 0, "Cannot record diagram with a null evaluator");
-
- for (Index i = bg; i != end; ++i)
- {
- AssertMsg(i->vine() != 0, "Cannot process a null vine in record_diagram");
- if (!i->sign()) continue;
- record_knee(i);
- }
+ rLog(rlVineyard, "Entered: record_diagram()");
+ AssertMsg(evaluator != 0, "Cannot record diagram with a null evaluator");
+
+ for (Iterator i = bg; i != end; ++i)
+ {
+ AssertMsg(i->vine() != 0, "Cannot process a null vine in record_diagram");
+ if (!i->sign()) continue;
+ record_knee(i);
+ }
}
-template<class FS>
-void
-Vineyard<FS>::
-save_edges(const std::string& filename) const
+template<class I, class It, class E>
+void
+Vineyard<I,It,E>::
+save_edges(const std::string& filename, bool skip_infinite) const
{
- for (unsigned int i = 0; i < vines_vector.size(); ++i)
- {
- std::ostringstream os; os << i;
- std::string fn = filename + os.str() + ".edg";
- std::ofstream out(fn.c_str());
- for (typename VineList::const_iterator vi = vines_vector[i]->begin(); vi != vines_vector[i]->end(); ++vi)
- for (typename VineS::const_iterator ki = vi->begin(), kiprev = ki++; ki != vi->end(); kiprev = ki++)
- {
- if (kiprev->is_infinite() || ki->is_infinite()) continue;
- out << kiprev->birth << ' ' << kiprev->death << ' ' << kiprev->time << std::endl;
- out << ki->birth << ' ' << ki->death << ' ' << ki->time << std::endl;
- }
- out.close();
- }
+ for (unsigned int i = 0; i < vines_vector.size(); ++i)
+ {
+ std::ostringstream os; os << i;
+ std::string fn = filename + os.str() + ".edg";
+ std::ofstream out(fn.c_str());
+ for (typename VineList::const_iterator vi = vines_vector[i]->begin(); vi != vines_vector[i]->end(); ++vi)
+ for (typename Vine::const_iterator ki = vi->begin(), kiprev = ki++; ki != vi->end(); kiprev = ki++)
+ {
+ if (skip_infinite && (kiprev->is_infinite() || ki->is_infinite()))
+ {
+ std::cerr << "Warning: skipping an infinite knee in save_edges() in dimension " << i << std::endl;
+ continue;
+ }
+ out << kiprev->birth << ' ' << kiprev->death << ' ' << kiprev->time << std::endl;
+ out << ki->birth << ' ' << ki->death << ' ' << ki->time << std::endl;
+ }
+ out.close();
+ }
+}
+
+template<class I, class It, class E>
+void
+Vineyard<I,It,E>::
+save_vines(const std::string& filename, bool skip_infinite) const
+{
+ for (unsigned int i = 0; i < vines_vector.size(); ++i)
+ {
+ std::ostringstream os; os << i;
+ std::string fn = filename + os.str() + ".vin";
+ std::ofstream out(fn.c_str());
+ for (typename VineList::const_iterator vi = vines_vector[i]->begin(); vi != vines_vector[i]->end(); ++vi)
+ {
+ for (typename Vine::const_iterator ki = vi->begin(); ki != vi->end(); ki++)
+ {
+ if (skip_infinite && ki->is_infinite())
+ {
+ std::cerr << "Warning: skipping an infinite knee in save_edges() in dimension " << i << std::endl;
+ continue;
+ }
+ out << ki->birth << ' ' << ki->death << ' ' << ki->time << " ";
+ }
+ out << std::endl;
+ }
+ out.close();
+ }
}
/// Records a knee for the given simplex
-template<class FS>
+template<class I, class It, class E>
+template<class Iter>
void
-Vineyard<FS>::
-record_knee(Index i)
+Vineyard<I,It,E>::
+record_knee(Iter i)
{
- rLog(rlVineyard, "Entered record_knee()");
- AssertMsg(evaluator != 0, "Cannot record knee with a null evaluator");
- AssertMsg(i->vine() != 0, "Cannot add a knee to a null vine");
- AssertMsg(i->sign(), "record_knee() must be called on a positive simplex");
-
- if (!i->is_paired())
- i->vine()->add(evaluator->value(*i), Infinity, evaluator->time());
- else
- {
- rLog(rlVineyard, "Creating knee");
- KneeS k(evaluator->value(*i), evaluator->value(*(i->pair())), evaluator->time());
- rLog(rlVineyard, "Knee created: %s", tostring(k).c_str());
+ rLog(rlVineyard, "Entered record_knee()");
+ AssertMsg(evaluator != 0, "Cannot record knee with a null evaluator");
+ AssertMsg(i->vine() != 0, "Cannot add a knee to a null vine");
+ AssertMsg(i->sign(), "record_knee() must be called on a positive simplex");
+
+ if (i->unpaired())
+ i->vine()->add((*evaluator)(i), Infinity, evaluator->time());
+ else
+ {
+ rLog(rlVineyard, "Creating knee");
+ Knee k((*evaluator)(i), (*evaluator)((i->pair)), evaluator->time());
+ rLog(rlVineyard, "Knee created: %s", tostring(k).c_str());
rLog(rlVineyard, "Vine: %s", tostring(*(i->vine())).c_str());
- if (!k.is_diagonal() || i->vine()->empty()) // non-diagonal k, or empty vine
- {
- rLog(rlVineyard, "Extending a vine");
- i->vine()->add(k);
- }
- else if (i->vine()->back().is_diagonal()) // last knee is diagonal
- {
- AssertMsg(i->vine()->size() == 1, "Only first knee may be diagonal for a live vine");
- rLog(rlVineyard, "Overwriting first diagonal knee");
- i->vine()->back() = k;
- } else // finish this vine
- {
- rLog(rlVineyard, "Finishing a vine");
- i->vine()->add(k);
- start_vine(i);
- i->vine()->add(k);
- }
- }
-
- i->vine()->back().set_cycle(resolve_cycle(i));
- rLog(rlVineyard, "Leaving record_knee()");
+ if (!k.is_diagonal() || i->vine()->empty()) // non-diagonal k, or empty vine
+ {
+ rLog(rlVineyard, "Extending a vine");
+ i->vine()->add(k);
+ }
+ else if (i->vine()->back().is_diagonal()) // last knee is diagonal
+ {
+ AssertMsg(i->vine()->size() == 1, "Only first knee may be diagonal for a live vine");
+ rLog(rlVineyard, "Overwriting first diagonal knee");
+ i->vine()->back() = k;
+ } else // finish this vine
+ {
+ rLog(rlVineyard, "Finishing a vine");
+ i->vine()->add(k);
+ start_vine(i);
+ i->vine()->add(k);
+ }
+ }
+
+ rLog(rlVineyard, "Leaving record_knee()");
}
-
-template<class FS>
-typename Vineyard<FS>::KneeS::SimplexList
-Vineyard<FS>::
-resolve_cycle(Index i) const
-{
- rLog(rlVineyard, "Entered resolve_cycle");
- const Cycle& cycle = i->cycle();
-
- // Resolve simplices
- typename KneeS::SimplexList lst;
- for (typename Cycle::const_iterator cur = cycle.begin(); cur != cycle.end(); ++cur)
- lst.push_back(**cur);
-
- return lst;
-}
--- a/include/utilities/indirect.h Thu Dec 03 13:21:23 2009 -0800
+++ b/include/utilities/indirect.h Sun Dec 20 08:27:01 2009 -0800
@@ -7,6 +7,9 @@
// TODO: write documentation
+/**************
+ * Comparison *
+ **************/
template<class Comparison_>
struct IndirectComparison
@@ -64,6 +67,27 @@
}
};
+template<class Evaluator_>
+class ThroughEvaluatorComparison
+{
+ public:
+ typedef Evaluator_ Evaluator;
+
+ ThroughEvaluatorComparison(const Evaluator& eval):
+ eval_(eval) {}
+
+ template<class T>
+ bool operator()(T a, T b) const { return (eval_(a) < eval_(b)); }
+
+ private:
+ const Evaluator& eval_;
+};
+
+
+/*************
+ * Iterators *
+ *************/
+
// Iterates over the difference of the two sorted sequences, dereferencing into the first sequence
template<class Iterator1, class Iterator2, class StrictWeakOrdering>
class difference_iterator: public boost::iterator_facade<difference_iterator<Iterator1, Iterator2, StrictWeakOrdering>,
--- a/include/utilities/property-maps.h Thu Dec 03 13:21:23 2009 -0800
+++ b/include/utilities/property-maps.h Sun Dec 20 08:27:01 2009 -0800
@@ -1,7 +1,7 @@
#ifndef __PROPERTY_MAPS_H__
#define __PROPERTY_MAPS_H__
-#include <boost/property_map.hpp>
+#include <boost/property_map/property_map.hpp>
#include <boost/iterator/iterator_traits.hpp>
#include <algorithm>
#include "utilities/log.h"
@@ -106,6 +106,46 @@
{ return OffsetMap<From_, To_>(bg_from, bg_to); }
+template<class From_, class To_, class FromIndex_, class ToIndex_>
+struct OffsetBeginMap
+{
+ typedef From_ From;
+ typedef To_ To;
+ typedef FromIndex_ key_type;
+ typedef ToIndex_ value_type;
+
+ OffsetBeginMap(const From& from, const To& to):
+ from_(from), to_(to) {}
+
+ value_type operator[](key_type i) const { return to_.begin() + (i - from_.begin()); }
+
+ const From& from() const { return from_; }
+ const To& to() const { return to_; }
+
+ private:
+ const From& from_;
+ const To& to_;
+};
+
+
+/* ChainMap */
+template<class Map1, class Map2>
+class ChainMap
+{
+ public:
+ typedef typename Map1::key_type key_type;
+ typedef typename Map2::value_type value_type;
+
+
+ ChainMap(const Map1& m1, const Map2& m2):
+ m1_(m1), m2_(m2) {}
+ value_type operator[](const key_type& k) const { return m2_[m1_[k]]; }
+
+ private:
+ const Map1& m1_;
+ const Map2& m2_;
+};
+
/* ThroughMap */
template<class Functor_, class Map_>
class ThroughMap
@@ -115,7 +155,7 @@
typedef Functor_ Functor;
typedef typename Functor::result_type result_type;
- typedef typename Map::key_type first_argument_type;
+ typedef typename Map_::key_type first_argument_type;
ThroughMap(const Map& map,
const Functor& functor):