Intermediate commit while converting to the new architecture dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Thu, 18 Dec 2008 16:43:42 -0800
branchdev
changeset 97 0a9bd3f34419
parent 96 f283106e8124
child 98 d81e460e267a
child 101 9efac05d629b
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
CMakeLists.txt
examples/CMakeLists.txt
examples/alphashapes/CMakeLists.txt
examples/alphashapes/alphashapes3d.cpp
examples/alphashapes/alphashapes3d.h
examples/alphashapes/alphashapes3d.hpp
examples/alphashapes/compare-diagrams.cpp
examples/ar-vineyard/CMakeLists.txt
examples/cech-complex/cech-complex.cpp
examples/fitness/CMakeLists.txt
examples/poincare/CMakeLists.txt
examples/triangle/CMakeLists.txt
examples/triangle/triangle.cpp
include/topology/chain.h
include/topology/chain.hpp
include/topology/complex-traits.h
include/topology/cycle.h
include/topology/cycle.hpp
include/topology/cycles.h
include/topology/dynamic-persistence.h
include/topology/dynamic-persistence.hpp
include/topology/filtration.h
include/topology/filtration.hpp
include/topology/order.h
include/topology/persistence-diagram.h
include/topology/persistence-diagram.hpp
include/topology/simplex.h
include/topology/simplex.hpp
include/topology/static-persistence.h
include/topology/static-persistence.hpp
include/utilities/containers.h
include/utilities/indirect.h
include/utilities/property-maps.h
include/utilities/types.h
--- 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__