Intermediate commit while converting to the new architecture
* Converted Filtration into StaticPersistence and DynamicPersistenceTrails
(both right now work with the underlying std::vector rather than std::list order)
* Filtration is now just an auxilliary glue (a map between Complex and Persistence)
* Whether chains are vectors or lists can be interchanged
* Added PersistenceDiagram with a simple bottleneck_distance() function
* Converted triangle, alphashapes3d, and cech-complex examples
* Lots of organizational changes
(factoring utilities out into containers.h, indirect.h, property-maps.h)
* Trying to document along the way with NaturalDocs-type comments
--- a/CMakeLists.txt Mon Sep 22 21:53:25 2008 -0700
+++ b/CMakeLists.txt Thu Dec 18 16:43:42 2008 -0800
@@ -9,7 +9,7 @@
option (use_synaps "Build examples that use SYNAPS" OFF)
# Find everything that's always required
-find_package (Boost REQUIRED)
+find_package (Boost REQUIRED COMPONENTS program_options serialization signals)
find_package (Doxygen)
if (use_dsrpdb)
find_library (dsrpdb_LIBRARY NAMES dsrpdb)
@@ -44,6 +44,7 @@
${gmp_LIBRARY}
${gmpxx_LIBRARY}
${m_LIBRARY})
+add_definitions (-DCGAL_NO_ASSERTIONS -DCGAL_NO_PRECONDITIONS)
# SYNAPS
if (use_synaps)
--- a/examples/CMakeLists.txt Mon Sep 22 21:53:25 2008 -0700
+++ b/examples/CMakeLists.txt Thu Dec 18 16:43:42 2008 -0800
@@ -1,6 +1,3 @@
-find_library (boost_program_options_LIBRARY NAME boost_program_options
- PATHS ${Boost_LIBRARY_DIR})
-
add_subdirectory (alphashapes)
add_subdirectory (ar-vineyard)
add_subdirectory (cech-complex)
--- a/examples/alphashapes/CMakeLists.txt Mon Sep 22 21:53:25 2008 -0700
+++ b/examples/alphashapes/CMakeLists.txt Thu Dec 18 16:43:42 2008 -0800
@@ -1,8 +1,13 @@
set (targets
alphashapes3d
alphashapes2d
- alpharadius)
+ alpharadius
+ compare-diagrams)
+set (libraries ${libraries}
+ ${Boost_SERIALIZATION_LIBRARY}
+ ${Boost_PROGRAM_OPTIONS_LIBRARY})
+
add_definitions (${cgal_cxxflags})
foreach (t ${targets})
add_executable (${t} ${t}.cpp)
--- a/examples/alphashapes/alphashapes3d.cpp Mon Sep 22 21:53:25 2008 -0700
+++ b/examples/alphashapes/alphashapes3d.cpp Thu Dec 18 16:43:42 2008 -0800
@@ -2,57 +2,102 @@
#include "alphashapes3d.h"
#include <topology/filtration.h>
+#include <topology/static-persistence.h>
+#include <topology/persistence-diagram.h>
#include <iostream>
+
#include <fstream>
+#include <boost/archive/binary_oarchive.hpp>
+#include <boost/serialization/map.hpp>
+
+#include <boost/program_options.hpp>
-typedef Filtration<AlphaSimplex3D> AlphaFiltration;
+typedef Filtration<AlphaSimplex3DVector> AlphaFiltration;
+typedef StaticPersistence<> Persistence;
+typedef PersistenceDiagram<> PDgm;
+
+namespace po = boost::program_options;
int main(int argc, char** argv)
{
#ifdef LOGGING
- rlog::RLogInit(argc, argv);
+ rlog::RLogInit(argc, argv);
- stdoutLog.subscribeTo( RLOG_CHANNEL("info") );
- stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
- //stdoutLog.subscribeTo( RLOG_CHANNEL("topology/filtration") );
- //stdoutLog.subscribeTo( RLOG_CHANNEL("topology/cycle") );
+ stdoutLog.subscribeTo( RLOG_CHANNEL("info") );
+ stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
+ //stdoutLog.subscribeTo( RLOG_CHANNEL("topology/persistence") );
+ //stdoutLog.subscribeTo( RLOG_CHANNEL("topology/chain") );
#endif
- SetFrequency(GetCounter("filtration/pair"), 10000);
- SetTrigger(GetCounter("filtration/pair"), GetCounter(""));
+ SetFrequency(GetCounter("persistence/pair"), 10000);
+ SetTrigger(GetCounter("persistence/pair"), GetCounter(""));
+
+ std::string infilename, outfilename;
+
+ po::options_description hidden("Hidden options");
+ hidden.add_options()
+ ("input-file", po::value<std::string>(&infilename), "Point set whose alpha shape filtration and persistence we want to compute")
+ ("output-file", po::value<std::string>(&outfilename), "Where to write the collection of persistence diagrams");
- // Read in the point set and compute its Delaunay triangulation
- std::istream& in = std::cin;
- double x,y,z;
- Delaunay Dt;
- while(in)
- {
- in >> x >> y >> z;
- Point p(x,y,z);
- Dt.insert(p);
- }
- rInfo("Delaunay triangulation computed");
-
- AlphaSimplex3DVector alpha_ordering;
- fill_alpha_order(Dt, alpha_ordering);
- rInfo("Simplices: %d", alpha_ordering.size());
+ po::positional_options_description pos;
+ pos.add("input-file", 1);
+ pos.add("output-file", 2);
+
+ po::options_description all; all.add(hidden);
+
+ po::variables_map vm;
+ po::store(po::command_line_parser(argc, argv).
+ options(all).positional(pos).run(), vm);
+ po::notify(vm);
+
+ if (!vm.count("input-file") || !vm.count("output-file"))
+ {
+ std::cout << "Usage: " << argv[0] << " input-file output-file" << std::endl;
+ std::cout << hidden << std::endl;
+ return 1;
+ }
+
- // Create the alpha-shape filtration
- AlphaFiltration af;
- for (AlphaSimplex3DVector::const_iterator cur = alpha_ordering.begin(); cur != alpha_ordering.end(); ++cur)
- af.append(*cur);
- af.fill_simplex_index_map();
- rInfo("Filled simplex-index map");
- af.pair_simplices(af.begin(), af.end(), false);
- rInfo("Simplices paired");
+ // Read in the point set and compute its Delaunay triangulation
+ std::ifstream in(infilename.c_str());
+ double x,y,z;
+ Delaunay Dt;
+ while(in)
+ {
+ in >> x >> y >> z;
+ Point p(x,y,z);
+ Dt.insert(p);
+ }
+ rInfo("Delaunay triangulation computed");
+
+ AlphaSimplex3DVector complex;
+ fill_complex(Dt, complex);
+ rInfo("Simplices: %d", complex.size());
- for (AlphaFiltration::Index i = af.begin(); i != af.end(); ++i)
- if (i->is_paired())
- {
- if (i->sign())
- std::cout << i->dimension() << " " << i->value() << " " << i->pair()->value() << std::endl;
- } //else std::cout << i->value() << std::endl;
+ // Create the alpha-shape filtration
+ AlphaFiltration af(complex.begin(), complex.end(), AlphaSimplex3D::AlphaOrder());
+ rInfo("Filtration initialized");
+
+ Persistence p(af);
+ rInfo("Persistence initializaed");
+
+ p.pair_simplices();
+ rInfo("Simplices paired");
+ 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())));
+#if 0
+ std::cout << 0 << std::endl << dgms[0] << std::endl;
+ std::cout << 1 << std::endl << dgms[1] << std::endl;
+ std::cout << 2 << std::endl << dgms[2] << std::endl;
+#endif
+
+ std::ofstream ofs(outfilename.c_str());
+ boost::archive::binary_oarchive oa(ofs);
+ oa << dgms;
}
-
--- a/examples/alphashapes/alphashapes3d.h Mon Sep 22 21:53:25 2008 -0700
+++ b/examples/alphashapes/alphashapes3d.h Thu Dec 18 16:43:42 2008 -0800
@@ -22,9 +22,6 @@
typedef Delaunay::Point Point;
typedef Delaunay::Vertex Vertex;
typedef Delaunay::Vertex_handle Vertex_handle;
-typedef Delaunay::Edge Edge;
-typedef Delaunay::Facet Facet;
-typedef Delaunay::Cell Cell;
typedef Delaunay::Cell_handle Cell_handle;
typedef K::FT RealValue;
@@ -35,13 +32,13 @@
typedef Delaunay::Facet_circulator Facet_circulator;
-class AlphaSimplex3D: public SimplexWithVertices<Vertex_handle>
+class AlphaSimplex3D: public Simplex<Vertex_handle>
{
public:
- typedef std::set<AlphaSimplex3D> SimplexSet;
- typedef SimplexWithVertices<Vertex_handle> Parent;
+ typedef Simplex<Vertex_handle> Parent;
+ typedef std::set<AlphaSimplex3D, Parent::VertexComparison> SimplexSet;
typedef Parent::VertexContainer VertexSet;
- typedef std::list<AlphaSimplex3D> Cycle;
+ typedef std::list<AlphaSimplex3D> Boundary;
public:
AlphaSimplex3D() {}
@@ -51,23 +48,31 @@
Parent(s) { attached_ = s.attached_; alpha_ = s.alpha_; }
AlphaSimplex3D(const ::Vertex& v);
- AlphaSimplex3D(const Edge& e);
- AlphaSimplex3D(const Edge& e, const SimplexSet& simplices, const Delaunay& Dt, Facet_circulator facet_bg);
+ AlphaSimplex3D(const Delaunay::Edge& e);
+ AlphaSimplex3D(const Delaunay::Edge& e, const SimplexSet& simplices, const Delaunay& Dt, Facet_circulator facet_bg);
- AlphaSimplex3D(const Facet& f);
- AlphaSimplex3D(const Facet& f, const SimplexSet& simplices, const Delaunay& Dt);
+ AlphaSimplex3D(const Delaunay::Facet& f);
+ AlphaSimplex3D(const Delaunay::Facet& f, const SimplexSet& simplices, const Delaunay& Dt);
- AlphaSimplex3D(const Cell& c);
+ AlphaSimplex3D(const Delaunay::Cell& c);
RealType value() const { return CGAL::to_double(alpha_); }
RealValue alpha() const { return alpha_; }
bool attached() const { return attached_; }
- Cycle boundary() const;
+ Boundary boundary() const;
// Ordering
struct AlphaOrder
{ bool operator()(const AlphaSimplex3D& first, const AlphaSimplex3D& second) const; };
-
+
+ struct AlphaValueEvaluator
+ {
+ typedef AlphaSimplex3D first_argument_type;
+ typedef RealType result_type;
+
+ RealType operator()(const AlphaSimplex3D& s) const { return s.value(); }
+ };
+
std::ostream& operator<<(std::ostream& out) const;
private:
@@ -76,8 +81,8 @@
};
typedef std::vector<AlphaSimplex3D> AlphaSimplex3DVector;
-void fill_alpha_order(const Delaunay& Dt,
- AlphaSimplex3DVector& alpha_order);
+void fill_complex(const Delaunay& Dt,
+ AlphaSimplex3DVector& alpha_order);
std::ostream& operator<<(std::ostream& out, const AlphaSimplex3D& s) { return s.operator<<(out); }
--- a/examples/alphashapes/alphashapes3d.hpp Mon Sep 22 21:53:25 2008 -0700
+++ b/examples/alphashapes/alphashapes3d.hpp Thu Dec 18 16:43:42 2008 -0800
@@ -7,7 +7,7 @@
}
AlphaSimplex3D::
-AlphaSimplex3D(const Edge& e)
+AlphaSimplex3D(const Delaunay::Edge& e)
{
Cell_handle c = e.first;
Parent::add(c->vertex(e.second));
@@ -15,7 +15,7 @@
}
AlphaSimplex3D::
-AlphaSimplex3D(const Edge& e, const SimplexSet& simplices, const Delaunay& Dt, Facet_circulator facet_bg)
+AlphaSimplex3D(const Delaunay::Edge& e, const SimplexSet& simplices, const Delaunay& Dt, Facet_circulator facet_bg)
{
Cell_handle c = e.first;
Parent::add(c->vertex(e.second));
@@ -26,14 +26,15 @@
SimplexSet::const_iterator cur_iter = simplices.find(AlphaSimplex3D(*cur));
RealValue min = cur_iter->alpha();
- VertexSet::const_iterator v = Parent::vertices().begin();
+ const VertexSet& vertices = static_cast<const Parent*>(this)->vertices();
+ VertexSet::const_iterator v = vertices.begin();
const Point& p1 = (*v++)->point();
const Point& p2 = (*v)->point();
attached_ = false;
if (facet_bg != 0) do
{
- VertexSet::const_iterator v = Parent::vertices().begin();
+ VertexSet::const_iterator v = vertices.begin();
int i0 = (*cur).first->index(*v++);
int i1 = (*cur).first->index(*v);
int i = 6 - i0 - i1 - (*cur).second;
@@ -56,7 +57,7 @@
}
AlphaSimplex3D::
-AlphaSimplex3D(const Facet& f)
+AlphaSimplex3D(const Delaunay::Facet& f)
{
Cell_handle c = f.first;
for (int i = 0; i < 4; ++i)
@@ -65,7 +66,7 @@
}
AlphaSimplex3D::
-AlphaSimplex3D(const Facet& f, const SimplexSet& simplices, const Delaunay& Dt)
+AlphaSimplex3D(const Delaunay::Facet& f, const SimplexSet& simplices, const Delaunay& Dt)
{
Cell_handle c = f.first;
for (int i = 0; i < 4; ++i)
@@ -75,7 +76,7 @@
Cell_handle o = c->neighbor(f.second);
int oi = o->index(c);
- VertexSet::const_iterator v = Parent::vertices().begin();
+ VertexSet::const_iterator v = static_cast<const Parent*>(this)->vertices().begin();
const Point& p1 = (*v++)->point();
const Point& p2 = (*v++)->point();
const Point& p3 = (*v)->point();
@@ -105,11 +106,11 @@
}
AlphaSimplex3D::
-AlphaSimplex3D(const Cell& c): attached_(false)
+AlphaSimplex3D(const Delaunay::Cell& c): attached_(false)
{
for (int i = 0; i < 4; ++i)
Parent::add(c.vertex(i));
- VertexSet::const_iterator v = Parent::vertices().begin();
+ VertexSet::const_iterator v = static_cast<const Parent*>(this)->vertices().begin();
Point p1 = (*v++)->point();
Point p2 = (*v++)->point();
Point p3 = (*v++)->point();
@@ -117,12 +118,12 @@
alpha_ = CGAL::squared_radius(p1, p2, p3, p4);
}
-AlphaSimplex3D::Cycle
+AlphaSimplex3D::Boundary
AlphaSimplex3D::boundary() const
{
- Cycle bdry;
- Parent::Cycle pbdry = Parent::boundary();
- for (Parent::Cycle::const_iterator cur = pbdry.begin(); cur != pbdry.end(); ++cur)
+ Boundary bdry;
+ Parent::Boundary pbdry = Parent::boundary();
+ for (Parent::Boundary::const_iterator cur = pbdry.begin(); cur != pbdry.end(); ++cur)
bdry.push_back(*cur);
return bdry;
}
@@ -150,7 +151,7 @@
}
-void fill_alpha_order(const Delaunay& Dt, AlphaSimplex3DVector& alpha_order)
+void fill_complex(const Delaunay& Dt, AlphaSimplex3DVector& alpha_order)
{
// Compute all simplices with their alpha values and attachment information
AlphaSimplex3D::SimplexSet simplices;
@@ -170,6 +171,7 @@
// 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::AlphaOrder());
+ std::sort(alpha_order.begin(), alpha_order.end(), AlphaSimplex3D::VertexComparison());
}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/compare-diagrams.cpp Thu Dec 18 16:43:42 2008 -0800
@@ -0,0 +1,56 @@
+#include <utilities/types.h>
+#include <topology/persistence-diagram.h>
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <boost/archive/binary_iarchive.hpp>
+#include <boost/serialization/map.hpp>
+
+#include <boost/program_options.hpp>
+
+
+typedef PersistenceDiagram<> PDgm;
+
+namespace po = boost::program_options;
+
+
+int main(int argc, char* argv[])
+{
+ std::string filename1, filename2;
+
+ po::options_description hidden("Hidden options");
+ hidden.add_options()
+ ("input-file1", po::value<std::string>(&filename1), "The first collection of persistence diagrams")
+ ("input-file2", po::value<std::string>(&filename2), "The second collection of persistence diagrams");
+
+ po::positional_options_description p;
+ p.add("input-file1", 1);
+ p.add("input-file2", 2);
+
+ po::options_description all; all.add(hidden);
+
+ po::variables_map vm;
+ po::store(po::command_line_parser(argc, argv).
+ options(all).positional(p).run(), vm);
+ po::notify(vm);
+
+ if (!vm.count("input-file1") || !vm.count("input-file2"))
+ {
+ std::cout << "Usage: " << argv[0] << " input-file1 input-file2" << std::endl;
+ std::cout << hidden << std::endl;
+ return 1;
+ }
+
+ std::ifstream ifs1(filename1.c_str()), ifs2(filename2.c_str());
+ boost::archive::binary_iarchive ia1(ifs1), ia2(ifs2);
+
+ std::map<Dimension, PDgm> dgms1, dgms2;
+
+ ia1 >> dgms1;
+ ia2 >> dgms2;
+
+ std::cout << "Distance between dimension 0: " << bottleneck_distance(dgms1[0], dgms2[0]) << std::endl;
+ std::cout << "Distance between dimension 1: " << bottleneck_distance(dgms1[1], dgms2[1]) << std::endl;
+ std::cout << "Distance between dimension 2: " << bottleneck_distance(dgms1[2], dgms2[2]) << std::endl;
+}
--- a/examples/ar-vineyard/CMakeLists.txt Mon Sep 22 21:53:25 2008 -0700
+++ b/examples/ar-vineyard/CMakeLists.txt Thu Dec 18 16:43:42 2008 -0800
@@ -1,15 +1,10 @@
set (targets
ar-vineyard)
-find_library (boost_signals_LIBRARY NAMES boost_signals
- PATHS ${Boost_LIBRARY_DIR})
-find_library (boost_program_options_LIBRARY NAME boost_program_options
- PATHS ${Boost_LIBRARY_DIR})
-
add_definitions (${cgal_cxxflags})
foreach (t ${targets})
add_executable (${t} ${t}.cpp ${external_sources})
target_link_libraries (${t} ${libraries} ${cgal_libraries}
- ${boost_signals_LIBRARY}
- ${boost_program_options_LIBRARY})
+ ${Boost_SIGNALS_LIBRARY}
+ ${Boost_PROGRAM_OPTIONS_LIBRARY})
endforeach (t)
--- a/examples/cech-complex/cech-complex.cpp Mon Sep 22 21:53:25 2008 -0700
+++ b/examples/cech-complex/cech-complex.cpp Thu Dec 18 16:43:42 2008 -0800
@@ -2,6 +2,9 @@
#include <topology/simplex.h>
#include <topology/filtration.h>
+#include <topology/static-persistence.h>
+//#include <topology/dynamic-persistence.h>
+#include <topology/persistence-diagram.h>
#include "Miniball_dynamic_d.h"
#include <algorithm>
@@ -9,112 +12,124 @@
#include <fstream>
-typedef std::vector<Point> PointContainer;
-typedef unsigned int PointIndex;
-typedef SimplexWithValue<PointIndex> Simplex;
-typedef std::vector<Simplex> SimplexVector;
-typedef Filtration<Simplex> CechFiltration;
-typedef DimensionValueComparison<Simplex> DimValComparison;
+typedef std::vector<Point> PointContainer;
+typedef unsigned int PointIndex;
+typedef Simplex<PointIndex, double> Smplx;
+typedef std::vector<Smplx> SimplexVector;
+typedef Filtration<SimplexVector> CechFiltration;
+typedef StaticPersistence<> Persistence;
+//typedef DynamicPersistenceTrails<> Persistence;
+typedef PersistenceDiagram<> PDgm;
int choose(int n, int k)
{
- if (k > n/2) k = n-k;
- int numerator = 1, denominator = 1;
- for (int i = 0; i < k; ++i)
- { numerator *= (n-i); denominator *= (i+1); }
- return numerator/denominator;
+ if (k > n/2) k = n-k;
+ int numerator = 1, denominator = 1;
+ for (int i = 0; i < k; ++i)
+ { numerator *= (n-i); denominator *= (i+1); }
+ return numerator/denominator;
}
void add_simplices(SimplexVector& sv, int d, const PointContainer& points)
{
- PointIndex indices[d+1];
- for (int i = 0; i < d+1; ++i)
- indices[i] = d - i;
+ PointIndex indices[d+1];
+ for (int i = 0; i < d+1; ++i)
+ indices[i] = d - i;
- while(indices[d] < points.size() - d)
- {
- // Add simplex
- Miniball mb(points[indices[0]].dim());
- Simplex s;
- for (int i = 0; i < d+1; ++i)
- {
- s.add(indices[i]);
- mb.check_in(points[indices[i]]);
- }
- mb.build();
- s.set_value(mb.squared_radius());
- sv.push_back(s);
+ while(indices[d] < points.size() - d)
+ {
+ // Add simplex
+ Miniball mb(points[indices[0]].dim());
+ Smplx s;
+ for (int i = 0; i < d+1; ++i)
+ {
+ s.add(indices[i]);
+ mb.check_in(points[indices[i]]);
+ }
+ mb.build();
+ s.data() = mb.squared_radius();
+ sv.push_back(s);
-
- // Advance indices
- for (int i = 0; i < d+1; ++i)
- {
- ++indices[i];
- if (indices[i] < points.size() - i)
- {
- for (int j = i-1; j >= 0; --j)
- indices[j] = indices[j+1] + 1;
- break;
- }
- }
- }
+
+ // Advance indices
+ for (int i = 0; i < d+1; ++i)
+ {
+ ++indices[i];
+ if (indices[i] < points.size() - i)
+ {
+ for (int j = i-1; j >= 0; --j)
+ indices[j] = indices[j+1] + 1;
+ break;
+ }
+ }
+ }
}
int main(int argc, char** argv)
{
- SetFrequency(GetCounter("filtration/pair"), 10000);
- SetTrigger(GetCounter("filtration/pair"), GetCounter(""));
+#ifdef LOGGING
+ rlog::RLogInit(argc, argv);
+
+ stdoutLog.subscribeTo( RLOG_CHANNEL("info") );
+ stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
+#endif
+
+ SetFrequency(GetCounter("persistence/pair"), 10000);
+ SetTrigger(GetCounter("persistence/pair"), GetCounter(""));
- // Read in the point set and compute its Delaunay triangulation
- std::istream& in = std::cin;
- int ambient_d, homology_d;
- in >> ambient_d >> homology_d;
-
- std::cout << "Ambient dimension: " << ambient_d << std::endl;
- std::cout << "Will compute PD up to dimension: " << homology_d << std::endl;
-
- // Read points
- PointContainer points;
- while(in)
- {
- Point p(ambient_d);
- for (int i = 0; i < ambient_d; ++i)
- in >> p[i];
- points.push_back(p);
- }
- std::cout << "Points read: " << points.size() << std::endl;
+ // Read in the point set and compute its Delaunay triangulation
+ std::istream& in = std::cin;
+ int ambient_d, homology_d;
+ in >> ambient_d >> homology_d;
+
+ rInfo("Ambient dimension: %d", ambient_d);
+ rInfo("Will compute PD up to dimension: %d", homology_d);
+
+ // Read points
+ PointContainer points;
+ while(in)
+ {
+ Point p(ambient_d);
+ for (int i = 0; i < ambient_d; ++i)
+ in >> p[i];
+ points.push_back(p);
+ }
+ rInfo("Points read: %d", points.size());
- // Compute Cech values
- CechFiltration cf;
- { // resource acquisition is initialization
- 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);
- std::cout << "Reserved SimplexVector of size: " << num_simplices << std::endl;
+ // 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);
+
+ 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());
+
+ // Set up the filtration
+ CechFiltration cf(sv.begin(), sv.end(), Smplx::DataDimensionComparison());
+ rInfo("Filtration initialized");
- for (int i = 0; i <= homology_d + 1; ++i)
- add_simplices(sv, i, points);
- std::cout << "Size of SimplexVector: " << sv.size() << std::endl;
-
- std::sort(sv.begin(), sv.end(), DimValComparison());
-
- for (SimplexVector::const_iterator cur = sv.begin(); cur != sv.end(); ++cur)
- cf.append(*cur);
- }
+ // Compute persistence
+ Persistence p(cf);
+ rInfo("Persistence initialized");
+ p.pair_simplices();
+ rInfo("Simplices paired");
- // Compute persistence
- cf.fill_simplex_index_map();
- cf.pair_simplices(cf.begin(), cf.end(), false);
- std::cout << "Simplices paired" << std::endl;
+ 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())));
- for (CechFiltration::Index i = cf.begin(); i != cf.end(); ++i)
- if (i->is_paired())
- {
- if (i->sign())
- std::cout << i->dimension() << " " << i->get_value() << " " << i->pair()->get_value() << std::endl;
- } //else std::cout << i->value() << std::endl;
-
+ for (int i = 0; i <= homology_d; ++i)
+ {
+ std::cout << i << std::endl << dgms[i] << std::endl;
+ }
}
--- a/examples/fitness/CMakeLists.txt Mon Sep 22 21:53:25 2008 -0700
+++ b/examples/fitness/CMakeLists.txt Thu Dec 18 16:43:42 2008 -0800
@@ -5,5 +5,5 @@
foreach (t ${targets})
add_executable (${t} ${t}.cpp)
- target_link_libraries (${t} ${libraries} ${boost_program_options_LIBRARY})
+ target_link_libraries (${t} ${libraries} ${Boost_PROGRAM_OPTIONS_LIBRARY})
endforeach (t ${targets})
--- a/examples/poincare/CMakeLists.txt Mon Sep 22 21:53:25 2008 -0700
+++ b/examples/poincare/CMakeLists.txt Thu Dec 18 16:43:42 2008 -0800
@@ -3,5 +3,5 @@
foreach (t ${targets})
add_executable (${t} ${t}.cpp)
- target_link_libraries (${t} ${libraries} ${boost_program_options_LIBRARY})
+ target_link_libraries (${t} ${libraries} ${Boost_PROGRAM_OPTIONS_LIBRARY})
endforeach (t ${targets})
--- a/examples/triangle/CMakeLists.txt Mon Sep 22 21:53:25 2008 -0700
+++ b/examples/triangle/CMakeLists.txt Thu Dec 18 16:43:42 2008 -0800
@@ -1,7 +1,9 @@
-set (targets
- triangle)
-
-foreach (t ${targets})
- add_executable (${t} ${t}.cpp)
- target_link_libraries (${t} ${libraries})
-endforeach (t ${targets})
+set (targets
+ triangle)
+
+set (libraries ${libraries} ${Boost_SERIALIZATION_LIBRARY})
+
+foreach (t ${targets})
+ add_executable (${t} ${t}.cpp)
+ target_link_libraries (${t} ${libraries})
+endforeach (t ${targets})
--- a/examples/triangle/triangle.cpp Mon Sep 22 21:53:25 2008 -0700
+++ b/examples/triangle/triangle.cpp Thu Dec 18 16:43:42 2008 -0800
@@ -1,70 +1,123 @@
#include <utilities/log.h>
+#include "topology/simplex.h"
#include "topology/filtration.h"
-#include "topology/simplex.h"
+//#include "topology/static-persistence.h"
+#include "topology/dynamic-persistence.h"
+#include "topology/persistence-diagram.h"
+#include <utilities/indirect.h>
+
#include <vector>
+#include <map>
#include <iostream>
-typedef SimplexWithValue<int> Simplex;
-typedef Filtration<Simplex> TriangleFiltration;
+
+#if 1
+#include <fstream>
+#include <boost/archive/text_oarchive.hpp>
+#include <boost/archive/text_iarchive.hpp>
+#include <boost/serialization/vector.hpp>
+#endif
-void fillTriangleSimplices(TriangleFiltration& f)
+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;
+
+void fillTriangleSimplices(Complex& c)
{
- typedef std::vector<int> VertexVector;
- VertexVector vertices(4);
- vertices[0] = 0; vertices[1] = 1; vertices[2] = 2;
- vertices[3] = 0;
-
- VertexVector::const_iterator bg = vertices.begin();
- VertexVector::const_iterator end = vertices.end();
- f.append(Simplex(bg, bg + 1, 0)); // 0 = A
- f.append(Simplex(bg + 1, bg + 2, 1)); // 1 = B
- f.append(Simplex(bg + 2, bg + 3, 2)); // 2 = C
- f.append(Simplex(bg, bg + 2, 2.5)); // AB
- f.append(Simplex(bg + 1, bg + 3, 2.9)); // BC
- f.append(Simplex(bg + 2, end, 3.5)); // CA
- f.append(Simplex(bg, bg + 3, 5)); // ABC
+ typedef std::vector<Vertex> VertexVector;
+ VertexVector vertices(4);
+ vertices[0] = 0; vertices[1] = 1; vertices[2] = 2;
+ vertices[3] = 0;
+
+ VertexVector::const_iterator bg = vertices.begin();
+ VertexVector::const_iterator end = vertices.end();
+ c.push_back(Smplx(bg, bg + 1, 0)); // 0 = A
+ c.push_back(Smplx(bg + 1, bg + 2, 1)); // 1 = B
+ c.push_back(Smplx(bg + 2, bg + 3, 2)); // 2 = C
+ c.push_back(Smplx(bg, bg + 2, 2.5)); // AB
+ 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)
{
#ifdef LOGGING
- rlog::RLogInit(argc, argv);
+ rlog::RLogInit(argc, argv);
+
+ stdoutLog.subscribeTo(RLOG_CHANNEL("topology/persistence"));
+ //stdoutLog.subscribeTo(RLOG_CHANNEL("topology/chain"));
+ //stdoutLog.subscribeTo(RLOG_CHANNEL("topology/vineyard"));
+#endif
+
+ Complex c;
+ fillTriangleSimplices(c);
+ std::cout << "Simplices filled" << std::endl;
+ for (Complex::const_iterator cur = c.begin(); cur != c.end(); ++cur)
+ std::cout << " " << *cur << std::endl;
- stdoutLog.subscribeTo(RLOG_CHANNEL("topology/filtration"));
- //stdoutLog.subscribeTo(RLOG_CHANNEL("topology/cycle"));
- //stdoutLog.subscribeTo(RLOG_CHANNEL("topology/vineyard"));
+#if 1 // testing serialization of Complex (really Simplex)
+ {
+ std::ofstream ofs("complex");
+ boost::archive::text_oarchive oa(ofs);
+ oa << c;
+ c.clear();
+ }
+
+ {
+ std::ifstream ifs("complex");
+ boost::archive::text_iarchive ia(ifs);
+ ia >> c;
+ }
#endif
- Evaluator<Simplex> e;
- TriangleFiltration::Vineyard v(&e);
- TriangleFiltration tf(&v);
- fillTriangleSimplices(tf);
-
- tf.fill_simplex_index_map();
- tf.pair_simplices();
- v.start_vines(tf.begin(), tf.end());
-
- std::cout << "Filtration size: " << tf.size() << std::endl;
- std::cout << tf << std::endl;
+ Fltr f(c.begin(), c.end(), Smplx::DataDimensionComparison());
+ std::cout << "Filtration initialized" << std::endl;
+ std::cout << f << std::endl;
+
+ Persistence p(f);
+ std::cout << "Persistence initialized" << std::endl;
+
+ p.pair_simplices();
+ std::cout << "Simplices paired" << std::endl;
+
+ std::map<Dimension, PDgm> dgms;
+ init_diagrams(dgms, p.begin(), p.end(),
+ evaluate_through_map(OffsetMap<Persistence::OrderIndex, Fltr::Index>(p.begin(), f.begin()),
+ evaluate_through_filtration(f, Smplx::DataEvaluator())),
+ evaluate_through_map(OffsetMap<Persistence::OrderIndex, Fltr::Index>(p.begin(), f.begin()),
+ evaluate_through_filtration(f, Smplx::DimensionExtractor())));
+
+ std::cout << 0 << std::endl << dgms[0] << std::endl;
+ std::cout << 1 << std::endl << dgms[1] << std::endl;
-#if 1
- Simplex A; A.add(0);
- std::cout << A << std::endl;
- std::cout << *tf.get_index(A) << std::endl;
+ // Transpositions
+ p.transpose(p.begin()); // transposition case 1.2 special
+
+#if 0
+ 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 << tf;
- Simplex BC; BC.add(1); BC.add(2);
- Simplex 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 << AB << std::endl;
- std::cout << *tf.get_index(AB) << std::endl;
- tf.transpose(tf.get_index(AB));
- std::cout << tf;
+ 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 << AB << std::endl;
+ std::cout << *tf.get_index(AB) << std::endl;
+ tf.transpose(tf.get_index(AB));
+ std::cout << tf;
#endif
}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/chain.h Thu Dec 18 16:43:42 2008 -0800
@@ -0,0 +1,120 @@
+/**
+ * Author: Dmitriy Morozov
+ * Department of Computer Science, Duke University, 2005-2008
+ */
+
+#ifndef __CHAIN_H__
+#define __CHAIN_H__
+
+//#include "utilities/types.h"
+#include <boost/serialization/access.hpp>
+#include <boost/serialization/serialization.hpp>
+#include <boost/iterator/iterator_traits.hpp>
+#include <boost/optional.hpp>
+#include <iterator>
+#include <string>
+
+#include <utilities/containers.h>
+
+/**
+ * Class: ChainWrapper
+ * Class storing a chain of simplices. Stores items in the order defined by ConsistencyCmp.
+ * The actual order of the elements is defined by OrderCmp. Instances of those
+ * classes are not stored in Chain for efficiency, and are passed as arguments to those methods
+ * that require them.
+ *
+ * \ingroup topology
+ *
+ * TODO: add field (defaulting to Z_2)
+ */
+template <class Container_>
+class ChainWrapper: public Container_,
+ public SizeStorage<Container_>
+{
+ public:
+ /// \name Template Parameters
+ /// @{
+ typedef Container_ Container;
+ typedef typename boost::iterator_value<typename Container::iterator>::type Item;
+ /// @}
+
+ typedef ChainWrapper Self;
+ typedef Container ChainRepresentation;
+
+ /// \name Accessor typedefs
+ /// @{
+ typedef typename ChainRepresentation::iterator iterator;
+ typedef typename ChainRepresentation::const_iterator const_iterator;
+ typedef typename ChainRepresentation::const_reference const_reference;
+ typedef typename ChainRepresentation::reference reference;
+ typedef typename ChainRepresentation::pointer pointer;
+ typedef typename ChainRepresentation::const_pointer const_pointer;
+ typedef Item value_type;
+ /// @}
+
+ public:
+ ChainWrapper();
+ ChainWrapper(const ChainWrapper& c);
+
+ /// \name Whole Chain operations
+ /// @{
+ /** Add c to *this assuming both Chains are sorted in increasing order according to cmp. */
+ template<class ConsistencyComparison>
+ Self& add(const Self& c, const ConsistencyComparison& cmp);
+
+ void swap(ChainWrapper& c); ///< Swaps the contents of c and *this (like STL's swap destroys c)
+
+ template<class ConsistencyComparison>
+ void sort(const ConsistencyComparison& cmp); ///< Sort elements to enforce consistency
+
+ size_t size() const { return SizeStorage<Container>::size(*this); }
+
+ using ChainRepresentation::empty;
+ using ChainRepresentation::clear;
+ /// @}
+
+ /// \name Modifiers
+ /// @{
+ using ChainRepresentation::erase;
+
+ template<class ConsistencyComparison>
+ void append(const_reference x, const ConsistencyComparison& cmp);
+ /// @}
+
+ /// \name Accessors
+ /// @{
+ using ChainRepresentation::begin;
+ using ChainRepresentation::end;
+
+ template<class OrderComparison>
+ 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
+ /// @}
+
+ /// \name Debugging
+ /// @{
+ template<class OrderComparison>
+ const_reference get_first(const OrderComparison& cmp) const; ///< First element in cmp order
+
+ template<class OutputMap>
+ std::string tostring(const OutputMap& outmap = OutputMap()) const;
+ /// @}
+
+ private:
+ using ChainRepresentation::push_back;
+ using ChainRepresentation::insert;
+
+ private:
+ // Serialization
+ typedef Container Parent;
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive& ar, boost::serialization::version_type );
+};
+
+#include "chain.hpp"
+
+#endif // __CHAIN_H__
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/chain.hpp Thu Dec 18 16:43:42 2008 -0800
@@ -0,0 +1,222 @@
+#include <algorithm>
+#include <vector>
+#include "utilities/containers.h"
+
+#include <boost/serialization/split_member.hpp>
+#include <boost/serialization/nvp.hpp>
+#include <boost/serialization/binary_object.hpp>
+#include <boost/utility.hpp>
+
+#include "utilities/log.h"
+#include "utilities/counter.h"
+
+using boost::serialization::make_nvp;
+using boost::serialization::make_binary_object;
+
+#ifdef LOGGING
+static rlog::RLogChannel* rlChain = DEF_CHANNEL( "topology/chain", rlog::Log_Debug);
+#endif // LOGGING
+
+#ifdef COUNTERS
+static Counter* cChainAddBasic = GetCounter("chain/add/basic");
+static Counter* cChainAddComparison = GetCounter("chain/add/comparison");
+#endif // COUNTERS
+
+template<class C>
+ChainWrapper<C>::
+ChainWrapper()
+{}
+
+template<class C>
+ChainWrapper<C>::
+ChainWrapper(const ChainWrapper& c): ChainRepresentation(c)
+{}
+
+template<class C>
+template<class ConsistencyCmp>
+void
+ChainWrapper<C>::
+append(const_reference x, const ConsistencyCmp& cmp)
+{
+ SizeStorage<Container>::increment();
+
+ // First try the special cases that x goes at the end
+ const_reference last = ChainRepresentation::back();
+ if (empty() || cmp(last, x))
+ {
+ push_back(x);
+ return;
+ }
+
+ for (iterator cur = begin(); cur != end(); ++cur)
+ if (cmp(x, *cur))
+ {
+ insert(cur, x);
+ return;
+ }
+}
+
+template<class C>
+template<class OrderComparison>
+typename ChainWrapper<C>::const_reference
+ChainWrapper<C>::
+top(const OrderComparison& cmp) const
+{
+ AssertMsg(!empty(), "Chain must not be empty for low()");
+ const_iterator min = begin();
+ for (const_iterator cur = ++begin(); cur != end(); ++cur)
+ if (cmp(*cur, *min))
+ min = cur;
+ return *min;
+}
+
+template<class C>
+void
+ChainWrapper<C>::
+swap(ChainWrapper& c)
+{
+ ChainRepresentation::swap(c);
+ SizeStorage<Container>::swap(c);
+}
+
+template<class C>
+template<class ConsistencyComparison>
+void
+ChainWrapper<C>::
+sort(const ConsistencyComparison& cmp)
+{
+ ContainerTraits<C,ConsistencyComparison>::sort(*this, cmp);
+}
+
+template<class C>
+boost::optional<typename ChainWrapper<C>::const_iterator>
+ChainWrapper<C>::
+contains(const_reference x) const
+{
+#if 0
+ for (const_iterator cur = begin(); cur != end(); ++cur)
+ if (!cmp(*cur, x) && !cmp(x, *cur))
+ return make_optional(cur);
+
+ return make_optional(false);
+#endif
+ const_iterator res = std::find(begin(), end(), x);
+ return make_optional(res != end(), res);
+}
+
+template<class C>
+boost::optional<typename ChainWrapper<C>::iterator>
+ChainWrapper<C>::
+contains(reference x)
+{
+ iterator res = std::find(begin(), end(), x);
+ return boost::make_optional(res != end(), res);
+}
+
+template<class C>
+template<class OutputMap>
+std::string
+ChainWrapper<C>::
+tostring(const OutputMap& outmap) const
+{
+ std::string str;
+ for (const_iterator cur = begin(); cur != end(); ++cur)
+ {
+ if (cur != begin()) str += ", ";
+ str += outmap(*cur);
+ }
+ // str += "(last: " + *last + ")"; // For debugging only
+ return str;
+}
+
+template<class C>
+template<class ConsistencyCmp>
+typename ChainWrapper<C>::Self&
+ChainWrapper<C>::
+add(const Self& c, const ConsistencyCmp& cmp)
+{
+ // TODO: tmp-based addition is necessary and useful for Containers that are vectors,
+ // however, I believe it creates costly overhead for Containers that are lists.
+ // Need to put some thought into this.
+ ChainRepresentation tmp;
+ SizeStorage<Container> size;
+
+ // FIXME: need to do something about the output
+ rLog(rlChain, "Adding chains"); //: ": %s + %s", tostring(*this).c_str(), tostring(c).c_str());
+
+ iterator cur1 = begin();
+ const_iterator cur2 = c.begin();
+
+ while (cur2 != c.end())
+ {
+ if (cur1 == end())
+ {
+ while (cur2 != c.end())
+ {
+ tmp.push_back(*cur2++);
+ size.increment();
+ Count(cChainAddBasic);
+ }
+ rLog(rlChain, "After addition"); //: %s", tostring(*this).c_str());
+
+ static_cast<ChainRepresentation*>(this)->swap(tmp);
+ static_cast<SizeStorage<Container>*>(this)->swap(size);
+ return *this;
+ }
+
+ // mod 2
+ int res = cmp.compare(*cur1, *cur2);
+ Count(cChainAddComparison);
+ rLog(rlChain, "Comparison result: %i", res);
+ if (res == 0) // *cur1 == *cur2
+ {
+ rLog(rlChain, "Equality");
+ //cur1 = erase(cur1); // erase cur1 --- as a result cur1 will be pointing at old_cur1++
+ ++cur1;
+ ++cur2;
+ } else if (res < 0) // *cur1 < *cur2
+ {
+ rLog(rlChain, "Less than");
+ //cur1++;
+ tmp.push_back(*cur1++);
+ size.increment();
+ } else if (res > 0) // *cur1 > *cur2
+ {
+ rLog(rlChain, "Greater than");
+ //insert(cur1, *cur2);
+ //++cur2;
+ tmp.push_back(*cur2++);
+ size.increment();
+ }
+ Count(cChainAddBasic);
+ }
+ while (cur1 != end())
+ {
+ tmp.push_back(*cur1++);
+ size.increment();
+ Count(cChainAddBasic);
+ }
+
+ rLog(rlChain, "After addition"); //: %s", tostring(*this).c_str());
+
+ static_cast<ChainRepresentation*>(this)->swap(tmp);
+ static_cast<SizeStorage<Container>*>(this)->swap(size);
+ return *this;
+}
+
+template<class C>
+template<class OrderComparison>
+typename ChainWrapper<C>::const_reference
+ChainWrapper<C>::
+get_first(const OrderComparison& cmp) const
+{ return top(cmp); }
+
+
+template<class C>
+template<class Archive>
+void
+ChainWrapper<C>::
+serialize(Archive& ar, boost::serialization::version_type )
+{
+ ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP(Parent);
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/complex-traits.h Thu Dec 18 16:43:42 2008 -0800
@@ -0,0 +1,71 @@
+#ifndef __COMPLEX_TRAITS_H__
+#define __COMPLEX_TRAITS_H__
+
+#include <vector>
+#include <map>
+#include "utilities/property-maps.h"
+
+// Class: ComplexTraits
+// Class that extracts various properties of a complex.
+//
+// Parameter:
+// Complex - the type whose traits we are extracting
+template<class Complex_>
+struct ComplexTraits
+{
+ /*
+ * Typedefs: Associated types
+ * Types extractable by the traits class.
+ *
+ * Complex - the type of the complex
+ * Simplex - the type of the simplex stored in the complex
+ * Index - the type used to refer to the elements of the complex
+ * IndexComparison - the type used to compare indices
+ */
+ typedef Complex_ Complex;
+ typedef typename Complex::Simplex Simplex;
+ typedef typename Complex::Index Index;
+ typedef std::less<Index> IndexComparison;
+
+ /* Typedefs: Maps
+ * Maps to go from Indicies to Simplices and back.
+ *
+ * SimplexIndexMap - the map from Simplex to the Index (could be implemented via a binary search, or a map)
+ */
+ typedef AssociativeMap<std::map<Simplex, Index,
+ typename Simplex::VertexComparison> > SimplexIndexMap;
+
+ // Function: simplex_index_map()
+ // Initializes an SimplexIndexMap given a Complex.
+ static SimplexIndexMap
+ simplex_index_map(const Complex& complex) { return SimplexIndexMap(complex.begin(),
+ complex.end()); }
+
+ static Index begin(const Complex& complex) { return complex.begin(); }
+ static Index end(const Complex& complex) { return complex.end(); }
+};
+
+
+// Specializations
+template<class Simplex_>
+struct ComplexTraits<std::vector<Simplex_> >
+{
+ typedef std::vector<Simplex_> Complex;
+ typedef Simplex_ Simplex;
+ typedef typename Complex::const_iterator Index;
+ typedef std::less<Index> IndexComparison;
+
+ typedef BinarySearchMap<Simplex, Index,
+ typename Simplex::VertexComparison> SimplexIndexMap;
+
+ static SimplexIndexMap
+ simplex_index_map(const Complex& complex) { return SimplexIndexMap(complex.begin(),
+ complex.end()); }
+ static SimplexIndexMap
+ simplex_index_map(Index bg, Index end) { return SimplexIndexMap(bg, end); }
+
+ static Index begin(const Complex& complex) { return complex.begin(); }
+ static Index end(const Complex& complex) { return complex.end(); }
+};
+
+#endif // __COMPLEX_TRAITS_H__
--- a/include/topology/cycle.h Mon Sep 22 21:53:25 2008 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,114 +0,0 @@
-/**
- * Author: Dmitriy Morozov
- * Department of Computer Science, Duke University, 2005-2006
- */
-
-#ifndef __CYCLE_H__
-#define __CYCLE_H__
-
-#include "utilities/types.h"
-#include "utilities/circular_list.h"
-#include <list>
-#include <boost/serialization/access.hpp>
-
-/**
- * Class storing a cycle of simplices. Stores items in the order defined by ConsistencyCmp.
- * The actual order of the elements is defined by OrderCmp. Instances of those
- * classes are not stored in Cycle for efficiency, and are passed as arguments to those methods
- * that require them.
- *
- * \ingroup topology
- */
-template <class Itm, class OrderCmp, class ConsistencyCmp = OrderCmp>
-class Cycle: public List<Itm>
-{
- public:
- /// \name Template Parameters
- /// @{
- typedef Itm Item;
- typedef OrderCmp OrderComparison;
- typedef ConsistencyCmp ConsistencyComparison;
- /// @}
-
- typedef Cycle<Item, OrderComparison, ConsistencyCmp> Self;
- typedef List<Item> CycleRepresentation;
-
- /// \name Accessor typedefs
- /// @{
- typedef typename CycleRepresentation::iterator iterator;
- typedef typename CycleRepresentation::const_iterator const_iterator;
- typedef typename CycleRepresentation::const_reference const_reference;
- typedef typename CycleRepresentation::reference reference;
- typedef typename CycleRepresentation::pointer pointer;
- typedef typename CycleRepresentation::const_pointer const_pointer;
- /// @}
-
- public:
- Cycle();
- Cycle(const Cycle& c);
-
- /// \name Whole Cycle operations
- /// @{
- /** Add c to *this assuming both Cycles are sorted in increasing order according to cmp. */
- Self& add(const Self& c, const ConsistencyComparison& cmp);
- void swap(Cycle& c); ///< Swaps the contents of c and *this (like STL's swap destroys c)
- //void insertion_sort(const Comparison& cmp); ///< Sort list[i] using insertion sort
- void sort(const ConsistencyComparison& cmp); ///< Sort elements to enforce consistency
- using CycleRepresentation::empty;
- using CycleRepresentation::clear;
- using CycleRepresentation::size;
- /// @}
-
- /// \name Modifiers
- /// @{
- using CycleRepresentation::erase;
- void append(const_reference x, const ConsistencyComparison& cmp);
- /// @}
-
- /// \name Accessors
- /// @{
- using CycleRepresentation::begin;
- using CycleRepresentation::end;
- const_reference top(const OrderComparison& cmp) const; ///< First element in cmp order
- iterator get_second(const OrderComparison& cmp) const; ///< Second element in cmp order
- /// @}
-
- /// \name Block access optimizations
- // Between operations used in optimization of transpositions for regular vertices. Maybe we don't need these? TODO
- /// @{
- /** Return first after i, but before or equal j; return i if no such element found */
- const_reference first_between(const_reference i, const_reference j, const OrderComparison& cmp);
- /// Add lists and remove all elements after i and before or equal to j
- const_reference add_and_first_between(const Self& c, const ConsistencyComparison& consistency_cmp,
- const_reference i, const_reference j, const OrderComparison& order_cmp);
- /// Erase everything after i, but before or equal to j
- void erase_between(const_reference i, const_reference j, const OrderComparison& cmp);
- /// @}
-
- /// \name Debugging
- /// @{
- const_reference get_first(const OrderComparison& cmp) const; ///< First element in cmp order
- std::ostream& operator<<(std::ostream& out) const;
- /// @}
-
- private:
- typedef std::list<Item> TemporaryCycleRepresenation;
-
- using CycleRepresentation::push_back;
- using CycleRepresentation::insert;
-
- private:
- size_t sz;
-
- private:
- // Serialization
- typedef List<Item> Parent;
- friend class boost::serialization::access;
- template<class Archive>
- void serialize(Archive& ar, version_type );
-};
-
-#include "cycle.hpp"
-
-#endif // __CYCLE_H__
-
--- a/include/topology/cycle.hpp Mon Sep 22 21:53:25 2008 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,294 +0,0 @@
-#include <algorithm>
-#include <vector>
-
-#include <boost/serialization/split_member.hpp>
-#include <boost/serialization/nvp.hpp>
-#include <boost/serialization/binary_object.hpp>
-#include <boost/utility.hpp>
-
-#include "utilities/log.h"
-#include "utilities/counter.h"
-
-using boost::serialization::make_nvp;
-using boost::serialization::make_binary_object;
-
-#ifdef LOGGING
-static rlog::RLogChannel* rlCycle = DEF_CHANNEL( "topology/cycle", rlog::Log_Debug);
-#endif // LOGGING
-
-#ifdef COUNTERS
-static Counter* cCycleAddBasic = GetCounter("cycle/add/basic");
-static Counter* cCycleAddComparison = GetCounter("cycle/add/comparison");
-#endif // COUNTERS
-
-template<class I, class OrderCmp, class ConsistencyCmp>
-Cycle<I,OrderCmp,ConsistencyCmp>::
-Cycle()
-{}
-
-template<class I, class OrderCmp, class ConsistencyCmp>
-Cycle<I,OrderCmp,ConsistencyCmp>::
-Cycle(const Cycle& c): CycleRepresentation(c)
-{}
-
-template<class I, class OrderCmp, class ConsistencyCmp>
-void
-Cycle<I,OrderCmp,ConsistencyCmp>::
-append(const_reference x, const ConsistencyCmp& cmp)
-{
- // First try the special cases that x goes at the end
- const_reference last = CycleRepresentation::back();
- if (empty() || cmp(last, x))
- {
- push_back(x);
- return;
- }
-
- for (iterator cur = begin(); cur != end(); ++cur)
- if (cmp(x, *cur))
- {
- insert(cur, x);
- return;
- }
-}
-
-template<class I, class OrderCmp, class ConsistencyCmp>
-typename Cycle<I,OrderCmp,ConsistencyCmp>::const_reference
-Cycle<I,OrderCmp,ConsistencyCmp>::
-top(const OrderComparison& cmp) const
-{
- AssertMsg(!empty(), "Cycle must not be empty for top()");
- const_iterator min = begin();
- for (const_iterator cur = ++begin(); cur != end(); ++cur)
- if (cmp(*cur, *min))
- min = cur;
- return *min;
-}
-
-template<class I, class OrderCmp, class ConsistencyCmp>
-void
-Cycle<I,OrderCmp,ConsistencyCmp>::
-swap(Cycle& c)
-{
- CycleRepresentation::swap(c);
-}
-
-template<class I, class OrderCmp, class ConsistencyCmp>
-void
-Cycle<I,OrderCmp,ConsistencyCmp>::
-sort(const ConsistencyComparison& cmp)
-{
- std::vector<Item> tmp(begin(), end());
- std::sort(tmp.begin(), tmp.end(), cmp);
- std::copy(tmp.begin(), tmp.end(), begin());
-}
-
-template<class I, class OrderCmp, class ConsistencyCmp>
-typename Cycle<I,OrderCmp,ConsistencyCmp>::iterator
-Cycle<I,OrderCmp,ConsistencyCmp>::
-get_second(const OrderComparison& cmp) const
-{
- AssertMsg(!empty(), "Cycle must not be empty for get_second()");
- if (size() < 2) return begin(); // Indicates that there is no second.
-
- rLog(rlCycle, "Looking for second");
- AssertMsg(size() >= 2, "Cycle must have at least two elements for get_second()");
- iterator min = begin();
- iterator second = ++begin();
- if (cmp(*second, *min))
- std::swap(min, second);
- for (iterator cur = boost::next(begin(),2); cur != end(); ++cur)
- {
- if (cmp(*cur, *min))
- {
- second = min;
- min = cur;
- } else if (cmp(*cur, *second))
- {
- second = cur;
- }
- }
-
- rLog(rlCycle, "Looked up: %s", tostring(**second).c_str());
- return second;
-}
-
-template<typename I, class OrderCmp, class ConsistencyCmp>
-typename Cycle<I,OrderCmp,ConsistencyCmp>::const_reference
-Cycle<I,OrderCmp,ConsistencyCmp>::
-first_between(const_reference i, const_reference j, const OrderComparison& cmp)
-{
- // Find the first element in ConsistencyComparison order (> i and <= j)
- const_pointer first = &i;
- iterator cur = begin();
- for (; cur != end(); ++cur)
- {
- if ((*cur == j) || (cmp(*cur, j) && cmp(i, *cur)))
- {
- first = &(*cur);
- break;
- }
- }
-
- // If no such element found, then we are done
- if (cur == end())
- return i;
-
- // Find first element in OrderComparison order (> i and <= j)
- for (++cur; cur != end(); ++cur)
- {
- if ((*cur == j) || (cmp(*cur, j) && cmp(i, *cur)))
- {
- if (cmp(*cur, *first))
- first = &(*cur);
- }
- }
- return *first;
-}
-
-template<typename I, class OrderCmp, class ConsistencyCmp>
-void
-Cycle<I,OrderCmp,ConsistencyCmp>::
-erase_between(const_reference i, const_reference j, const OrderComparison& cmp)
-{
- for (iterator cur = begin(); cur != end(); ++cur)
- while ((cur != end()) && ((*cur == j) || (cmp(*cur, j) && cmp(i, *cur))))
- {
- rLog(rlCycle, "Iteration of the erase while loop");
- cur = erase(cur);
- }
-}
-
-template<typename I, class OrderCmp, class ConsistencyCmp>
-std::ostream&
-Cycle<I,OrderCmp,ConsistencyCmp>::
-operator<<(std::ostream& out) const
-{
- for (const_iterator cur = begin(); cur != end(); ++cur)
- {
- if (cur != begin()) out << ", ";
- out << **cur;
- }
- // out << "(last: " << *last << ")"; // For debugging only
- return out;
-}
-
-template<typename I, class OrderCmp, class ConsistencyCmp>
-std::ostream&
-operator<<(std::ostream& out, const Cycle<I, OrderCmp, ConsistencyCmp>& c)
-{
- return c.operator<<(out);
-}
-
-template<typename I, class OrderCmp, class ConsistencyCmp>
-typename Cycle<I, OrderCmp, ConsistencyCmp>::Self&
-Cycle<I, OrderCmp, ConsistencyCmp>::
-add(const Self& c, const ConsistencyCmp& cmp)
-{
- rLog(rlCycle, "Adding cycles: %s + %s", tostring(*this).c_str(), tostring(c).c_str());
-
- iterator cur1 = begin();
- const_iterator cur2 = c.begin();
-
- while (cur2 != c.end())
- {
- if (cur1 == end())
- {
- while (cur2 != c.end())
- {
- push_back(*cur2++);
- Count(cCycleAddBasic);
- }
- rLog(rlCycle, "After addition: %s", tostring(*this).c_str());
- return *this;
- }
-
- // mod 2
- int res = cmp.compare(*cur1, *cur2);
- Count(cCycleAddComparison);
- rLog(rlCycle, "Comparison result: %i", res);
- if (res == 0) // *cur1 == *cur2
- {
- rLog(rlCycle, "Equality");
- cur1 = erase(cur1); // erase cur1 --- as a result cur1 will be pointing at old_cur1++
- ++cur2;
- } else if (res < 0) // *cur1 < *cur2
- {
- rLog(rlCycle, "Less than");
- cur1++;
- } else if (res > 0) // *cur1 > *cur2
- {
- rLog(rlCycle, "Greater than");
- insert(cur1, *cur2);
- ++cur2;
- }
- Count(cCycleAddBasic);
- }
-
- rLog(rlCycle, "After addition: %s", tostring(*this).c_str());
- return *this;
-}
-
-template<typename I, class OrderCmp, class ConsistencyCmp>
-typename Cycle<I,OrderCmp,ConsistencyCmp>::const_reference
-Cycle<I,OrderCmp,ConsistencyCmp>::
-add_and_first_between(const Self& c, const ConsistencyComparison& consistency_cmp,
- const_reference i, const_reference j, const OrderComparison& order_cmp)
-{
- add(c, consistency_cmp);
- return first_between(i,j, order_cmp);
-}
-
-template<typename I, class OrderCmp, class ConsistencyCmp>
-typename Cycle<I,OrderCmp,ConsistencyCmp>::const_reference
-Cycle<I,OrderCmp,ConsistencyCmp>::
-get_first(const OrderComparison& cmp) const
-{ return top(cmp); }
-
-
-template<typename I, class OrderCmp, class ConsistencyCmp>
-template<class Archive>
-void
-Cycle<I,OrderCmp,ConsistencyCmp>::
-serialize(Archive& ar, version_type )
-{
- ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP(Parent);
-}
-
-
-/*
-template<typename I, class Cmp>
-void
-Cycle<I, Cmp>::
-insertion_sort(const Comparison& cmp)
-{
- TemporaryCycleRepresenation tmp;
-
- // Insertion sort into the temporary list
- for (const_iterator cur = begin(); cur != end(); ++cur)
- {
- typename TemporaryCycleRepresenation::iterator tmp_cur = tmp.end();
- typename TemporaryCycleRepresenation::iterator tmp_next = tmp_cur--;
-
- while (tmp_next != tmp.begin())
- {
- if (cmp(*cur, *tmp_cur))
- tmp_next = tmp_cur--;
- else
- break;
- }
- tmp.insert(tmp_next, *cur);
- }
-
- // Copy temporary list back into ourselves
- iterator cur = begin();
- typename TemporaryCycleRepresenation::const_iterator tmp_cur = tmp.begin();
- while(tmp_cur != tmp.end())
- {
- *cur = *tmp_cur;
- ++cur; ++tmp_cur;
- }
-}
-*/
-
-
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/cycles.h Thu Dec 18 16:43:42 2008 -0800
@@ -0,0 +1,44 @@
+#ifndef __CYCLES_H__
+#define __CYCLES_H__
+
+#include "chain.h"
+#include <vector>
+#include "utilities/circular_list.h"
+
+template<class OrderIndex_ = int>
+struct VectorChains
+{
+ typedef OrderIndex_ OrderIndex;
+ typedef ChainWrapper<std::vector<OrderIndex> > Chain;
+ typedef Chain Cycle;
+
+ Cycle cycle;
+
+ VectorChains() {}
+ VectorChains(Cycle z): cycle(z) {}
+
+ bool sign() const { return cycle.empty(); }
+
+ template<class U> struct rebind
+ { typedef VectorChains<U> other; };
+};
+
+template<class OrderIndex_ = int>
+struct ListChains
+{
+ typedef OrderIndex_ OrderIndex;
+ typedef ChainWrapper<List<OrderIndex> > Chain;
+ typedef Chain Cycle;
+
+ Cycle cycle;
+
+ ListChains() {}
+ ListChains(Cycle z): cycle(z) {}
+
+ bool sign() const { return cycle.empty(); }
+
+ template<class U> struct rebind
+ { typedef ListChains<U> other; };
+};
+
+#endif // __CYCLES_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/dynamic-persistence.h Thu Dec 18 16:43:42 2008 -0800
@@ -0,0 +1,144 @@
+#ifndef __DYNAMIC_PERSISTENCE_H__
+#define __DYNAMIC_PERSISTENCE_H__
+
+#include "static-persistence.h"
+#include <utilities/types.h>
+
+#ifdef COUNTERS
+static Counter* cTrailLength = GetCounter("persistence/pair/traillength"); // the size of matrix U in RU decomposition
+#endif // COUNTERS
+
+template<class Data, class OrderDescriptor_, class ConsistencyIndex>
+struct TrailData_: public Data
+{
+ typedef TrailData_<Data, OrderDescriptor_, ConsistencyIndex> Self;
+
+ typedef typename OrderTraits<typename OrderDescriptor_::
+ template RebindData<Self>::
+ other>::Index OrderIndex;
+ typedef typename OrderDescriptor_::
+ Chains::template rebind<OrderIndex>::other::Chain Trail;
+
+ template<class Comparison>
+ struct ConsistencyComparison
+ {
+ typedef const OrderIndex& first_argument_type;
+ typedef const OrderIndex& second_argument_type;
+ typedef bool result_type;
+
+ ConsistencyComparison(const Comparison& cmp = Comparison()): cmp_(cmp) {}
+
+ bool operator()(const OrderIndex& a, const OrderIndex& b) const { return cmp_(a->consistency, b->consistency); }
+
+ Comparison cmp_;
+ };
+
+ Trail trail;
+ ConsistencyIndex consistency;
+};
+
+/**
+ * Class: DynamicPersistenceTrails
+ * Derives from StaticPersistence and allows one to update persistence
+ * after a transposition of two contiguous simplices in a filtration.
+ * In addition to reduced cycles, it stores each OrderElement's trails,
+ * i.e. in addition to matrix R, it stores matrix U in vineyard notation.
+ *
+ * Template parameters:
+ * Data_ - auxilliary contents to store with each OrderElement
+ * OrderDescriptor_ - class describing how the order is stored; it defaults to <VectorOrderDescriptor>
+ * which serves as a prototypical class
+ */
+// TODO: perhaps Consistency should be wrapped into a ConsistencyDescriptor that somehow knows how to initialize it.
+// That way one could provide a simple consistency descriptor that just stored some integers describing the original
+// position, or one could provide consistency that is references into the complex
+template<class Data_ = Empty,
+ class OrderDescriptor_ = VectorOrderDescriptor<>,
+ class ConsistencyIndex_ = size_t,
+ class ConsistencyComparison_ = std::less<ConsistencyIndex_> >
+class DynamicPersistenceTrails:
+ public StaticPersistence<TrailData_<Data_, OrderDescriptor_, ConsistencyIndex_>,
+ OrderDescriptor_>
+{
+ public:
+ typedef Data_ Data;
+ typedef TrailData_<Data_, OrderDescriptor_, ConsistencyIndex_> TrailData;
+ typedef StaticPersistence<TrailData, OrderDescriptor_> Parent;
+
+ typedef typename Parent::Traits Traits;
+ typedef typename Parent::OrderDescriptor OrderDescriptor;
+ typedef typename Parent::OrderComparison OrderComparison;
+ typedef typename Parent::OrderIndex OrderIndex;
+ typedef ConsistencyIndex_ ConsistencyIndex;
+ typedef ThreeOutcomeCompare<
+ typename TrailData::
+ template ConsistencyComparison<ConsistencyComparison_> > ConsistencyComparison;
+
+ /**
+ * Constructor: DynamicPersistenceTrails()
+ * TODO: write a description
+ *
+ * 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_());
+
+ 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 Visitor>
+ bool transpose(OrderIndex i, const Visitor& visitor = Visitor());
+
+ using Parent::begin;
+ using Parent::end;
+
+ // Struct: TranspositionVisitor
+ //
+ // For example, a VineardVisitor could implement this archetype.
+ struct TranspositionVisitor
+ {
+ // Function: transpose(i)
+ // 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 {}
+
+ // 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 {}
+ };
+
+ protected:
+ using Parent::order;
+
+ private:
+ void swap(OrderIndex i, OrderIndex j);
+ void pairing_switch(OrderIndex i, OrderIndex j);
+
+ struct PairingTrailsVisitor: public Parent::PairVisitor
+ {
+ // TODO: this is specialized for std::vector
+ PairingTrailsVisitor(OrderIndex bg, ConsistencyComparison ccmp):
+ bg_(bg), 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); }
+
+ OrderIndex bg_;
+ ConsistencyComparison ccmp_;
+ };
+
+ ConsistencyComparison ccmp_;
+};
+
+#include "dynamic-persistence.hpp"
+
+#endif // __DYNAMIC_PERSISTENCE_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/dynamic-persistence.hpp Thu Dec 18 16:43:42 2008 -0800
@@ -0,0 +1,258 @@
+#include <utilities/log.h>
+
+#ifdef LOGGING
+static rlog::RLogChannel* rlTranspositions = DEF_CHANNEL("topology/persistence/transpositions", rlog::Log_Debug);
+#endif // LOGGING
+
+#ifdef COUNTERS
+static Counter* cTransposition = GetCounter("persistence/transposition");
+static Counter* cTranspositionDiffDim = GetCounter("persistence/transposition/diffdim");
+static Counter* cTranspositionCase12 = GetCounter("persistence/transposition/case/1/2");
+static Counter* cTranspositionCase12s = GetCounter("persistence/transposition/case/1/2/special");
+static Counter* cTranspositionCase112 = GetCounter("persistence/transposition/case/1/1/2");
+static Counter* cTranspositionCase111 = GetCounter("persistence/transposition/case/1/1/1");
+static Counter* cTranspositionCase22 = GetCounter("persistence/transposition/case/2/2");
+static Counter* cTranspositionCase212 = GetCounter("persistence/transposition/case/2/1/2");
+static Counter* cTranspositionCase211 = GetCounter("persistence/transposition/case/2/1/1");
+static Counter* cTranspositionCase32 = GetCounter("persistence/transposition/case/3/2");
+static Counter* cTranspositionCase31 = GetCounter("persistence/transposition/case/3/1");
+static Counter* cTranspositionCase4 = GetCounter("persistence/transposition/case/4");
+#endif // COUNTERS
+
+
+template<class D, class OD, class CI, class CC>
+template<class Filtration>
+DynamicPersistenceTrails<D,OD,CI,CC>::
+DynamicPersistenceTrails(const Filtration& f, const OrderComparison& ocmp, const CC& ccmp):
+ Parent(f, ocmp), ccmp_(ccmp)
+{}
+
+template<class D, class OD, class CI, class CC>
+void
+DynamicPersistenceTrails<D,OD,CI,CC>::
+pair_simplices()
+{
+ Parent::pair_simplices(begin(), end(), PairingTrailsVisitor(begin(), ccmp_));
+}
+
+template<class D, class OD, class CI, class CC>
+template<class Visitor>
+bool
+DynamicPersistenceTrails<D,OD,CI,CC>::
+transpose(OrderIndex i, const Visitor& visitor)
+{
+#if LOGGING
+ typename Traits::OutputMap outmap(order());
+#endif
+
+ Count(cTransposition);
+ typedef OrderIndex Index;
+ typedef typename TrailData::Trail::iterator TrailIterator;
+
+ visitor.transpose(i);
+
+ Index i_prev = i++;
+
+#if 0 // Persistence no longer has the notion of dimension
+ if (i_prev->dimension() != i->dimension())
+ {
+ swap(i_prev, i);
+ rLog(rlTranspositions, "Different dimension");
+ Count(cTranspositionDiffDim);
+ return false;
+ }
+#endif
+
+ bool si = i_prev->sign(), sii = i->sign();
+ if (si && sii)
+ {
+ 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)
+ {
+ rLog(rlTranspositions, "Case 1, U[i,i+1] = 1");
+ i_prev->trail.erase(*i_in_i_prev);
+ }
+
+ Index k = i_prev->pair;
+ Index l = i->pair;
+
+ // Explicit treatment of unpaired simplex
+ if (l == i)
+ {
+ swap(i_prev, i);
+ rLog(rlTranspositions, "Case 1.2 --- unpaired");
+ rLog(rlTranspositions, "%s", outmap(i_prev).c_str());
+ Count(cTranspositionCase12);
+ return false;
+ } else if (k == i_prev)
+ {
+ if (!(l->cycle.contains(i_prev)))
+ {
+ // Case 1.2
+ swap(i_prev, i);
+ rLog(rlTranspositions, "Case 1.2 --- unpaired");
+ rLog(rlTranspositions, outmap(i_prev).c_str());
+ Count(cTranspositionCase12);
+ return false;
+ } else
+ {
+ // Case 1.2 --- special version (plain swap, but pairing switches)
+ swap(i_prev, i);
+ pairing_switch(i_prev, i);
+ visitor.switched(i_prev, Case12);
+ rLog(rlTranspositions, "Case 1.2 --- unpaired (pairing switch)");
+ rLog(rlTranspositions, outmap(i_prev).c_str());
+ Count(cTranspositionCase12s);
+ return true;
+ }
+ }
+
+ rLog(rlTranspositions, "l cycle: %s", l->cycle.tostring(outmap).c_str());
+ if (!(l->cycle.contains(i_prev)))
+ {
+ // Case 1.2
+ swap(i_prev, i);
+ rLog(rlTranspositions, "Case 1.2");
+ Count(cTranspositionCase12);
+ return false;
+ } else
+ {
+ // Case 1.1
+ if (not2(ccmp_)(k,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
+ rLog(rlTranspositions, "Case 1.1.1");
+ Count(cTranspositionCase111);
+ return false;
+ } else
+ {
+ // 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
+ pairing_switch(i_prev, i);
+ visitor.switched(i_prev, Case112);
+ rLog(rlTranspositions, "Case 1.1.2");
+ Count(cTranspositionCase112);
+ return true;
+ }
+ }
+ } else if (!si && !sii)
+ {
+ // Case 2
+ if (!(i_prev->trail.contains(i)))
+ {
+ // Case 2.2
+ swap(i_prev, i);
+ rLog(rlTranspositions, "Case 2.2");
+ Count(cTranspositionCase22);
+ return false;
+ } 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
+ swap(i_prev, i);
+ if (not2(ccmp_)(low_ii, 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
+ pairing_switch(i_prev, i);
+ visitor.switched(i_prev, Case212);
+ rLog(rlTranspositions, "Case 2.1.2");
+ Count(cTranspositionCase212);
+ return true;
+ }
+
+ // Case 2.1.1
+ rLog(rlTranspositions, "Case 2.1.1");
+ Count(cTranspositionCase211);
+ return false;
+ }
+ } else if (!si && sii)
+ {
+ // Case 3
+ if (!(i_prev->trail.contains(i)))
+ {
+ // Case 3.2
+ swap(i_prev, i);
+ rLog(rlTranspositions, "Case 3.2");
+ Count(cTranspositionCase32);
+ return false;
+ } 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
+ 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
+ pairing_switch(i_prev, i);
+ visitor.switched(i_prev, Case31);
+ rLog(rlTranspositions, "Case 3.1");
+ Count(cTranspositionCase31);
+ return true;
+ }
+ } else if (si && !sii)
+ {
+ // Case 4
+ boost::optional<TrailIterator> i_in_i_prev = i_prev->trail.contains(i);
+ if (i_in_i_prev)
+ {
+ rLog(rlTranspositions, "Case 4, U[i,i+1] = 1");
+ i_prev->trail.erase(*i_in_i_prev);
+ }
+ swap(i_prev, i);
+ rLog(rlTranspositions, "Case 4");
+ Count(cTranspositionCase4);
+ return false;
+ }
+
+ return false; // to avoid compiler complaints; we should never reach this point
+}
+
+template<class D, class OD, class CI, class CC>
+void
+DynamicPersistenceTrails<D,OD,CI,CC>::
+swap(OrderIndex i, OrderIndex 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);
+}
+
+template<class D, class OD, class CI, class CC>
+void
+DynamicPersistenceTrails<D,OD,CI,CC>::
+pairing_switch(OrderIndex i, OrderIndex j)
+{
+ OrderIndex i_pair = i->pair;
+ OrderIndex j_pair = j->pair;
+
+ if (i_pair == i)
+ j->pair = j;
+ else
+ {
+ j->pair = i_pair;
+ i_pair->pair = j;
+ }
+
+ if (j_pair == j)
+ i->pair = i;
+ else
+ {
+ i->pair = j_pair;
+ j_pair->pair = i;
+ }
+}
--- a/include/topology/filtration.h Mon Sep 22 21:53:25 2008 -0700
+++ b/include/topology/filtration.h Thu Dec 18 16:43:42 2008 -0800
@@ -1,137 +1,99 @@
-/*
- * Author: Dmitriy Morozov
- * Department of Computer Science, Duke University, 2005 -- 2006
- */
-
#ifndef __FILTRATION_H__
#define __FILTRATION_H__
-#include "utilities/log.h"
+#include <vector>
+#include <iostream>
+#include "complex-traits.h"
+#include "utilities/indirect.h"
+#include "utilities/property-maps.h"
-#include "filtrationcontainer.h"
-#include "filtrationsimplex.h"
-#include "vineyard.h"
-
-#include <map>
-#include <vector>
-
-#include <boost/serialization/access.hpp>
-/**
- * Filtration class. Serves as an (ordered) container for the simplices,
- * and provides pair_simplices() method that computes the RU-decomposition
- * for the simplex order stored in the filtration. Iterators remain valid
- * through all the operations.
- *
- * \ingroup topology
- */
-template<class Smplx, class FltrSmplx = FiltrationSimplex<Smplx>, class Vnrd = Vineyard<FltrSmplx> >
-class Filtration: public FltrSmplx::Container
+// Class: Filtration
+//
+// 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_> >
+class Filtration
{
- public:
- typedef Smplx Simplex;
- typedef FltrSmplx FiltrationSimplex;
- typedef Vnrd Vineyard;
-
- /// \name Container Types
- /// @{
- /** The actual container type (which is the parent of the Filtration) */
- typedef typename FiltrationSimplex::Container FiltrationContainer;
- typedef typename FiltrationContainer::Index Index;
- typedef typename FiltrationContainer::const_Index const_Index;
- typedef Index iterator;
- typedef const_Index const_iterator;
- /// @}
-
- /// \name Cycles and Trails
- /// @{
- typedef typename FiltrationContainer::GreaterThanComparison CyclesComparator;
- typedef typename FiltrationContainer::LessThanComparison TrailsComparator;
- typedef typename FiltrationContainer::ConsistencyComparison ConsistencyComparator;
- typedef typename FiltrationContainer::Cycle Cycle;
- typedef typename FiltrationContainer::Trail Trail;
- typedef typename Cycle::iterator CycleIterator;
- typedef typename Trail::iterator TrailIterator;
- /// @}
-
- typedef Filtration<Simplex, FiltrationSimplex, Vineyard> Self;
- typedef FiltrationContainer Parent;
+ public:
+ // Typedefs: Template parameters
+ typedef Index_ IntermediateIndex;
+ typedef Complex_ Complex;
+ typedef ComplexTraits_ ComplexTraits;
+
+ // Typedefs: Complex
+ typedef typename ComplexTraits::Index ComplexIndex;
+ typedef typename ComplexTraits::Simplex Simplex;
+ typedef typename ComplexTraits::SimplexIndexMap SimplexIndexMap;
+ typedef typename Simplex::Boundary SimplexBoundary;
+ typedef std::vector<IntermediateIndex> IndexBoundary;
+
+ // Typedefs: Order
+ typedef std::vector<ComplexIndex> Order;
+ typedef typename Order::const_iterator Index;
+ typedef std::vector<IntermediateIndex> ReverseOrder;
+ typedef typename ReverseOrder::const_iterator ReverseOrderIndex;
+
+ // 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; }
+
+ // 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;
- public:
- Filtration(Vineyard* vineyard);
-
- /// \name Core Functionality
- /// @{
- /// Computes RU decomposition of the simplices in [bg, end) range, assuming that everything before bg has been paired
- void pair_simplices(Index bg, Index end, bool store_trails = true);
- void pair_simplices(bool store_trails = true) { pair_simplices(begin(), end(), store_trails); }
- bool transpose(Index i, bool maintain_lazy = false);
- bool is_paired() const;
- Index append(const Simplex& s); ///< Appends s to the filtration
- Index insert(Index prior, const Simplex& s); ///< Inserts s after prior
- const_Index get_index(const Simplex& s) const; /**< Returns the iterator pointing to s
- (end() if s not in filtration) */
- Index get_index(const Simplex& s); ///< \overload
- void fill_simplex_index_map(); ///< Fills the mapping for get_index()
- /// @}
-
- /// \name Accessors
- /// @{
- Vineyard* vineyard() { return vineyard_; }
- const Vineyard* vineyard() const { return vineyard_; }
- /// @}
-
- protected:
- using Parent::swap;
- bool transpose_simplices(Index i, bool maintain_lazy);
+ Index begin() const { return order_.begin(); }
+ Index end() const { return order_.end(); }
+ size_t size() const { return order_.size(); }
+
+ std::ostream& operator<<(std::ostream& out) const;
+
+ private:
+ Order order_;
+ ReverseOrder reverse_order_;
+ OffsetMap<ComplexIndex,
+ ReverseOrderIndex>
+ complex_order_map_;
+ SimplexIndexMap simplex_index_map_;
+};
+
+template<class C, class I, class CT>
+std::ostream&
+operator<<(std::ostream& out, const Filtration<C,I,CT>& f) { return f.operator<<(out); }
+
- public:
- /// \name Container Operations
- /// @{
- using Parent::size;
- using Parent::begin;
- using Parent::end;
- /// @}
-
- std::ostream& operator<<(std::ostream& out) const;
+template<class Functor_, class Filtration_>
+class ThroughFiltration
+{
+ public:
+ typedef Filtration_ Filtration;
+ typedef Functor_ Functor;
- public: // doesn't really need to be public, except for assertions
- /// \name Comparator accessors
- /// @{
- const ConsistencyComparator& get_consistency_cmp() const { return consistency_cmp; }
- const CyclesComparator& get_cycles_cmp() const { return cycles_cmp; }
- const TrailsComparator& get_trails_cmp() const { return trails_cmp; }
- /// @}
-
- private:
- typedef std::map<Simplex, Index> SimplexMap;
+ typedef typename Functor::result_type result_type;
+ typedef typename Filtration::Index first_argument_type;
- /// Initializes the cycle with the indices of the simplices in the boundary, and the trail with the index of this simplex
- void init_cycle_trail(Index j);
- void pairing_switch(Index i, Index j);
-
- bool paired;
- bool trails_stored;
- SimplexMap inverse_simplices;
+ ThroughFiltration(const Filtration& filtration,
+ const Functor& functor):
+ filtration_(filtration),
+ functor_(functor) {}
- Vineyard* vineyard_;
+ result_type operator()(first_argument_type a) const { return functor_(filtration_.simplex(a)); }
- CyclesComparator cycles_cmp;
- TrailsComparator trails_cmp;
- ConsistencyComparator consistency_cmp;
+ private:
+ const Filtration& filtration_;
+ const Functor& functor_;
+};
- private:
- /* Serialization */
- friend class boost::serialization::access;
-
- typedef std::map<const_Index, SizeType, ConsistencyComparator> IndexIntMap;
- typedef std::vector<Index> IndexVector;
-
- template<class Archive> void save(Archive& ar, version_type ) const;
- template<class Archive> void load(Archive& ar, version_type );
- BOOST_SERIALIZATION_SPLIT_MEMBER()
-};
+template<class Filtration, class Functor>
+ThroughFiltration<Functor, Filtration>
+evaluate_through_filtration(const Filtration& filtration, const Functor& functor)
+{ return ThroughFiltration<Functor, Filtration>(filtration, functor); }
#include "filtration.hpp"
-#endif // __FILTRATION_H__
+#endif // __FILTRATION_H__
--- a/include/topology/filtration.hpp Mon Sep 22 21:53:25 2008 -0700
+++ b/include/topology/filtration.hpp Thu Dec 18 16:43:42 2008 -0800
@@ -1,508 +1,46 @@
-#include "utilities/types.h"
-#include "utilities/counter.h"
-#include <algorithm>
-
-#include <boost/utility.hpp>
-#include <boost/serialization/vector.hpp>
-#include <boost/serialization/nvp.hpp>
-#include <boost/serialization/list.hpp>
-
-using boost::serialization::make_nvp;
-
-/* Filtration Public */
-
-#ifdef LOGGING
-static rlog::RLogChannel* rlFiltration = DEF_CHANNEL("topology/filtration", rlog::Log_Debug);
-static rlog::RLogChannel* rlFiltrationTranspositions = DEF_CHANNEL("topology/filtration/transpositions", rlog::Log_Debug);
-#endif // LOGGING
-
-#ifdef COUNTERS
-static Counter* cFiltrationPair = GetCounter("filtration/pair");
-static Counter* cFiltrationPairBoundaries = GetCounter("filtration/pair/boundaries");
-static Counter* cFiltrationPairCycleLength = GetCounter("filtration/pair/cyclelength");
-static Counter* cFiltrationPairTrailLength = GetCounter("filtration/pair/traillength");
-static Counter* cFiltrationTransposition = GetCounter("filtration/transposition");
-static Counter* cFiltrationTranspositionDiffDim = GetCounter("filtration/transposition/diffdim");
-static Counter* cFiltrationTranspositionCase12 = GetCounter("filtration/transposition/case/1/2");
-static Counter* cFiltrationTranspositionCase12s = GetCounter("filtration/transposition/case/1/2/special");
-static Counter* cFiltrationTranspositionCase112 = GetCounter("filtration/transposition/case/1/1/2");
-static Counter* cFiltrationTranspositionCase111 = GetCounter("filtration/transposition/case/1/1/1");
-static Counter* cFiltrationTranspositionCase22 = GetCounter("filtration/transposition/case/2/2");
-static Counter* cFiltrationTranspositionCase212 = GetCounter("filtration/transposition/case/2/1/2");
-static Counter* cFiltrationTranspositionCase211 = GetCounter("filtration/transposition/case/2/1/1");
-static Counter* cFiltrationTranspositionCase32 = GetCounter("filtration/transposition/case/3/2");
-static Counter* cFiltrationTranspositionCase31 = GetCounter("filtration/transposition/case/3/1");
-static Counter* cFiltrationTranspositionCase4 = GetCounter("filtration/transposition/case/4");
-#endif // COUNTERS
-
-
-template<class S, class FS, class V>
-Filtration<S, FS, V>::
-Filtration(Vineyard* vnrd = 0): paired(false), vineyard_(vnrd)
-{}
-
-template<class S, class FS, class V>
-void
-Filtration<S, FS, V>::
-pair_simplices(Index bg, Index end, bool store_trails)
-{
- if (!is_paired())
- trails_stored = store_trails;
- else
- trails_stored &= store_trails;
-
- rLog(rlFiltration, "Entered: compute_pairing");
- for (Index j = bg; j != end; ++j)
- {
- rLog(rlFiltration, "Simplex %s", tostring(*j).c_str());
- init_cycle_trail(j);
- Cycle& bdry = j->cycle();
- rLog(rlFiltration, " has boundary: %s", tostring(bdry).c_str());
-
- CountNum(cFiltrationPairBoundaries, j->dimension());
- Count(cFiltrationPair);
-
- while(!bdry.empty())
- {
- Index i = bdry.top(cycles_cmp);
- rLog(rlFiltration, "%s: %s", tostring(*i).c_str(), tostring(*(i->pair())).c_str());
- AssertMsg(!cycles_cmp(i, j),
- "Simplices in the cycle must precede current simplex: (%s in cycle of %s)",
- tostring(*i).c_str(), tostring(*j).c_str());
-
- // i is not paired, so we pair j with i
- if (i->pair() == i)
- {
- rLog(rlFiltration, "Pairing %s and %s with cycle %s",
- tostring(*i).c_str(), tostring(*j).c_str(),
- tostring(j->cycle()).c_str());
- i->set_pair(j);
- j->set_pair(i);
- CountNum(cFiltrationPairCycleLength, j->cycle().size());
- CountBy(cFiltrationPairCycleLength, j->cycle().size());
- break;
- }
-
- rLog(rlFiltration, " Adding: [%s] + [%s]",
- tostring(bdry).c_str(), tostring(i->pair()->cycle()).c_str());
- bdry.add(i->pair()->cycle(), get_consistency_cmp());
- if (store_trails) i->pair()->trail().append(j, get_consistency_cmp());
- Count(cFiltrationPairTrailLength);
- rLog(rlFiltration, "After addition: %s", tostring(bdry).c_str());
- }
- rLog(rlFiltration, "Finished with %s: %s",
- tostring(*j).c_str(), tostring(*(j->pair())).c_str());
- }
- paired = true;
-}
+#include "utilities/containers.h"
-template<class S, class FS, class V>
-bool
-Filtration<S, FS, V>::
-is_paired() const
-{ return paired; }
-
-/**
- * Transposes simplices at i and i+1, and records the knee in the vineyard if there is a change in pairing.
- * Returns true if the pairing changed.
- */
-template<class S, class FS, class V>
-bool
-Filtration<S,FS,V>::
-transpose(Index i, bool maintain_lazy)
+template<class C, class I, class CT>
+template<class Comparison>
+Filtration<C,I,CT>::
+Filtration(ComplexIndex bg, ComplexIndex end, const Comparison& cmp):
+ order_(RecursiveIterator<ComplexIndex>(bg),
+ RecursiveIterator<ComplexIndex>(end)),
+ reverse_order_(order_.size()),
+ complex_order_map_(bg, reverse_order_.begin()),
+ simplex_index_map_(bg, end)
{
- AssertMsg(vineyard() != 0, "We must have a vineyard for transpositions");
- AssertMsg(trails_stored, "We must have trails (matrix U) to perform transpositions");
-
- Index i_orig = i++;
-
- rLog(rlFiltrationTranspositions, "Transposing:");
- rLog(rlFiltrationTranspositions, " %s: (%s | %s)",
- tostring(*i_orig).c_str(), tostring(i_orig->cycle()).c_str(),
- tostring(i_orig->trail()).c_str());
- rLog(rlFiltrationTranspositions, " and");
- rLog(rlFiltrationTranspositions, " %s: (%s | %s)",
- tostring(*i).c_str(), tostring(i->cycle()).c_str(),
- tostring(i->trail()).c_str());
-
- AssertMsg(i_orig->pair() != i, "Transposing simplices must not be paired");
- bool result = transpose_simplices(i_orig, maintain_lazy);
- AssertMsg(i_orig == boost::next(i), "Wrong indices after transposition");
-
- if (result) vineyard()->switched(i, i_orig);
- return result;
-}
-
-template<class S, class FS, class V>
-typename Filtration<S, FS, V>::Index
-Filtration<S, FS, V>::
-append(const Simplex& s)
-{
- Index i = push_back(FiltrationSimplex(s));
- return i;
-}
-
-template<class S, class FS, class V>
-typename Filtration<S, FS, V>::Index
-Filtration<S, FS, V>::
-insert(Index prior, const Simplex& s)
-{
- Index i = Parent::insert(prior, FiltrationSimplex(s));
- paired = false;
-
- return i;
-}
-
-template<class S, class FS, class V>
-typename Filtration<S, FS, V>::const_Index
-Filtration<S, FS, V>::
-get_index(const Simplex& s) const
-{
- typename SimplexMap::const_iterator i = inverse_simplices.find(s);
- if (i == inverse_simplices.end())
- return end();
- else
- return i->second;
-}
-
-template<class S, class FS, class V>
-typename Filtration<S, FS, V>::Index
-Filtration<S, FS, V>::
-get_index(const Simplex& s)
-{
- typename SimplexMap::const_iterator i = inverse_simplices.find(s);
- if (i == inverse_simplices.end())
- return end();
- else
- return i->second;
-}
-
-template<class S, class FS, class V>
-void
-Filtration<S, FS, V>::
-fill_simplex_index_map()
-{
- for (Index i = begin(); i != end(); ++i)
- inverse_simplices[*i] = i;
-}
-
-template<class S, class FS, class V>
-std::ostream&
-Filtration<S, FS, V>::
-operator<<(std::ostream& out) const
-{
- out << "Pairing: " << std::endl;
- for (const_Index i = begin(); i != end(); ++i)
- {
- out << "(" << *i << ", " << *(i->pair()) << "): ";
- out << i->cycle() << std::endl;
- }
- out << std::endl << std::endl;
-
- return out;
+ std::sort(order_.begin(), order_.end(), IndirectComparison<ComplexIndex, Comparison>(cmp));
+ for (Index obg = order_.begin(), cur = obg; cur != order_.end(); ++cur)
+ reverse_order_[*cur - bg] = cur - obg;
}
-
-/* Filtration Protected */
-/// Transposes simplices at i and i+1. Returns true if the pairing switched.
-template<class S, class FS, class V>
-bool
-Filtration<S,FS,V>::
-transpose_simplices(Index i, bool maintain_lazy)
+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(is_paired(), "Pairing must be computed before transpositions");
- Count(cFiltrationTransposition);
-
- Index i_prev = i++;
-
- if (i_prev->dimension() != i->dimension())
- {
- swap(i_prev, i);
- rLog(rlFiltrationTranspositions, "Different dimension");
- Count(cFiltrationTranspositionDiffDim);
- return false;
- }
-
- bool si = i_prev->sign(), sii = i->sign();
- if (si && sii)
- {
- rLog(rlFiltrationTranspositions, "Trail prev: %s", tostring(i_prev->trail()).c_str());
-
- // Case 1
- TrailIterator i_in_i_prev = std::find(i_prev->trail().begin(), i_prev->trail().end(), i);
- if (i_in_i_prev != i_prev->trail().end())
- {
- rLog(rlFiltrationTranspositions, "Case 1, U[i,i+1] = 1");
- i_prev->trail().erase(i_in_i_prev);
- }
-
- Index k = i_prev->pair();
- Index l = i->pair();
+ AssertMsg(bdry.empty(), "We are initializing the boundary from scratch");
+ SimplexBoundary simplex_bdry = (*i)->boundary();
+ ContainerTraits<Cycle>::reserve(bdry, simplex_bdry.size());
+ typename Map::template rebind_from<IntermediateIndex>::other order_bdry_map(0, map.to());
- // Explicit treatment of unpaired simplex
- if (l == i)
- {
- swap(i_prev, i);
- rLog(rlFiltrationTranspositions, "Case 1.2 --- unpaired");
- rLog(rlFiltrationTranspositions, "%s", tostring(*i_prev).c_str());
- Count(cFiltrationTranspositionCase12);
- return false;
- } else if (k == i_prev)
- {
- if (std::find(l->cycle().begin(), l->cycle().end(), i_prev) == l->cycle().end())
- {
- // Case 1.2
- swap(i_prev, i);
- rLog(rlFiltrationTranspositions, "Case 1.2 --- unpaired");
- rLog(rlFiltrationTranspositions, tostring(*i_prev).c_str());
- Count(cFiltrationTranspositionCase12);
- return false;
- } else
- {
- // Case 1.2 --- special version (plain swap, but pairing switches)
- swap(i_prev, i);
- pairing_switch(i_prev, i);
- rLog(rlFiltrationTranspositions, "Case 1.2 --- unpaired (pairing switch)");
- rLog(rlFiltrationTranspositions, tostring(*i_prev).c_str());
- Count(cFiltrationTranspositionCase12s);
- return true;
- }
- }
-
- rLog(rlFiltrationTranspositions, "l cycle: %s", tostring(l->cycle()).c_str());
- if (std::find(l->cycle().begin(), l->cycle().end(), i_prev) == l->cycle().end())
- {
- // Case 1.2
- if (maintain_lazy)
- {
- TrailIterator k_in_l = std::find(l->trail().begin(), l->trail().end(), k);
- if (k_in_l != l->trail().end())
- {
- l->trail().add(k->trail(), Filtration::get_consistency_cmp()); // Add row k to l
- k->cycle().add(l->cycle(), Filtration::get_consistency_cmp()); // Add column l to k
- }
- }
- swap(i_prev, i);
- rLog(rlFiltrationTranspositions, "Case 1.2");
- Count(cFiltrationTranspositionCase12);
- return false;
- } else
- {
- // Case 1.1
- if (trails_cmp(k,l))
- {
- // Case 1.1.1
- swap(i_prev, i);
- l->cycle().add(k->cycle(), Filtration::get_consistency_cmp()); // Add column k to l
- k->trail().add(l->trail(), Filtration::get_consistency_cmp()); // Add row l to k
- rLog(rlFiltrationTranspositions, "Case 1.1.1");
- Count(cFiltrationTranspositionCase111);
- return false;
- } else
- {
- // Case 1.1.2
- swap(i_prev, i);
- k->cycle().add(l->cycle(), Filtration::get_consistency_cmp()); // Add column l to k
- l->trail().add(k->trail(), Filtration::get_consistency_cmp()); // Add row k to l
- pairing_switch(i_prev, i);
- rLog(rlFiltrationTranspositions, "Case 1.1.2");
- Count(cFiltrationTranspositionCase112);
- return true;
- }
- }
- } else if (!si && !sii)
- {
- // Case 2
- if (std::find(i_prev->trail().begin(), i_prev->trail().end(), i) == i_prev->trail().end())
- {
- // Case 2.2
- swap(i_prev, i);
- rLog(rlFiltrationTranspositions, "Case 2.2");
- Count(cFiltrationTranspositionCase22);
- return false;
- } else
- {
- // Case 2.1
- Index low_i = i_prev->pair();
- Index low_ii = i->pair();
- i_prev->trail().add(i->trail(), Filtration::get_consistency_cmp()); // Add row i to i_prev
- i->cycle().add(i_prev->cycle(), Filtration::get_consistency_cmp()); // Add column i_prev to i
- swap(i_prev, i);
- if (Filtration::get_trails_cmp()(low_ii, low_i))
- {
- // Case 2.1.2
- i_prev->cycle().add(i->cycle(), Filtration::get_consistency_cmp()); // Add column i to i_prev (after transposition)
- i->trail().add(i_prev->trail(), Filtration::get_consistency_cmp()); // Add row i to i_prev
- pairing_switch(i_prev, i);
- rLog(rlFiltrationTranspositions, "Case 2.1.2");
- Count(cFiltrationTranspositionCase212);
- return true;
- }
-
- // Case 2.1.1
- rLog(rlFiltrationTranspositions, "Case 2.1.1");
- Count(cFiltrationTranspositionCase211);
- return false;
- }
- } else if (!si && sii)
- {
- // Case 3
- if (std::find(i_prev->trail().begin(), i_prev->trail().end(), i) == i_prev->trail().end())
- {
- // Case 3.2
- swap(i_prev, i);
- rLog(rlFiltrationTranspositions, "Case 3.2");
- Count(cFiltrationTranspositionCase32);
- return false;
- } else
- {
- // Case 3.1
- i_prev->trail().add(i->trail(), Filtration::get_consistency_cmp()); // Add row i to i_prev
- i->cycle().add(i_prev->cycle(), Filtration::get_consistency_cmp()); // Add column i_prev to i
- swap(i_prev, i);
- i_prev->cycle().add(i->cycle(), Filtration::get_consistency_cmp()); // Add column i_prev to i (after transposition)
- i->trail().add(i_prev->trail(), Filtration::get_consistency_cmp()); // Add row i to i_prev
- pairing_switch(i_prev, i);
- rLog(rlFiltrationTranspositions, "Case 3.1");
- Count(cFiltrationTranspositionCase31);
- return true;
- }
- } else if (si && !sii)
- {
- // Case 4
- TrailIterator i_in_i_prev = std::find(i_prev->trail().begin(), i_prev->trail().end(), i);
- if (i_in_i_prev != i_prev->trail().end())
- {
- rLog(rlFiltrationTranspositions, "Case 4, U[i,i+1] = 1");
- i_prev->trail().erase(i_in_i_prev);
- }
- swap(i_prev, i);
- rLog(rlFiltrationTranspositions, "Case 4");
- Count(cFiltrationTranspositionCase4);
- return false;
- }
-
- return false; // to avoid compiler complaints, should never reach this point
+ for (typename SimplexBoundary::const_iterator cur = simplex_bdry.begin(); cur != simplex_bdry.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;
+ }
}
-
-/* Filtration Private */
-template<class S, class FS, class V>
-void
-Filtration<S, FS, V>::
-init_cycle_trail(Index j)
-{
- typename Simplex::Cycle bdry = j->boundary();
-
- for (typename Simplex::Cycle::const_iterator cur = bdry.begin(); cur != bdry.end(); ++cur)
- {
- rLog(rlFiltration, "Appending in init_cycle_trail(): %s",
- tostring(*cur).c_str());
- AssertMsg(get_index(*cur) != end(), "Non-existent simplex in the cycle");
- j->cycle().append(get_index(*cur), get_consistency_cmp());
- }
- j->trail().append(j, get_consistency_cmp());
- j->set_pair(j);
-}
-
-/// Update the pairing, so that whoever was paired with i is now paired with j and vice versa.
-template<class S, class FS, class V>
-void
-Filtration<S,FS,V>::
-pairing_switch(Index i, Index j)
-{
- Index i_pair = i->pair();
- Index j_pair = j->pair();
-
- if (i_pair == i)
- j->set_pair(j);
- else
- {
- j->set_pair(i_pair);
- i_pair->set_pair(j);
- }
-
- if (j_pair == j)
- i->set_pair(i);
- else
- {
- i->set_pair(j_pair);
- j_pair->set_pair(i);
- }
-}
-
-/* Serialization */
-template<class S, class FS, class V>
-template<class Archive>
-void
-Filtration<S, FS, V>::
-save(Archive& ar, version_type ) const
+template<class C, class I, class CT>
+std::ostream&
+Filtration<C,I,CT>::
+operator<<(std::ostream& out) const
{
- ar << BOOST_SERIALIZATION_NVP(paired);
- ar << BOOST_SERIALIZATION_NVP(cycles_cmp);
- ar << BOOST_SERIALIZATION_NVP(trails_cmp);
- ar << BOOST_SERIALIZATION_NVP(consistency_cmp);
-
- SizeType sz = size();
- ar << make_nvp("size", sz);
- rLog(rlFiltration, "Size: %i", sz);
-
- /* Record integer indices */
- IndexIntMap index_map; SizeType i = 0;
- for (const_Index cur = begin(); cur != end(); ++cur)
- { index_map[cur] = i++; }
-
- /* Save the simplices */
- int count = 0;
- for (const_Index cur = begin(); cur != end(); ++cur)
- {
- count++;
- // FIXME
- //FiltrationSimplexSerialization simplex = FiltrationSimplexSerialization(*cur, index_map);
- //ar << make_nvp("FiltrationSimplex", simplex);
- }
- rLog(rlFiltration, "%i simplices serialized", count);
+ for (Index i = begin(); i != end(); ++i)
+ out << **i << std::endl;
+ return out;
}
-
-template<class S, class FS, class V>
-template<class Archive>
-void
-Filtration<S, FS, V>::
-load(Archive& ar, version_type )
-{
- rLog(rlFiltration, "Starting to read filtration");
- ar >> BOOST_SERIALIZATION_NVP(paired);
- ar >> BOOST_SERIALIZATION_NVP(cycles_cmp);
- ar >> BOOST_SERIALIZATION_NVP(trails_cmp);
- ar >> BOOST_SERIALIZATION_NVP(consistency_cmp);
- rLog(rlFiltration, "In Filtration: first block read");
-
- SizeType sz;
- ar >> make_nvp("size", sz);
- rLog(rlFiltration, "In Filtration: size read %i", sz);
-
- IndexVector index_vector(sz);
- for (SizeType i = 0; i < sz; ++i)
- {
- index_vector[i] = append(Simplex());
- }
-
- int count = 0;
- for (SizeType i = 0; i < sz; ++i)
- {
- // FIXME
- //FiltrationSimplexSerialization simplex;
- //ar >> make_nvp("FiltrationSimplex", simplex);
- count++;
- rLog(rlFiltration, "In Filtration: simplex read (%i)", count);
- //simplex.set_filtration_simplex(*index_vector[i], index_vector);
- }
- rLog(rlFiltration, "In Filtration: simplices read");
-}
-
-template<class S, class FS, class V>
-std::ostream&
-operator<<(std::ostream& out, const Filtration<S, FS, V>& f)
-{ return f.operator<<(out); }
-
-
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/order.h Thu Dec 18 16:43:42 2008 -0800
@@ -0,0 +1,156 @@
+#ifndef __ORDER_H__
+#define __ORDER_H__
+
+#include "cycles.h"
+#include "utilities/types.h"
+
+#include <vector>
+#include <list>
+
+//#include <iostream>
+#include <sstream>
+#include <string>
+
+template<class OrderDescriptor>
+struct OrderTraits {};
+
+/**
+ * Class: VectorOrderDescriptor
+ * Class that stores a dense order of simplices in a vector.
+ * That prevents it from performing efficient insertions.
+ */
+template<class Chains_ = VectorChains<>,
+ class Data_ = Empty>
+struct VectorOrderDescriptor:
+ public Chains_::template rebind<typename OrderTraits<VectorOrderDescriptor<Chains_, Data_> >::Index>::other,
+ public Data_
+{
+ typedef VectorOrderDescriptor<Chains_, Data_> Self;
+
+ typedef typename OrderTraits<Self>::Index OrderIndex;
+ typedef typename Chains_::template rebind<OrderIndex>::other Chains;
+ typedef Data_ Data;
+
+ template<class OtherData_> struct RebindData
+ { typedef VectorOrderDescriptor<Chains_, OtherData_> other; };
+
+
+ // TODO: Do we need any of these constructors?
+ VectorOrderDescriptor(const Data& d = Data()):
+ Data(d) {}
+ VectorOrderDescriptor(typename Chains::Cycle z, const Data& d = Data()):
+ Chains(z), Data(d) {}
+ VectorOrderDescriptor(OrderIndex i, const Data& d = Data()):
+ Data(d), pair(i) {}
+ VectorOrderDescriptor(OrderIndex i, typename Chains::Cycle z, const Data& d = Data()):
+ Chains(z), Data(d), pair(i) {}
+
+ OrderIndex pair;
+};
+
+// Specialization for VectorOrderDescriptor
+template<class Chains, class Data>
+struct OrderTraits<VectorOrderDescriptor<Chains, Data> >
+{
+ typedef VectorOrderDescriptor<Chains, Data> Descriptor;
+
+ typedef std::vector<Descriptor> Order;
+ typedef Descriptor Element;
+ typedef typename Order::iterator Index;
+ typedef ThreeOutcomeCompare<std::greater<Index> > Comparison;
+
+ class OutputMap
+ {
+ public:
+ OutputMap(const Order& order):
+ bg_(order.begin()) {}
+
+ // returns a string with (i - bg_)
+ std::string operator()(Index i) const
+ {
+ std::stringstream s; s << (i - bg_);
+ return s.str();
+ }
+
+ private:
+ typename Order::const_iterator bg_;
+ };
+};
+
+#if 0
+template<class StoragePolicy_ = ListRV<> >
+struct ListOrderDescriptor: public StoragePolicy_::template RebindOrder<typename std::list<ListOrderDescriptor<StoragePolicy_> >::iterator>::other
+{
+ typedef ListOrderDescriptor<StoragePolicy_> Self;
+ typedef typename std::list<Self>::iterator OrderIndex;
+ OrderIndex pair;
+
+ typedef typename StoragePolicy_::template RebindOrder<OrderIndex>::other StoragePolicy;
+ typedef typename StoragePolicy::ComplexIndex ComplexIndex;
+
+ ListOrderDescriptor(ComplexIndex i):
+ StoragePolicy(i) {}
+
+ // Acts as a rebind
+ template<class OtherStoragePolicy_> struct Order
+ { typedef std::list<ListOrderDescriptor<OtherStoragePolicy_> > type; };
+};
+
+
+template<class T>
+struct ConsistencyCmp
+{
+ int compare(T a, T b) const { if (a < b) return -1; if (a == b) return 0; return 1; }
+ bool operator()(T a, T b) const { return compare(a,b) == -1; }
+};
+
+// Traits
+template<class Order_>
+struct OrderTraits
+{
+ typedef Order_ Order;
+ typedef void StoragePolicy;
+ typedef void Index;
+ typedef void Element;
+};
+
+// Specialization
+template<class T>
+struct OrderTraits<std::vector<T> >
+{
+ typedef std::vector<T> Order;
+ typedef typename T::StoragePolicy StoragePolicy;
+ typedef typename T::OrderIndex Index;
+ typedef T Element;
+
+ typedef std::less<Index> OrderComparison;
+ typedef ConsistencyCmp<Index> ConsistencyComparison;
+
+ static Index begin(Order& order) { return order.begin(); }
+ static Index end(Order& order) { return order.end(); }
+
+ template<class Comparison>
+ static void sort(Order& order, const Comparison& cmp) { std::sort(order.begin(), order.end(), cmp); }
+};
+
+template<class T>
+struct OrderTraits<std::list<T> >
+{
+ typedef std::list<T> Order;
+ typedef typename T::StoragePolicy StoragePolicy;
+ typedef typename T::OrderIndex Index;
+ typedef T Element;
+
+ // FIXME
+ typedef std::less<Index> OrderComparison;
+ typedef ConsistencyCmp<Index> ConsistencyComparison;
+
+ static Index begin(Order& order) { return order.begin(); }
+ static Index end(Order& order) { return order.end(); }
+
+ template<class Comparison>
+ static void sort(Order& order, const Comparison& cmp) { order.sort(cmp); }
+};
+#endif
+
+#endif // __ORDER_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/persistence-diagram.h Thu Dec 18 16:43:42 2008 -0800
@@ -0,0 +1,160 @@
+#ifndef __PERSISTENCE_DIAGRAM_H__
+#define __PERSISTENCE_DIAGRAM_H__
+
+#include <utilities/types.h>
+
+#include <vector>
+#include <iostream>
+#include <cmath>
+
+#include <boost/compressed_pair.hpp>
+#include <boost/serialization/access.hpp>
+
+/**
+ * Class: PDPoint
+ *
+ * Stores birth-death pair plus any additional information provided by `Data` template parameter.
+ */
+template<class Data_>
+class PDPoint
+{
+ public:
+ typedef Data_ Data;
+
+ PDPoint(RealType x = 0, RealType y = 0, const Data& data = Data());
+
+ RealType x() const { return point_.first().first; }
+ RealType y() const { return point_.first().second; }
+ const Data& data() const { return point_.second(); }
+ Data& data() { return point_.second(); }
+
+ std::ostream& operator<<(std::ostream& out) const { return (out << x() << " " << y() << " " << data()); }
+
+ struct Visitor
+ {
+ template<class Iterator>
+ void point(Iterator i, PDPoint& p) const {}
+ };
+
+ private:
+ RealType& x() { return point_.first().first; }
+ RealType& y() { return point_.first().second; }
+
+ private:
+ boost::compressed_pair<std::pair<RealType, RealType>, Data> point_;
+
+ private:
+ /* Serialization */
+ friend class boost::serialization::access;
+
+ template<class Archive>
+ void serialize(Archive& ar, version_type );
+};
+
+template<class Data>
+std::ostream& operator<<(std::ostream& out, const PDPoint<Data>& point)
+{ return (point.operator<<(out)); }
+
+template<class Point, class Iterator, class Evaluator, class Visitor>
+Point make_point(Iterator i, const Evaluator& evaluator, const Visitor& visitor);
+
+template<class Point, class Iterator, class Evaluator>
+Point make_point(Iterator i, const Evaluator& evaluator) { return make_point<Point>(i, evaluator, Point::Visitor()); }
+
+
+/**
+ * Class: PersistenceDiagram
+ *
+ * Stores birth-death pairs, i.e. points in the extended plane. Each point can also store
+ * additional information described by `Data_` template parameter.
+ */
+template<class Data_ = Empty>
+class PersistenceDiagram
+{
+ public:
+ typedef Data_ Data;
+ typedef PDPoint<Data> Point;
+ typedef std::vector<Point> PointVector;
+ typedef typename PointVector::const_iterator const_iterator;
+
+ PersistenceDiagram() {}
+
+ template<class Iterator, class Evaluator>
+ PersistenceDiagram(Iterator bg, Iterator end,
+ const Evaluator& eval = Evaluator());
+
+ template<class Iterator, class Evaluator, class Visitor>
+ PersistenceDiagram(Iterator bg, Iterator end,
+ const Evaluator& eval = Evaluator(),
+ const Visitor& visitor = Visitor());
+
+ template<class Iterator, class Evaluator, class Visitor>
+ void init(Iterator bg, Iterator end,
+ const Evaluator& eval = Evaluator(),
+ const Visitor& visitor = Visitor());
+
+ const_iterator begin() const { return points_.begin(); }
+ const_iterator end() const { return points_.end(); }
+ size_t size() const { return points_.size(); }
+
+ void push_back(const Point& point) { points_.push_back(point); }
+
+ std::ostream& operator<<(std::ostream& out) const;
+
+ private:
+ PointVector points_;
+
+ private:
+ /* Serialization */
+ friend class boost::serialization::access;
+
+ template<class Archive>
+ void serialize(Archive& ar, version_type );
+};
+
+template<class Data>
+std::ostream& operator<<(std::ostream& out, const PersistenceDiagram<Data>& pd)
+{ return (pd.operator<<(out)); }
+
+// Function: init_diagram_vector()
+template<class Diagrams, class Iterator, class Evaluator, class DimensionExtractor>
+void init_diagrams(Diagrams& diagrams,
+ Iterator bg, Iterator end,
+ const Evaluator& evaluator = Evaluator(),
+ const DimensionExtractor& dimension = DimensionExtractor());
+
+template<class Diagrams, class Iterator, class Evaluator, class DimensionExtractor, class Visitor>
+void init_diagrams(Diagrams& diagrams,
+ Iterator bg, Iterator end,
+ const Evaluator& evaluator = Evaluator(),
+ const DimensionExtractor& dimension = DimensionExtractor(),
+ const Visitor& visitor = Visitor());
+
+// Class: Linfty
+// Functor that computes L infinity norm between two points
+template<class Point1, class Point2>
+struct Linfty
+{
+ RealType operator()(const Point1& p1, const Point2& p2) const { return std::max(std::abs(p1.x() - p2.x()),
+ std::abs(p1.y() - p2.y())); }
+
+ template<class Point>
+ RealType diagonal(const Point& p) const { return std::abs(p.y() - p.x()); }
+};
+
+// Function: bottleneck_distance(dgm1, dgm2)
+// Computes bottleneck distance between the two diagrams.
+template<class Diagram1,
+ class Diagram2,
+ class Norm>
+RealType bottleneck_distance(const Diagram1& dgm1, const Diagram2& dgm2, const Norm& norm = Norm());
+
+template<class Diagram1,
+ class Diagram2>
+RealType bottleneck_distance(const Diagram1& dgm1, const Diagram2& dgm2)
+{ return bottleneck_distance(dgm1, dgm2, Linfty<typename Diagram1::Point, typename Diagram2::Point>()); }
+
+
+#include "persistence-diagram.hpp"
+
+#endif // __PERSISTENCE_DIAGRAM_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/persistence-diagram.hpp Thu Dec 18 16:43:42 2008 -0800
@@ -0,0 +1,225 @@
+#include <boost/serialization/vector.hpp>
+#include <boost/serialization/nvp.hpp>
+
+using boost::serialization::make_nvp;
+
+template<class D>
+PDPoint<D>::
+PDPoint(RealType x, RealType y, const Data& data)
+{
+ point_.first().first = x;
+ point_.first().second = y;
+ point_.second() = data;
+}
+
+template<class D>
+template<class Iterator, class Evaluator>
+PersistenceDiagram<D>::
+PersistenceDiagram(Iterator bg, Iterator end, const Evaluator& eval)
+{
+ init(bg, end, eval, Point::Visitor());
+}
+
+template<class D>
+template<class Iterator, class Evaluator, class Visitor>
+PersistenceDiagram<D>::
+PersistenceDiagram(Iterator bg, Iterator end, const Evaluator& eval, const Visitor& visitor)
+{
+ init(bg, end, eval, visitor);
+}
+
+template<class D>
+template<class Iterator, class Evaluator, class Visitor>
+void
+PersistenceDiagram<D>::
+init(Iterator bg, Iterator end, const Evaluator& evaluator, const Visitor& visitor)
+{
+ for (Iterator cur = bg; cur != end; ++cur)
+ if (cur->sign())
+ push_back(make_point(cur, evaluator, visitor));
+}
+
+template<class Point, class Iterator, class Evaluator, class Visitor>
+Point
+make_point(Iterator i, const Evaluator& evaluator, const Visitor& visitor)
+{
+ RealType x = evaluator(i);
+ RealType y = Infinity;
+ if (i->pair != i)
+ y = evaluator(i->pair);
+
+ Point p(x,y);
+ visitor.point(i, p);
+
+ return p;
+}
+
+template<class Diagrams, class Iterator, class Evaluator, class DimensionExtractor>
+void init_diagrams(Diagrams& diagrams,
+ Iterator bg, Iterator end,
+ const Evaluator& evaluator,
+ const DimensionExtractor& dimension)
+{
+ // FIXME: this is specialized for Diagrams that is std::map
+ typedef typename Diagrams::mapped_type PDiagram;
+
+ init_diagrams(diagrams, bg, end, evaluator, dimension, typename PDiagram::Point::Visitor());
+}
+
+template<class Diagrams, class Iterator, class Evaluator, class DimensionExtractor, class Visitor>
+void init_diagrams(Diagrams& diagrams,
+ Iterator bg, Iterator end,
+ const Evaluator& evaluator,
+ const DimensionExtractor& dimension,
+ const Visitor& visitor)
+{
+ // FIXME: this is specialized for Diagrams that is std::map
+ typedef typename Diagrams::mapped_type PDiagram;
+
+ for (Iterator cur = bg; cur != end; ++cur)
+ if (cur->sign())
+ diagrams[dimension(cur)].push_back(make_point<typename PDiagram::Point>(cur, evaluator, visitor));
+}
+
+template<class D>
+std::ostream&
+PersistenceDiagram<D>::
+operator<<(std::ostream& out) const
+{
+ for (const_iterator cur = begin(); cur != end(); ++cur)
+ out << *cur << std::endl;
+ return out;
+}
+
+template<class D>
+template<class Archive>
+void
+PDPoint<D>::
+serialize(Archive& ar, version_type )
+{
+ ar & make_nvp("x", x());
+ ar & make_nvp("y", y());
+ ar & make_nvp("data", data());
+}
+
+template<class D>
+template<class Archive>
+void
+PersistenceDiagram<D>::
+serialize(Archive& ar, version_type )
+{
+ ar & make_nvp("points", points_);
+}
+
+
+/**
+ * Some structures to compute bottleneck distance between two persistence diagrams (in bottleneck_distance() function below)
+ * by setting up bipartite graphs, and finding maximum cardinality matchings in them using Boost Graph Library.
+ */
+#include <utilities/indirect.h>
+#include <boost/graph/adjacency_list.hpp>
+#include <boost/graph/max_cardinality_matching.hpp>
+
+struct Edge: public std::pair<unsigned, unsigned>
+{
+ typedef std::pair<unsigned, unsigned> Parent;
+
+ Edge(unsigned v1, unsigned v2, RealType d):
+ Parent(v1, v2), distance(d) {}
+
+ bool operator<(const Edge& other) const { return distance < other.distance; }
+
+ RealType distance;
+};
+typedef std::vector<Edge> EdgeVector;
+typedef EdgeVector::const_iterator EV_const_iterator;
+
+struct CardinaliyComparison
+{
+ typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS> Graph;
+ typedef std::vector<boost::graph_traits<Graph>::vertex_descriptor> MatchingVector;
+
+ CardinaliyComparison(unsigned size, EV_const_iterator begin):
+ max_size(size), bg(begin), last(bg), g(2*max_size), mates(2*max_size)
+ { boost::add_edge(bg->first, bg->second, g); }
+
+ bool operator()(EV_const_iterator i1, EV_const_iterator i2)
+ {
+ //std::cout << "Max size: " << max_size << std::endl;
+ //std::cout << "Comparing: (" << i1->first << ", " << i1->second << ") and "
+ // << "(" << i2->first << ", " << i2->second << ")" << std::endl;
+
+ // FIXME: the matching is being recomputed from scratch every time, this should be fixed
+ if (i2 > last)
+ do
+ {
+ ++last;
+ boost::add_edge(last->first, last->second, g);
+ } while (last != i2);
+ else
+ do
+ {
+ boost::remove_edge(last->first, last->second, g);
+ } while (--last != i2);
+
+ edmonds_maximum_cardinality_matching(g, &mates[0]);
+ //std::cout << "Found matching of size: " << matching_size(g, &mates[0]) << std::endl;
+ return matching_size(g, &mates[0]) == max_size;
+ }
+
+ unsigned max_size;
+ EV_const_iterator bg;
+ EV_const_iterator last;
+ Graph g;
+ MatchingVector mates;
+};
+
+// Bottleneck distance
+template<class Diagram1, class Diagram2, class Norm>
+RealType bottleneck_distance(const Diagram1& dgm1, const Diagram2& dgm2, const Norm& norm)
+{
+ typedef typename Diagram1::const_iterator Citer1;
+ typedef typename Diagram2::const_iterator Citer2;
+
+ const unsigned max_size = dgm1.size() + dgm2.size();
+
+ // Compute all the edges and sort them by distance
+ EdgeVector edges;
+
+ // Connect all diagonal points to each other
+ for (unsigned i = dgm1.size(); i < max_size; ++i)
+ for (unsigned j = max_size + dgm2.size(); j < 2*max_size; ++j)
+ edges.push_back(Edge(i, j, 0));
+
+ // Edges between real points
+ unsigned i = 0;
+ for (Citer1 cur1 = dgm1.begin(); cur1 != dgm1.end(); ++cur1)
+ {
+ unsigned j = max_size;
+ for (Citer2 cur2 = dgm2.begin(); cur2 != dgm2.end(); ++cur2)
+ edges.push_back(Edge(i,j++, norm(*cur1, *cur2)));
+
+ ++i;
+ }
+
+ // Edges between real points and their corresponding diagonal points
+ i = 0;
+ for (Citer1 cur1 = dgm1.begin(); cur1 != dgm1.end(); ++cur1, ++i)
+ edges.push_back(Edge(i, max_size + dgm2.size() + i, norm.diagonal(*cur1)));
+ i = max_size;
+ for (Citer2 cur2 = dgm2.begin(); cur2 != dgm2.end(); ++cur2, ++i)
+ edges.push_back(Edge(dgm1.size() + (i - max_size), i, norm.diagonal(*cur2)));
+
+
+ std::sort(edges.begin(), edges.end());
+ //for (i = 0; i < edges.size(); ++i)
+ // std::cout << "Edge: " << edges[i].first << " " << edges[i].second << " " << edges[i].distance << std::endl;
+
+ // Perform cardinality based binary search
+ RecursiveIterator<EV_const_iterator> bdistance = std::upper_bound(RecursiveIterator<EV_const_iterator>(edges.begin()),
+ RecursiveIterator<EV_const_iterator>(edges.end()),
+ edges.begin(),
+ CardinaliyComparison(max_size, edges.begin()));
+
+ return (*bdistance)->distance;
+}
--- a/include/topology/simplex.h Mon Sep 22 21:53:25 2008 -0700
+++ b/include/topology/simplex.h Thu Dec 18 16:43:42 2008 -0800
@@ -1,6 +1,6 @@
/*
* Author: Dmitriy Morozov
- * Department of Computer Science, Duke University, 2005 -- 2007
+ * Department of Computer Science, Duke University, 2005 -- 2008
*/
#ifndef __SIMPLEX_H__
@@ -13,68 +13,124 @@
#include "utilities/types.h"
+#include <boost/compressed_pair.hpp>
#include <boost/serialization/access.hpp>
/**
- * SimplexWithVertices is a basic simplex class. It stores vertices of a given type,
- * and knows how to compute its own boundary. It should probably be used as a base
- * class for any explicit simplex representation.
+ * Class: Simplex
+ * Basic simplex class. It stores vertices and data of the given types in
+ * a `boost::compressed_pair` (so if data is an Empty class which is the default,
+ * it incurs no space overhead). The class knows how to compute its own <boundary()>.
+ *
+ * Parameter:
+ * V - vertex type
+ * T - data type
*
* \ingroup topology
*/
-template<class V>
-class SimplexWithVertices
+template<class V, class T = Empty>
+class Simplex
{
public:
+ /* Typedefs: Public Types
+ *
+ * Vertex - vertex type (template parameter V)
+ * Data - data type (template parameter T)
+ * Self -
+ * Boundary - type in which the boundary is stored
+ */
typedef V Vertex;
- typedef SimplexWithVertices<Vertex> Self;
-
- typedef std::vector<Vertex> VertexContainer;
- typedef std::list<Self> Cycle;
+ typedef T Data;
+ typedef Simplex<Vertex, Data> Self;
+ typedef std::list<Self> Boundary;
+ /* Typedefs: Internal representation
+ *
+ * VertexContainer - internal representation of the vertices
+ * VerticesDataPair - `compressed_pair` of VertexContainer and Data
+ */
+ typedef std::vector<Vertex> VertexContainer;
+ typedef boost::compressed_pair<VertexContainer, Data> VerticesDataPair;
+
+ // TODO
+
/// \name Constructors
/// @{
- SimplexWithVertices() {}
- SimplexWithVertices(const Self& s):
- vertices_(s.vertices_) {}
+ //
+ /// Constructor: Simplex()
+ /// Default constructor
+ Simplex() {}
+ /// Constructor: Simplex(Self other)
+ /// Copy constructor
+ Simplex(const Self& other):
+ vdpair_(other.vertices(), other.data()) {}
+ /// Constructor: Simplex(Iterator bg, Iterator end)
+ /// Initialize simplex by copying vertices in range [bg, end)
template<class Iterator>
- SimplexWithVertices(Iterator bg, Iterator end):
- vertices_(bg, end) { std::sort(vertices_.begin(), vertices_.end()); }
- SimplexWithVertices(const VertexContainer& v):
- vertices_(v) { std::sort(vertices_.begin(), vertices_.end()); }
- SimplexWithVertices(Dimension d, Vertex v):
- vertices_() { vertices_.reserve(d+1); add(v); }
- SimplexWithVertices(Dimension d):
- vertices_(d+1) {}
+ Simplex(Iterator bg, Iterator end, const Data& d = Data()) { join(bg, end); data() = d; }
+ /// Constructor: Simplex(VertexContainer v)
+ /// Initialize simplex by copying the given VertexContainer
+ Simplex(const VertexContainer& v, const Data& d = Data()):
+ vdpair_(v, d) { std::sort(vertices().begin(), vertices().end()); }
+ /// Constructor: Simplex(Dimension d, Vertex v)
+ /// Initialize simplex of dimension d, and set its first vertex to v
+ Simplex(Dimension d, Vertex v) { vertices().reserve(d+1); add(v); }
+ /// Constructor: Simplex(Dimension d)
+ /// Initialize simplex of dimension d
+ Simplex(Dimension d) { vertices().resize(d+1); }
/// @}
/// \name Core
/// @{
- Cycle boundary() const;
- Dimension dimension() const { return vertices_.size()-1; }
+ ///
+ /// Function: boundary()
+ /// Returns the boundary of the simplex (of type Boundary)
+ Boundary boundary() const;
+ /// Function: dimension()
+ /// Returns the dimension of the simplex
+ Dimension dimension() const { return vertices().size() - 1; }
/// @}
+ const Data& data() const { return vdpair_.second(); }
+ Data& data() { return vdpair_.second(); }
+ const VertexContainer& vertices() const { return vdpair_.first(); }
+
/// \name Vertex manipulation
/// @{
- bool contains(const Vertex& v) const;
- const VertexContainer& vertices() const { return vertices_; }
+ bool contains(const Vertex& v) const;
void add(const Vertex& v);
- void join(const Self& other);
+ template<class Iterator>
+ void join(Iterator bg, Iterator end);
+ void join(const Self& other) { join(other.vertices.begin(), other.vertices().end()); }
/// @}
- /// \name Assignment and comparison
- /// Gives an ordering on simplices (for example, so that simplices can be used as keys for std::map)
- /// @{
- const Self& operator=(const Self& s) { vertices_ = s.vertices_; return *this; }
- bool operator==(const Self& s) const { return vertices_ == s.vertices_; }
- bool operator<(const Self& s) const { return vertices_ < s.vertices_; }
- /// @}
+ const Self& operator=(const Self& s) { vdpair_ = s.vdpair_; return *this; }
std::ostream& operator<<(std::ostream& out) const;
+
+ /* Classes: Comparisons
+ *
+ * VertexComparison - compare simplices based on the lexicographic ordering of their <vertices()>
+ * DataComparison - compare simplices based on their <data()>
+ * DataDimensionComparison - compare simplices based on their <data()> within each <dimension()>
+ */
+
+ /* Classes: Functors
+ * DataEvaluator - return data given a simplex
+ * DimensionExtractor - return dimesnion given a simplex
+ */
+ struct VertexComparison;
+ struct DataComparison;
+ struct DataDimensionComparison;
+
+ struct DataEvaluator;
+ struct DimensionExtractor;
private:
- VertexContainer vertices_;
+ VertexContainer& vertices() { return vdpair_.first(); }
+
+ VerticesDataPair vdpair_;
private:
/* Serialization */
@@ -84,109 +140,67 @@
void serialize(Archive& ar, version_type );
};
-/**
- * SimplexWithValue explicitly adds a RealType value to the SimplexWithVertices.
- *
- * \ingroup topology
- */
-template<class Vert>
-class SimplexWithValue: public SimplexWithVertices<Vert>
+
+template<class V, class T>
+struct Simplex<V,T>::VertexComparison
{
- public:
- typedef Vert Vertex;
- typedef RealType Value;
- typedef SimplexWithValue<Vertex> Self;
- typedef SimplexWithVertices<Vertex> Parent;
+ typedef Self first_argument_type;
+ typedef Self second_argument_type;
+ typedef bool result_type;
- typedef typename Parent::VertexContainer VertexContainer;
-
- /// \name Constructors
- /// @{
- SimplexWithValue(Value value = 0): val(value) {}
- SimplexWithValue(const Self& s):
- Parent(s), val(s.val) {}
- SimplexWithValue(const Parent& s, Value value = 0):
- Parent(s), val(value) {}
- template<class Iterator>
- SimplexWithValue(Iterator bg, Iterator end, Value value = 0):
- Parent(bg, end), val(value) {}
- SimplexWithValue(const VertexContainer& v, Value value = 0):
- Parent(v), val(value) {}
- /// @}
+ bool operator()(const Self& a, const Self& b) const { return a.vertices() < b.vertices(); }
+};
- /// \name Core
- /// @{
- void set_value(Value value) { val = value; }
- Value get_value() const { return val; }
- /// @}
-
- const Self& operator=(const Self& s);
- std::ostream& operator<<(std::ostream& out) const;
+template<class V, class T>
+struct Simplex<V,T>::DataComparison
+{
+ typedef Self first_argument_type;
+ typedef Self second_argument_type;
+ typedef bool result_type;
- private:
- Value val;
-
- /* Serialization */
- friend class boost::serialization::access;
-
- template<class Archive>
- void serialize(Archive& ar, version_type );
+ bool operator()(const Self& a, const Self& b) const { return a.data() < b.data(); }
};
-/**
- * SimplexWithAttachment stores the vertex to which the simplex is attached (meant for lower-star filtrations)
- *
- * \ingroup topology
- */
-template<typename V>
-class SimplexWithAttachment: public SimplexWithVertices<V>
+template<class V, class T>
+struct Simplex<V,T>::DataDimensionComparison
{
- public:
- typedef V VertexIndex;
- typedef SimplexWithVertices<VertexIndex> Parent;
-
- /// \name Constructors
- /// @{
- SimplexWithAttachment():
- attachment(VertexIndex()) {}
- template<class Iterator>
- SimplexWithAttachment(Iterator bg, Iterator end):
- Parent(bg, end) {}
- SimplexWithAttachment(const Parent& s):
- Parent(s) {}
- SimplexWithAttachment(Dimension d, VertexIndex vi):
- Parent(d, vi), attachment(vi) {}
- /// @}
+ typedef Self first_argument_type;
+ typedef Self second_argument_type;
+ typedef bool result_type;
- void set_attachment(VertexIndex v) { attachment = v; }
- VertexIndex get_attachment() const { return attachment; }
-
- private:
- VertexIndex attachment;
-
- private:
- // Serialization
- friend class boost::serialization::access;
+ bool operator()(const Self& a, const Self& b) const
+ {
+ if (a.dimension() == b.dimension())
+ return a.data() < b.data();
+ else
+ return a.dimension() < b.dimension();
+ }
+};
+
+template<class V, class T>
+struct Simplex<V,T>::DataEvaluator
+{
+ typedef Self first_argument_type;
+ typedef Data result_type;
- template<class Archive>
- void serialize(Archive& ar, version_type );
+ result_type operator()(const first_argument_type& s) const { return s.data(); }
};
-template<class Simplex_>
-class DimensionValueComparison
+template<class V, class T>
+struct Simplex<V,T>::DimensionExtractor
{
- public:
- typedef Simplex_ Simplex;
+ typedef Self first_argument_type;
+ typedef Dimension result_type;
- bool operator()(const Simplex& s1, const Simplex& s2) const
- {
- if (s1.dimension() == s2.dimension())
- return s1.get_value() < s2.get_value();
- else
- return s1.dimension() < s2.dimension();
- }
+ result_type operator()(const first_argument_type& s) const { return s.dimension(); }
};
+
+// TODO: class DirectSimplex - class which stores indices of the simplices in its boundary
+// TODO: class CompactSimplex<V, T, N> - uses arrays instead of vectors to store simplices
+// (dimension N must be known at compile time)
+
+
#include "simplex.hpp"
#endif // __SIMPLEX_H__
--- a/include/topology/simplex.hpp Mon Sep 22 21:53:25 2008 -0700
+++ b/include/topology/simplex.hpp Thu Dec 18 16:43:42 2008 -0800
@@ -2,117 +2,80 @@
#include <boost/serialization/set.hpp>
#include <boost/serialization/nvp.hpp>
+#include <boost/iterator/filter_iterator.hpp>
+#include <functional>
/* Implementations */
-/* SimplexWithVertices */
-template<class V>
-typename SimplexWithVertices<V>::Cycle
-SimplexWithVertices<V>::
+/* Simplex */
+template<class V, class T>
+typename Simplex<V,T>::Boundary
+Simplex<V,T>::
boundary() const
{
- Cycle bdry;
- if (dimension() == 0) return bdry;
+ typedef std::not_equal_to<Vertex> NotEqualVertex;
- for (typename VertexContainer::const_iterator cur = vertices_.begin(); cur != vertices_.end(); ++cur)
- {
- bdry.push_back(Self(dimension() - 1));
- Self& s = bdry.back();
- std::remove_copy(vertices_.begin(), vertices_.end(), s.vertices_.begin(), *cur);
- }
-
- return bdry;
+ Boundary bdry;
+ if (dimension() == 0) return bdry;
+
+ for (typename VertexContainer::const_iterator cur = vertices().begin(); cur != vertices().end(); ++cur)
+ bdry.push_back(Self(boost::make_filter_iterator(std::bind2nd(NotEqualVertex(), *cur), vertices().begin(), vertices().end()),
+ boost::make_filter_iterator(std::bind2nd(NotEqualVertex(), *cur), vertices().end(), vertices().end())));
+
+ return bdry;
}
-template<class V>
+template<class V, class T>
bool
-SimplexWithVertices<V>::
+Simplex<V,T>::
contains(const Vertex& v) const
{
- typename VertexContainer::const_iterator location = std::lower_bound(vertices_.begin(), vertices_.end(), v);
- return ((location == vertices_.end()) || (*location == v));
+ // TODO: would std::find() be faster? (since most simplices we deal with are low dimensional)
+ typename VertexContainer::const_iterator location = std::lower_bound(vertices().begin(), vertices().end(), v);
+ return ((location != vertices().end()) && (*location == v));
}
-template<class V>
+template<class V, class T>
void
-SimplexWithVertices<V>::
+Simplex<V,T>::
add(const Vertex& v)
-{ vertices_.push_back(v); std::sort(vertices_.begin(), vertices_.end()); }
+{
+ // TODO: would find() or lower_bound() followed by insert be faster?
+ vertices().push_back(v); std::sort(vertices().begin(), vertices().end());
+}
-template<class V>
+template<class V, class T>
+template<class Iterator>
void
-SimplexWithVertices<V>::
-join(const Self& other)
+Simplex<V,T>::
+join(Iterator bg, Iterator end)
{
- vertices_.insert(vertices_.end(), other.vertices_.begin(), other.vertices_.end());
- std::sort(vertices_.begin(), vertices_.end());
+ vertices().insert(vertices().end(), bg, end);
+ std::sort(vertices().begin(), vertices().end());
}
-template<class V>
+template<class V, class T>
std::ostream&
-SimplexWithVertices<V>::
+Simplex<V,T>::
operator<<(std::ostream& out) const
{
- for (typename VertexContainer::const_iterator cur = vertices_.begin(); cur != vertices_.end(); ++cur)
- out << *cur << ' ';
-
+ for (typename VertexContainer::const_iterator cur = vertices().begin(); cur != vertices().end(); ++cur)
+ out << *cur;
+ out << " [" << data() << "] ";
+
return out;
}
-template<class V>
+template<class V, class T>
template<class Archive>
void
-SimplexWithVertices<V>::
+Simplex<V,T>::
serialize(Archive& ar, version_type )
-{ ar & BOOST_SERIALIZATION_NVP(vertices_); }
-
-template<class V>
-std::ostream& operator<<(std::ostream& out, const SimplexWithVertices<V>& s)
-{ return s.operator<<(out); }
-
-
-/* SimplexWithValue */
-template<class V>
-std::ostream&
-SimplexWithValue<V>::
-operator<<(std::ostream& out) const
-{
- Parent::operator<<(out);
- out << "(val = " << val << ")";
- return out;
+{
+ ar & make_nvp("vertices", vertices());
+ ar & make_nvp("data", data());
}
-template<class V>
-const typename SimplexWithValue<V>::Self&
-SimplexWithValue<V>::
-operator=(const Self& s)
-{
- Parent::operator=(s);
- val = s.val;
- return *this;
-}
-
-template<class V>
-template<class Archive>
-void
-SimplexWithValue<V>::
-serialize(Archive& ar, version_type )
-{
- ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP(Parent);
- ar & BOOST_SERIALIZATION_NVP(val);
-}
-
-template<typename V>
-template<class Archive>
-void
-SimplexWithAttachment<V>::
-serialize(Archive& ar, version_type )
-{
- ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP(Parent);
- ar & BOOST_SERIALIZATION_NVP(attachment);
-}
-
-
-template<class V>
-std::ostream& operator<<(std::ostream& out, const SimplexWithValue<V>& s)
+template<class V, class T>
+std::ostream& operator<<(std::ostream& out, const Simplex<V,T>& s)
{ return s.operator<<(out); }
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/static-persistence.h Thu Dec 18 16:43:42 2008 -0800
@@ -0,0 +1,101 @@
+#ifndef __STATIC_PERSISTENCE_H__
+#define __STATIC_PERSISTENCE_H__
+
+#include "order.h"
+#include "filtration.h"
+
+/**
+ * Class: StaticPersistence
+ * The class that encapsulates order and pairing information as well as
+ * implements methods to compute and maintain the pairing. Its most
+ * significant method is <pair_simplices()>.
+ *
+ * Template parameters:
+ * Data_ - auxilliary contents to store with each OrderElement
+ * OrderDescriptor_ - class describing how the order is stored; it defaults to <VectorOrderDescriptor>
+ * which serves as a prototypical class
+ */
+template<class Data_ = Empty,
+ class OrderDescriptor_ = VectorOrderDescriptor<> >
+class StaticPersistence
+{
+ public:
+ // Typedef: Data
+ // The data type stored in each order element
+ typedef Data_ Data;
+
+ /* Typedefs: Order
+ * Types associated with the order
+ *
+ * OrderDescriptor - template parameter rebound with `Data` carrying all the necessary information
+ * Order - the container class itself
+ * OrderIndex - iterator into the Order
+ * OrderComparison - a class for comparing OrderIndices
+ */
+ typedef typename OrderDescriptor_::
+ template RebindData<Data>::
+ other OrderDescriptor;
+ typedef OrderTraits<OrderDescriptor> Traits;
+ typedef typename Traits::Order Order;
+ typedef typename Traits::Index OrderIndex;
+ typedef typename Traits::Element OrderElement;
+ typedef typename Traits::Comparison OrderComparison;
+
+ /* Constructor: StaticPersistence()
+ * TODO: write a description
+ *
+ * Template parameters:
+ * Filtration - filtration of the complex whose persistence we are computing
+ */
+ template<class Filtration> StaticPersistence(const Filtration& f,
+ const OrderComparison& ocmp = OrderComparison());
+
+ // Function: pair_simplices()
+ // Compute persistence of the filtration
+ void pair_simplices() { pair_simplices<PairVisitor>(begin(), end()); }
+
+ // Functions: Accessors
+ // 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(); }
+ size_t size() const { return order_.size(); }
+ const OrderComparison& order_comparison() const { return ocmp_; }
+
+ 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());
+
+ // Struct: PairVisitor
+ // Acts as an archetype and if necessary a base class for visitors passed to <pair_simplices(bg, end, visitor)>.
+ struct PairVisitor
+ {
+ // 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 {}
+
+ // 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 {}
+
+ // Function: finished(j)
+ // Called after the processing of `j` is finished.
+ void finished(OrderIndex j) const {}
+ };
+
+ const Order& order() const { return order_; }
+
+ private:
+ Order order_;
+ OrderComparison ocmp_;
+};
+
+#include "static-persistence.hpp"
+
+#endif // __STATIC_PERSISTENCE_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/static-persistence.hpp Thu Dec 18 16:43:42 2008 -0800
@@ -0,0 +1,99 @@
+#include <utilities/log.h>
+#include <utilities/containers.h>
+#include <boost/utility/enable_if.hpp>
+#include <boost/utility.hpp>
+#include <utilities/property-maps.h>
+
+#ifdef LOGGING
+static rlog::RLogChannel* rlPersistence = DEF_CHANNEL("topology/persistence", rlog::Log_Debug);
+#endif // LOGGING
+
+#ifdef COUNTERS
+static Counter* cPersistencePair = GetCounter("persistence/pair");
+static Counter* cPersistencePairBoundaries = GetCounter("persistence/pair/boundaries");
+static Counter* cPersistencePairCycleLength = GetCounter("persistence/pair/cyclelength");
+#endif // COUNTERS
+
+template<class D, class OD>
+template<class Filtration>
+StaticPersistence<D, OD>::
+StaticPersistence(const Filtration& filtration,
+ const OrderComparison& ocmp):
+ order_(filtration.size()),
+ ocmp_(ocmp)
+{
+ OrderIndex ocur = begin();
+ OffsetMap<size_t, OrderIndex> om(0, ocur); // TODO: this is customized for std::vector Order
+ for (typename Filtration::Index cur = filtration.begin(); cur != filtration.end(); ++cur, ++ocur)
+ {
+ // Convert the Filtration::IndexBoundary into a Cycle, and
+ typedef typename OrderDescriptor::Chains::Cycle Cycle;
+
+ OrderElement& oe = *ocur;
+ oe.pair = ocur;
+
+ filtration.boundary(cur, oe.cycle, om);
+ oe.cycle.sort(ocmp_);
+ }
+}
+
+template<class D, class OD>
+template<class Visitor>
+void
+StaticPersistence<D, OD>::
+pair_simplices(OrderIndex bg, OrderIndex end, const Visitor& visitor)
+{
+#if LOGGING
+ typename Traits::OutputMap outmap(order_);
+#endif
+
+ // FIXME: need sane output for logging
+ rLog(rlPersistence, "Entered: pair_simplices");
+ for (OrderIndex j = bg; j != end; ++j)
+ {
+ visitor.init(j);
+ rLog(rlPersistence, "Simplex %s", outmap(j).c_str());
+
+ OrderElement& oe = *j;
+ typedef typename OrderDescriptor::Chains::Cycle Cycle;
+ Cycle& z = oe.cycle;
+ rLog(rlPersistence, " has boundary: %s", z.tostring(outmap).c_str());
+
+ CountNum(cPersistencePairBoundaries, oe.cycle.size());
+ Count(cPersistencePair);
+
+ while(!z.empty())
+ {
+ 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),
+ "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)
+ {
+ 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());
+ break;
+ }
+
+ // update element
+ z.add(i->pair->cycle, ocmp_);
+ visitor.update(j, i);
+ rLog(rlPersistence, " new cycle: %s", z.tostring(outmap).c_str());
+ }
+ visitor.finished(j);
+ rLog(rlPersistence, "Finished with %s: %s",
+ outmap(j).c_str(), outmap(oe.pair).c_str());
+ }
+}
+
+
+/* Private */
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/utilities/containers.h Thu Dec 18 16:43:42 2008 -0800
@@ -0,0 +1,89 @@
+#ifndef __CONTAINERS_H__
+#define __CONTAINERS_H__
+
+#include <vector>
+#include "circular_list.h"
+
+// TODO: write documentation
+
+template<class Container_, class Comparison_ = std::less<typename Container_::value_type> >
+struct ContainerTraits
+{
+ typedef Container_ Container;
+ typedef typename Container::value_type Item;
+ typedef Comparison_ Comparison;
+
+ static void reserve(Container& c, size_t sz) {}
+ static void sort(Container& c, const Comparison& cmp = Comparison()) { c.sort(cmp); }
+};
+
+/**
+ * Class: SizeStorage
+ *
+ * This class expresses how size should be stored for various containers to be able
+ * to deduce it in constant time. By default it is stored explicitly, so that if the
+ * Container is std::list everything works. However, specialization is available for
+ * std::vector which uses its builtin size() function.
+ */
+template<class Container_>
+class SizeStorage
+{
+ public:
+ typedef Container_ Container;
+
+ SizeStorage(): size_(0) {}
+
+ void increment(size_t inc = 1) { size_ += inc; }
+ void decrement(size_t dec = 1) { size_ -= dec; }
+ size_t size(const Container& c) const { return size_; }
+ void swap(SizeStorage& other) { std::swap(size_, other.size_); }
+
+ private:
+ size_t size_;
+};
+
+
+/* Specializations */
+
+template<class T, class Comparison_>
+struct ContainerTraits<std::vector<T>, Comparison_>
+{
+ typedef T Item;
+ typedef std::vector<T> Container;
+ typedef Comparison_ Comparison;
+
+ static void reserve(Container& c, size_t sz) { c.reserve(sz); }
+ static void sort(Container& c, const Comparison& cmp = Comparison()) { std::sort(c.begin(), c.end(), cmp); }
+};
+
+template<class T, class Comparison_>
+struct ContainerTraits<List<T>, Comparison_>
+{
+ typedef T Item;
+ typedef List<T> Container;
+ typedef Comparison_ Comparison;
+
+ static void reserve(Container& c, size_t sz) { }
+ static void sort(Container& c, const Comparison& cmp = Comparison())
+ {
+ std::vector<Item> tmp(c.begin(), c.end());
+ std::sort(tmp.begin(), tmp.end(), cmp);
+ std::copy(tmp.begin(), tmp.end(), c.begin());
+ }
+};
+
+// TODO: specialize for List (singly-linked list)
+
+template<class T>
+class SizeStorage<std::vector<T> >
+{
+ public:
+ typedef std::vector<T> Container;
+
+ void increment(size_t inc = 1) {}
+ void decrement(size_t dec = 1) {}
+ size_t size(const Container& c) const { return c.size(); }
+ void swap(SizeStorage& other) {}
+};
+
+#endif // __CONTAINERS_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/utilities/indirect.h Thu Dec 18 16:43:42 2008 -0800
@@ -0,0 +1,66 @@
+#ifndef __INDIRECT_H__
+#define __INDIRECT_H__
+
+#include <boost/iterator/iterator_adaptor.hpp>
+
+// TODO: write documentation
+
+
+template<class Iterator_, class Comparison_>
+struct IndirectComparison
+{
+ typedef Iterator_ Iterator;
+ typedef Comparison_ Comparison;
+
+ IndirectComparison(const Comparison& cmp):
+ cmp_(cmp)
+ {}
+
+ bool operator()(Iterator a, Iterator b) const
+ { return cmp_(*a, *b); }
+
+ const Comparison& cmp_;
+};
+
+template<class Comparison>
+struct ThreeOutcomeCompare: public Comparison
+{
+ typedef typename Comparison::first_argument_type first_argument_type;
+ typedef typename Comparison::second_argument_type second_argument_type;
+
+ ThreeOutcomeCompare(const Comparison& cmp = Comparison()): Comparison(cmp) {}
+
+ int compare(const first_argument_type& a, const second_argument_type& b) const
+ { if (operator()(a,b)) return -1;
+ else if (operator()(b,a)) return 1;
+ else return 0;
+ }
+};
+
+template<class Iterator_>
+class RecursiveIterator: public boost::iterator_adaptor<RecursiveIterator<Iterator_>, // Derived
+ Iterator_, // Base
+ Iterator_> // Value
+{
+ private:
+ struct enabler {};
+
+ public:
+ typedef Iterator_ Iterator;
+ typedef boost::iterator_adaptor<RecursiveIterator<Iterator>,
+ Iterator,
+ Iterator> Parent;
+
+ RecursiveIterator() {}
+ explicit RecursiveIterator(Iterator iter):
+ Parent(iter) {}
+
+ private:
+ friend class boost::iterator_core_access;
+ typename Parent::reference
+ dereference() const { return const_cast<typename Parent::reference>(this->base()); }
+ // FIXME: I dislike to const_cast, but it's not obvious how to get rid of it
+};
+
+
+#endif // __INDIRECT_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/utilities/property-maps.h Thu Dec 18 16:43:42 2008 -0800
@@ -0,0 +1,137 @@
+#ifndef __PROPERTY_MAPS_H__
+#define __PROPERTY_MAPS_H__
+
+#include <boost/property_map.hpp>
+#include <boost/iterator/iterator_traits.hpp>
+#include <algorithm>
+#include "utilities/log.h"
+
+
+/* Associative Map */ // FIXME: this needs to have more thought put into it
+template<class UniquePairAssociativeContainer_>
+class AssociativeMap: public boost::associative_property_map<UniquePairAssociativeContainer_>
+{
+ public:
+ typedef boost::associative_property_map<UniquePairAssociativeContainer_> Parent;
+
+ public:
+ // FIXME: takes begin, end, and initializes with the reverse
+ AssociativeMap():
+ container(), Parent(container) {}
+
+ private:
+ UniquePairAssociativeContainer_ container;
+};
+
+
+/* Dereference Map */
+template<class Iterator_>
+class DereferenceMap
+{
+ public:
+ typedef Iterator_ Iterator;
+ typedef Iterator key_type;
+ typedef typename boost::iterator_value<Iterator>::type value_type;
+ typedef boost::readable_property_map_tag category;
+
+ public:
+ value_type operator[](const key_type& k) const { return *k; }
+};
+
+
+/* Binary Search Map */
+template<class Query_, class Index_, class Comparison_ = std::less<Query_> >
+class BinarySearchMap
+{
+ public:
+ typedef Query_ Query;
+ typedef Index_ Index;
+ typedef Comparison_ Comparison;
+
+ typedef Query key_type;
+ typedef Index value_type;
+ typedef boost::readable_property_map_tag category;
+
+ public:
+ BinarySearchMap(Index bg, Index end,
+ Comparison cmp = Comparison()):
+ bg_(bg), end_(end), cmp_(cmp) {}
+
+ value_type operator[](const key_type& k) const
+ {
+ value_type res = std::lower_bound(bg_, end_, k, cmp_);
+ AssertMsg(!cmp_(*res, k) && !cmp_(k, *res), "Query must always be found");
+ return res;
+ }
+
+ private:
+ Index bg_;
+ Index end_;
+ Comparison cmp_;
+
+};
+
+/* Offset Map */
+template<class From_, class To_>
+struct OffsetMap
+{
+ typedef From_ From;
+ typedef To_ To;
+ typedef From key_type;
+ typedef To value_type;
+
+ OffsetMap(From bg_from, To bg_to):
+ bg_from_(bg_from), bg_to_(bg_to) {}
+
+ To operator[](From i) const { return bg_to_ + (i - bg_from_); }
+
+ From from() const { return bg_from_; }
+ To to() const { return bg_to_; }
+
+
+ template<class NewFrom_> struct rebind_from
+ { typedef OffsetMap<NewFrom_, To_> other; };
+ template<class NewTo_> struct rebind_to
+ { typedef OffsetMap<From_, NewTo_> other; };
+
+
+ private:
+ From bg_from_;
+ To bg_to_;
+};
+
+template<class From_, class To_>
+OffsetMap<From_, To_>
+make_offset_map(From_ bg_from, To_ bg_to)
+{ return OffsetMap<From_, To_>(bg_from, bg_to); }
+
+
+/* ThroughMap */
+template<class Functor_, class Map_>
+class ThroughMap
+{
+ public:
+ typedef Map_ Map;
+ typedef Functor_ Functor;
+
+ typedef typename Functor::result_type result_type;
+ typedef typename Map::key_type first_argument_type;
+
+ ThroughMap(const Map& map,
+ const Functor& functor):
+ map_(map),
+ functor_(functor) {}
+
+ result_type operator()(first_argument_type a) const { return functor_(map_[a]); }
+
+ private:
+ const Map& map_;
+ const Functor& functor_;
+};
+
+template<class Map, class Functor>
+ThroughMap<Functor, Map>
+evaluate_through_map(const Map& map, const Functor& functor)
+{ return ThroughMap<Functor, Map>(map, functor); }
+
+#endif // __PROPERTY_MAPS_H__
--- a/include/utilities/types.h Mon Sep 22 21:53:25 2008 -0700
+++ b/include/utilities/types.h Thu Dec 18 16:43:42 2008 -0800
@@ -2,6 +2,7 @@
#define __TYPES_H__
#include <limits>
+#include <iostream>
/* Types */
typedef bool Sign;
@@ -15,4 +16,33 @@
typedef const unsigned int& version_type;
+struct Empty {};
+std::ostream& operator<<(std::ostream& out, Empty e) { return out; }
+
+enum SwitchType
+{
+ DiffDim = 0,
+ Case1 = 0x4,
+ Case12 = 0x5,
+ Case112 = 0x6,
+ Case2 = 0x8,
+ Case212 = 0x9,
+ Case3 = 0x10,
+ Case31 = 0x11,
+ Case4 = 0x20,
+};
+
+
+// Nothing to do for serializing Empty, but still need to provide this function
+namespace boost {
+namespace serialization {
+
+template<class Archive>
+void serialize(Archive & ar, Empty&, const unsigned int )
+{}
+
+} // namespace serialization
+} // namespace boost
+
+
#endif // __TYPES_H__