# HG changeset patch
# User Dmitriy Morozov <dmitriy@mrzv.org>
# Date 1261006746 28800
# Node ID d15c6d14464532a67dcbfac4942ea8bf713b3ed1
# Parent a99fdaafa31a3bc8b37d507b4d70d4960271a102
Resurrected vineyard code:
* Switched StaticPersistence + (serializable) Filtration to Boost.MultiIndex
* Updated DynamicPersistenceTrails to work with the new MultiIndex way
* Created LSVineyard class, and fixed the grid examples
diff -r a99fdaafa31a -r d15c6d144645 examples/CMakeLists.txt
--- a/examples/CMakeLists.txt Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/CMakeLists.txt Wed Dec 16 15:39:06 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)
diff -r a99fdaafa31a -r d15c6d144645 examples/alphashapes/alphashapes2d.cpp
--- a/examples/alphashapes/alphashapes2d.cpp Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/alphashapes/alphashapes2d.cpp Wed Dec 16 15:39:06 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;
diff -r a99fdaafa31a -r d15c6d144645 examples/alphashapes/alphashapes2d.h
--- a/examples/alphashapes/alphashapes2d.h Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/alphashapes/alphashapes2d.h Wed Dec 16 15:39:06 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); }
diff -r a99fdaafa31a -r d15c6d144645 examples/alphashapes/alphashapes2d.hpp
--- a/examples/alphashapes/alphashapes2d.hpp Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/alphashapes/alphashapes2d.hpp Wed Dec 16 15:39:06 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);
}
diff -r a99fdaafa31a -r d15c6d144645 examples/alphashapes/alphashapes3d.cpp
--- a/examples/alphashapes/alphashapes3d.cpp Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/alphashapes/alphashapes3d.cpp Wed Dec 16 15:39:06 2009 -0800
@@ -13,7 +13,7 @@
#include <boost/program_options.hpp>
-typedef Filtration<AlphaSimplex3DVector> AlphaFiltration;
+typedef Filtration<AlphaSimplex3D> AlphaFiltration;
typedef StaticPersistence<> Persistence;
typedef PersistenceDiagram<> PDgm;
@@ -71,12 +71,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);
@@ -85,12 +85,11 @@
p.pair_simplices();
rInfo("Simplices paired");
+ 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;
diff -r a99fdaafa31a -r d15c6d144645 examples/alphashapes/alphashapes3d.h
--- a/examples/alphashapes/alphashapes3d.h Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/alphashapes/alphashapes3d.h Wed Dec 16 15:39:06 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); }
diff -r a99fdaafa31a -r d15c6d144645 examples/alphashapes/alphashapes3d.hpp
--- a/examples/alphashapes/alphashapes3d.hpp Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/alphashapes/alphashapes3d.hpp Wed Dec 16 15:39:06 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);
}
diff -r a99fdaafa31a -r d15c6d144645 examples/cech-complex/cech-complex.cpp
--- a/examples/cech-complex/cech-complex.cpp Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/cech-complex/cech-complex.cpp Wed Dec 16 15:39:06 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)
{
diff -r a99fdaafa31a -r d15c6d144645 examples/fitness/avida-distance.cpp
--- a/examples/fitness/avida-distance.cpp Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/fitness/avida-distance.cpp Wed Dec 16 15:39:06 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));
diff -r a99fdaafa31a -r d15c6d144645 examples/fitness/avida-rips-distance.cpp
--- a/examples/fitness/avida-rips-distance.cpp Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/fitness/avida-rips-distance.cpp Wed Dec 16 15:39:06 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));
diff -r a99fdaafa31a -r d15c6d144645 examples/grid/CMakeLists.txt
--- a/examples/grid/CMakeLists.txt Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/grid/CMakeLists.txt Wed Dec 16 15:39:06 2009 -0800
@@ -1,5 +1,6 @@
set (targets
test-grid2D
+ test-grid2D-vineyard
combustion-vineyard)
if (use_dsrpdb)
@@ -7,9 +8,13 @@
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})
+endif (CGAL_FOUND)
diff -r a99fdaafa31a -r d15c6d144645 examples/grid/combustion-vineyard.cpp
--- a/examples/grid/combustion-vineyard.cpp Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/grid/combustion-vineyard.cpp Wed Dec 16 15:39:06 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");
}
diff -r a99fdaafa31a -r d15c6d144645 examples/grid/grid2D.h
--- a/examples/grid/grid2D.h Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/grid/grid2D.h Wed Dec 16 15:39:06 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 value_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__
diff -r a99fdaafa31a -r d15c6d144645 examples/grid/grid2D.hpp
--- a/examples/grid/grid2D.hpp Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/grid/grid2D.hpp Wed Dec 16 15:39:06 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);
+ }
+}
diff -r a99fdaafa31a -r d15c6d144645 examples/grid/grid2Dvineyard.h
--- a/examples/grid/grid2Dvineyard.h Tue Jul 28 13:27:13 2009 -0700
+++ /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__
diff -r a99fdaafa31a -r d15c6d144645 examples/grid/grid2Dvineyard.hpp
--- a/examples/grid/grid2Dvineyard.hpp Tue Jul 28 13:27:13 2009 -0700
+++ /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_);
-}
-
diff -r a99fdaafa31a -r d15c6d144645 examples/grid/lsvineyard.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/grid/lsvineyard.h Wed Dec 16 15:39:06 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::value_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:
+ 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:
+ 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<Dimension, iterator>
+{
+ 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__
diff -r a99fdaafa31a -r d15c6d144645 examples/grid/lsvineyard.hpp
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/grid/lsvineyard.hpp Wed Dec 16 15:39:06 2009 -0800
@@ -0,0 +1,251 @@
+#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_
+ std::cout << "Vertex order:" << std::endl;
+ for (VertexIndex cur = vertices_.begin(); cur != vertices_.end(); ++cur)
+ std::cout << " " << cur->vertex() << std::endl;
+
+ // 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());
+ 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_;
+};
diff -r a99fdaafa31a -r d15c6d144645 examples/grid/pdbdistance-vineyard.cpp
--- a/examples/grid/pdbdistance-vineyard.cpp Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/grid/pdbdistance-vineyard.cpp Wed Dec 16 15:39:06 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);
diff -r a99fdaafa31a -r d15c6d144645 examples/grid/test-grid2D-vineyard.cpp
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/grid/test-grid2D-vineyard.cpp Wed Dec 16 15:39:06 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");
+}
diff -r a99fdaafa31a -r d15c6d144645 examples/poincare/poincare.cpp
--- a/examples/poincare/poincare.cpp Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/poincare/poincare.cpp Wed Dec 16 15:39:06 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;
diff -r a99fdaafa31a -r d15c6d144645 examples/rips/rips-pairwise.cpp
--- a/examples/rips/rips-pairwise.cpp Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/rips/rips-pairwise.cpp Wed Dec 16 15:39:06 2009 -0800
@@ -21,9 +21,8 @@
typedef Rips<PairDistances> Generator;
typedef Generator::Simplex Smplx;
-typedef std::vector<Smplx> SimplexVector;
-typedef Filtration<SimplexVector, unsigned> Fltr;
-//typedef StaticPersistence<> Persistence;
+typedef Filtration<Smplx> Fltr;
+// typedef StaticPersistence<> Persistence;
typedef DynamicPersistenceChains<> Persistence;
typedef PersistenceDiagram<> PDgm;
@@ -43,44 +42,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;
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() != 1) continue;
// std::cout << "Unpaired birth: " << size(b) << std::endl;
- // cycle = cur->chain;
+ // cycle = cur->chain; // TODO
if (b.dimension() >= skeleton) continue;
std::cout << b.dimension() << " " << size(b) << " inf" << std::endl;
}
@@ -89,13 +89,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");
}
diff -r a99fdaafa31a -r d15c6d144645 examples/rips/rips.cpp
--- a/examples/rips/rips.cpp Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/rips/rips.cpp Wed Dec 16 15:39:06 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");
diff -r a99fdaafa31a -r d15c6d144645 examples/triangle/triangle.cpp
--- a/examples/triangle/triangle.cpp Tue Jul 28 13:27:13 2009 -0700
+++ b/examples/triangle/triangle.cpp Wed Dec 16 15:39:06 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;
}
diff -r a99fdaafa31a -r d15c6d144645 include/topology/chain.h
--- a/include/topology/chain.h Tue Jul 28 13:27:13 2009 -0700
+++ b/include/topology/chain.h Wed Dec 16 15:39:06 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
diff -r a99fdaafa31a -r d15c6d144645 include/topology/chain.hpp
--- a/include/topology/chain.hpp Tue Jul 28 13:27:13 2009 -0700
+++ b/include/topology/chain.hpp Wed Dec 16 15:39:06 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>::
diff -r a99fdaafa31a -r d15c6d144645 include/topology/complex-traits.h
--- a/include/topology/complex-traits.h Tue Jul 28 13:27:13 2009 -0700
+++ b/include/topology/complex-traits.h Wed Dec 16 15:39:06 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__
diff -r a99fdaafa31a -r d15c6d144645 include/topology/dynamic-persistence.h
--- a/include/topology/dynamic-persistence.h Tue Jul 28 13:27:13 2009 -0700
+++ b/include/topology/dynamic-persistence.h Wed Dec 16 15:39:06 2009 -0800
@@ -6,42 +6,39 @@
#include <boost/progress.hpp>
+#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;
};
/**
@@ -61,28 +58,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()
@@ -91,22 +91,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
@@ -118,65 +119,71 @@
// 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):
- bg_(bg), ccmp_(ccmp) {}
+ PairingTrailsVisitor(Order& order, ConsistencyComparison ccmp):
+ 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 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)
- 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;
};
/**
@@ -195,28 +202,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()
@@ -225,9 +235,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();
@@ -239,10 +247,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
@@ -254,17 +264,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);
@@ -272,15 +287,14 @@
struct PairingChainsVisitor: public Parent::PairVisitor
{
- // TODO: this is specialized for std::vector
- PairingChainsVisitor(OrderIndex bg, ConsistencyComparison ccmp, unsigned size):
- bg_(bg), ccmp_(ccmp), show_progress(size) {}
+ PairingChainsVisitor(Order& order, ConsistencyComparison ccmp, unsigned size):
+ order_(order), ccmp_(ccmp), show_progress(size) {}
- 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 { CountBy(cChainLength, i->chain.size()); ++show_progress; }
+ 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 { CountBy(cChainLength, i->chain.size()); ++show_progress; }
- OrderIndex bg_;
+ Order& order_;
ConsistencyComparison ccmp_;
mutable boost::progress_display
diff -r a99fdaafa31a -r d15c6d144645 include/topology/dynamic-persistence.hpp
--- a/include/topology/dynamic-persistence.hpp Tue Jul 28 13:27:13 2009 -0700
+++ b/include/topology/dynamic-persistence.hpp Wed Dec 16 15:39:06 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(), PairingTrailsVisitor(begin(), ccmp_));
+ Parent::pair_simplices(begin(), end(), PairingTrailsVisitor(order(), ccmp_));
}
-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(), PairingChainsVisitor(begin(), ccmp_, size()));
+ Parent::pair_simplices(begin(), end(), PairingChainsVisitor(order(), ccmp_, size()));
}
-
diff -r a99fdaafa31a -r d15c6d144645 include/topology/filtration.h
--- a/include/topology/filtration.h Tue Jul 28 13:27:13 2009 -0700
+++ b/include/topology/filtration.h Wed Dec 16 15:39:06 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__
diff -r a99fdaafa31a -r d15c6d144645 include/topology/filtration.hpp
--- a/include/topology/filtration.hpp Tue Jul 28 13:27:13 2009 -0700
+++ b/include/topology/filtration.hpp Wed Dec 16 15:39:06 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
diff -r a99fdaafa31a -r d15c6d144645 include/topology/lowerstarfiltration.h
--- a/include/topology/lowerstarfiltration.h Tue Jul 28 13:27:13 2009 -0700
+++ b/include/topology/lowerstarfiltration.h Wed Dec 16 15:39:06 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__
diff -r a99fdaafa31a -r d15c6d144645 include/topology/order.h
--- a/include/topology/order.h Tue Jul 28 13:27:13 2009 -0700
+++ b/include/topology/order.h Wed Dec 16 15:39:06 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__
diff -r a99fdaafa31a -r d15c6d144645 include/topology/persistence-diagram.hpp
--- a/include/topology/persistence-diagram.hpp Tue Jul 28 13:27:13 2009 -0700
+++ b/include/topology/persistence-diagram.hpp Wed Dec 16 15:39:06 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);
}
}
diff -r a99fdaafa31a -r d15c6d144645 include/topology/simplex.h
--- a/include/topology/simplex.h Tue Jul 28 13:27:13 2009 -0700
+++ b/include/topology/simplex.h Wed Dec 16 15:39:06 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>
diff -r a99fdaafa31a -r d15c6d144645 include/topology/simplex.hpp
--- a/include/topology/simplex.hpp Tue Jul 28 13:27:13 2009 -0700
+++ b/include/topology/simplex.hpp Wed Dec 16 15:39:06 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
diff -r a99fdaafa31a -r d15c6d144645 include/topology/static-persistence.h
--- a/include/topology/static-persistence.h Tue Jul 28 13:27:13 2009 -0700
+++ b/include/topology/static-persistence.h Wed Dec 16 15:39:06 2009 -0800
@@ -5,18 +5,21 @@
#include "cycles.h"
#include "filtration.h"
+#include <boost/ref.hpp>
+#include <boost/lambda/lambda.hpp>
+namespace bl = boost::lambda;
+
#include <utilities/types.h>
-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;
@@ -26,6 +29,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;
@@ -44,9 +51,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:
@@ -59,15 +67,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
@@ -75,8 +84,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
@@ -86,16 +94,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, const Visitor& visitor = Visitor());
+ void pair_simplices(iterator bg, iterator end, 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)>.
@@ -104,20 +131,27 @@
// 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 {}
+ void finished(iterator j) const {}
};
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_;
diff -r a99fdaafa31a -r d15c6d144645 include/topology/static-persistence.hpp
--- a/include/topology/static-persistence.hpp Tue Jul 28 13:27:13 2009 -0700
+++ b/include/topology/static-persistence.hpp Wed Dec 16 15:39:06 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>
#ifdef LOGGING
static rlog::RLogChannel* rlPersistence = DEF_CHANNEL("topology/persistence", rlog::Log_Debug);
@@ -14,31 +16,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, const Visitor& visitor)
+StaticPersistence<D, CT, OT, E, Cmp>::
+pair_simplices(iterator bg, iterator end, const Visitor& visitor)
{
#if LOGGING
typename ContainerTraits::OutputMap outmap(order_);
@@ -46,16 +51,16 @@
// 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());
- CountNum(cPersistencePairBoundaries, oe.cycle.size());
+ CountNum(cPersistencePairBoundaries, z.size());
Count(cPersistencePair);
while(!z.empty())
@@ -63,33 +68,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 */
diff -r a99fdaafa31a -r d15c6d144645 include/topology/vineyard.h
--- a/include/topology/vineyard.h Tue Jul 28 13:27:13 2009 -0700
+++ b/include/topology/vineyard.h Wed Dec 16 15:39:06 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,141 @@
#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;
- 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"
diff -r a99fdaafa31a -r d15c6d144645 include/topology/vineyard.hpp
--- a/include/topology/vineyard.hpp Tue Jul 28 13:27:13 2009 -0700
+++ b/include/topology/vineyard.hpp Wed Dec 16 15:39:06 2009 -0800
@@ -6,160 +6,158 @@
#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();
+ }
}
/// 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;
-}
diff -r a99fdaafa31a -r d15c6d144645 include/utilities/indirect.h
--- a/include/utilities/indirect.h Tue Jul 28 13:27:13 2009 -0700
+++ b/include/utilities/indirect.h Wed Dec 16 15:39:06 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>,
diff -r a99fdaafa31a -r d15c6d144645 include/utilities/property-maps.h
--- a/include/utilities/property-maps.h Tue Jul 28 13:27:13 2009 -0700
+++ b/include/utilities/property-maps.h Wed Dec 16 15:39:06 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):