Merged in changes from master ar
authorDmitriy Morozov <morozov@cs.duke.edu>
Tue, 19 Feb 2008 06:56:49 -0500
branchar
changeset 60 bb79215d1f93
parent 59 4eb311c4d0e7 (current diff)
parent 40 122e4a1fa117 (diff)
child 61 f905b57dd7ab
Merged in changes from master
CMakeLists.txt
examples/ar-vineyard/CMakeLists.txt
examples/ar-vineyard/ar-vineyard.h
examples/ar-vineyard/ar-vineyard.hpp
include/topology/filtration.h
include/utilities/debug.h
include/utilities/sys.h
src/debug.cpp
--- a/CMakeLists.txt	Tue Feb 19 06:51:23 2008 -0500
+++ b/CMakeLists.txt	Tue Feb 19 06:56:49 2008 -0500
@@ -1,7 +1,8 @@
 project						(Dionysus)
 
+option						(logging 			"Build Dionysus with logging on" 		OFF)
+option						(counters			"Build Dionysus with counters on" 		OFF)
 option						(debug				"Build Dionysus with debugging on" 		OFF)
-option						(counters			"Build Dionysus with counters on" 		OFF)
 option						(optimize			"Build Dionysus with optimization"		ON)
 
 # Find everything that's always required
@@ -10,9 +11,6 @@
 find_library				(dsrpdb_LIBRARY 			NAMES dsrpdb)
 find_path					(dsrpdb_INCLUDE_DIR 		dsrpdb/Protein.h)
 
-set							(libraries 					${libraries}
-														${dsrpdb_LIBRARY})
-
 #CGAL
 execute_process				(COMMAND ${CMAKE_MAKE_PROGRAM} -f ${CMAKE_CURRENT_SOURCE_DIR}/FindCGAL.Makefile libpaths
 							 OUTPUT_VARIABLE cgal_libpaths)
@@ -26,7 +24,6 @@
 #string						(REPLACE "\n" "" cgal_ldflags 	${cgal_ldflags})
 string						(REPLACE "\n" "" cgal_cxxflags 	${cgal_cxxflags})
 string						(REPLACE "\n" "" cgal_libpath 	${cgal_libpath})
-add_definitions				(${cgal_cxxflags})
 find_library				(cgal_LIBRARY				NAMES CGAL
 														PATHS ${cgal_libpath})
 find_library				(core_LIBRARY				NAMES CGALcore++
@@ -52,36 +49,42 @@
 
 # Debugging
 if							(debug)
-	find_library			(cwd_LIBRARY				NAMES cwd)
-	find_path				(cwd_INCLUDE_DIR			libcwd/debug.h)
-	set 					(cwd_INCLUDE_DIR			${cwd_INCLUDE_DIR}/libcwd)
-	add_definitions			(-DCWDEBUG -g)
-	set						(external_sources			${CMAKE_CURRENT_SOURCE_DIR}/src/debug.cpp)
-	set						(libraries 					${libraries} ${cwd_LIBRARY})
+	if 						(optimize)
+			set				(cxx_flags					${CMAKE_CXX_FLAGS_RELWITHDEBINFO})
+	else					(optimize)
+			set				(cxx_flags					${CMAKE_CXX_FLAGS_DEBUG})
+	endif					(optimize)
 else						(debug)
-	add_definitions			(-DNDEBUG)
+	if 						(optimize)
+			set				(cxx_flags					${CMAKE_CXX_FLAGS_RELEASE})
+	else					(optimize)
+			set				(cxx_flags					${CMAKE_CXX_FLAGS})
+	endif					(optimize)
 endif						(debug)
+add_definitions				(${cxx_flags})
+
+# Logging
+if 							(logging)
+	find_library			(rlog_LIBRARY				NAMES rlog)
+	find_path				(rlog_INCLUDE_DIR			rlog/rlog.h)
+	set						(rlog_INCLUDE_DIR			${rlog_INCLUDE_DIR})
+	add_definitions			(-DLOGGING -DRLOG_COMPONENT=dionysus)
+	set						(libraries					${libraries} ${rlog_LIBRARY})
+endif						(logging)
 
 # Counters
 if							(counters)
 	add_definitions			(-DCOUNTERS)
 endif						(counters)
 
-# Optimization
-if							(optimize GREATER 0)
-	add_definitions			(-O${optimize})
-elseif						(optimize)
-	add_definitions			(-O)
-endif						(optimize)
-endif						(optimize GREATER 0)
-
 
 # Set includes
 include_directories			(${CMAKE_CURRENT_BINARY_DIR}
 							 ${CMAKE_CURRENT_SOURCE_DIR}/include
 							 ${Boost_INCLUDE_DIR}
 							 ${dsrpdb_INCLUDE_DIR}
-							 ${cwd_INCLUDE_DIR})
+							 ${cwd_INCLUDE_DIR}
+							 ${rlog_INCLUDE_DIR})
 
 # Doxygen (FIXME)
 if							(DOXYGEN_FOUND)
@@ -90,6 +93,6 @@
 #							DEPENDS Doxyfile)
 endif						(DOXYGEN_FOUND)
 
-# Set external sources
+# Process subdirectories
 add_subdirectory			(examples)
 add_subdirectory			(tests)
--- a/README	Tue Feb 19 06:51:23 2008 -0500
+++ b/README	Tue Feb 19 06:56:49 2008 -0500
@@ -1,15 +1,15 @@
 Dependencies
   CGAL-3.3 -    for alpha-shapes and kinetic data structures
   DSR-PDB -     for reading in PDB files
-  cmake -       for controlling the build process
+  CMake -       for controlling the build process
   boost -       great set of C++ libraries
   Doxygen -     for building documentation
-  libcwd -      for debugging only (is not needed by default)
-  synaps -      for solving polynomial (for kinetic kernel), which in turn requires GMP
+  rlog -        for logging only (is not needed by default)
+  SYNAPS -      for solving polynomials (for kinetic kernel), which in turn requires GMP
 
 Configuration
   The path to CGAL's Makefile is expected to be set in $CGAL_MAKEFILE, the rest
-  is just usual CMake configuration
+  is just the usual CMake configuration
 
 Building
   To build examples, create a directory build (to keep everything in one place),
@@ -20,9 +20,10 @@
   make
   
   In the cmake line you can provide -Ddebug:bool=on to turn on debugging,
-  -Dcounters:bool=on to turn on counters, -Doptimize:int=3 would set
-  optimization to -O3. All of this can be set using a text user interface by
-  running ccmake instead of cmake.
+  -Dcounters:bool=on to turn on counters, -Doptimize:bool=on to turn on
+  optimization.  Depending on the combination of debugging and optimization, a
+  particular CMAKE_CXX_FLAGS* is chosen.  All of this can be set using a text
+  user interface by running ccmake instead of cmake.
 
 Author
   Dmitriy Morozov <morozov@cs.duke.edu>
--- a/examples/alphashapes/CMakeLists.txt	Tue Feb 19 06:51:23 2008 -0500
+++ b/examples/alphashapes/CMakeLists.txt	Tue Feb 19 06:56:49 2008 -0500
@@ -1,8 +1,10 @@
 set							(targets						
 							 alphashapes3d
+							 alphashapes2d
 							 alpharadius)
 							 
+add_definitions				(${cgal_cxxflags})
 foreach 					(t ${targets})
-	add_executable			(${t} ${t}.cpp ${external_sources})
+	add_executable			(${t} ${t}.cpp)
 	target_link_libraries	(${t} ${libraries} ${cgal_libraries})
 endforeach 					(t ${targets})
--- a/examples/alphashapes/alpharadius.cpp	Tue Feb 19 06:51:23 2008 -0500
+++ b/examples/alphashapes/alpharadius.cpp	Tue Feb 19 06:56:49 2008 -0500
@@ -1,5 +1,4 @@
-#include <utilities/sys.h>
-#include <utilities/debug.h>
+#include <utilities/log.h>
 
 #include "alphashapes3d.h"
 #include <topology/filtration.h>
@@ -48,11 +47,13 @@
 
 int main(int argc, char** argv) 
 {
-#ifdef CWDEBUG
-	Debug(dc::filtration.on());
-	//Debug(dc::cycle.on());
+#ifdef LOGGING
+	rlog::RLogInit(argc, argv);
 
-	dionysus::debug::init();
+	stdoutLog.subscribeTo( RLOG_CHANNEL("info") );
+	stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
+	stdoutLog.subscribeTo( RLOG_CHANNEL("topology/filtration") );
+	//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/cycle") );
 #endif
 
 	std::istream& in = std::cin;
@@ -66,7 +67,7 @@
 		Point p(x,y,z);
 		Dt.insert(p);
 	}
-	std::cout << "Delaunay triangulation computed" << std::endl;
+	rInfo("Delaunay triangulation computed");
  
 	AlphaSimplex3DVector alpha_ordering;
 	fill_alpha_order(Dt, alpha_ordering);
@@ -83,17 +84,17 @@
 			continue;
 		
 		double current_alpha = CGAL::to_double(cur->value());
-		std::cout << "Current alpha: " << current_alpha << std::endl;
+		rInfo("Current alpha: %f", current_alpha);
 		std::sort(radius_ordering.begin(), radius_ordering.end(), ro);
-		std::cout << "Radius ordering size: " << radius_ordering.size() << std::endl;
+		rInfo("Radius ordering size: %i", radius_ordering.size());
 
 		RadiusFiltration rf;
 		for (SimplexVector::const_iterator cur = radius_ordering.begin(); cur != radius_ordering.end(); ++cur)
 			rf.append(*cur);
 		rf.fill_simplex_index_map();
-		std::cout << "Simplex index map filled" << std::endl;
+		rInfo("Simplex index map filled");
 		rf.pair_simplices(rf.begin(), rf.end());
-		std::cout << "Pairing computed" << std::endl;
+		rInfo("Pairing computed");
 	
 		for (RadiusFiltration::const_Index cur = rf.begin(); cur != rf.end(); ++cur)
 		{
@@ -101,7 +102,7 @@
 	
 			RealValue d1 = cur->distance;
 			//if (cur == cur->pair())
-			//	std::cout << "Unpaired " << cur->dimension() << ' ' << CGAL::to_double(d1) << std::endl;
+			//	rInfo("Unpaired %d %f", cur->dimension(), CGAL::to_double(d1));
 			
 			RealValue d2 = cur->pair()->distance;
 			if (d1 == d2)	continue;
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/alphashapes2d.cpp	Tue Feb 19 06:56:49 2008 -0500
@@ -0,0 +1,59 @@
+#include <utilities/log.h>
+
+#include "alphashapes2d.h"
+#include <topology/filtration.h>
+#include <iostream>
+#include <fstream>
+
+
+typedef Filtration<AlphaSimplex2D>				AlphaFiltration;
+
+int main(int argc, char** argv) 
+{
+#ifdef LOGGING
+	rlog::RLogInit(argc, argv);
+
+	stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
+	stdoutLog.subscribeTo( RLOG_CHANNEL("info") );
+	//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/filtration") );
+	//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/cycle") );
+#endif
+
+	SetFrequency(GetCounter("filtration/pair"), 10000);
+	SetTrigger(GetCounter("filtration/pair"), GetCounter(""));
+
+	// Read in the point set and compute its Delaunay triangulation
+	std::istream& in = std::cin;
+	double x,y;
+	Delaunay Dt;
+	while(in)
+	{
+		in >> x >> y;
+		Point p(x,y);
+		Dt.insert(p);
+	}
+	rInfo("Delaunay triangulation computed");
+   
+	AlphaSimplex2DVector alpha_ordering;
+	fill_alpha_order(Dt, alpha_ordering);
+	rInfo("Simplices: %i", alpha_ordering.size());
+
+	// Create the alpha-shape filtration
+	AlphaFiltration af;
+	for (AlphaSimplex2DVector::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");
+
+	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;
+
+}
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/alphashapes2d.h	Tue Feb 19 06:56:49 2008 -0500
@@ -0,0 +1,80 @@
+/**
+ * Author: Dmitriy Morozov
+ * Department of Computer Science, Duke University, 2007
+ */
+
+#ifndef __ALPHASHAPES2D_H__
+#define __ALPHASHAPES2D_H__
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Delaunay_triangulation_2.h>
+
+#include <topology/simplex.h>
+#include <utilities/types.h>
+
+#include <vector>
+#include <set>
+#include <iostream>
+
+struct K: CGAL::Exact_predicates_exact_constructions_kernel {};
+
+typedef CGAL::Delaunay_triangulation_2<K>    		Delaunay;
+typedef Delaunay::Point                				Point;
+typedef Delaunay::Vertex            				Vertex;
+typedef Delaunay::Vertex_handle            			Vertex_handle;
+typedef Delaunay::Edge								Edge;
+typedef Delaunay::Face								Face;
+typedef Delaunay::Face_handle						Face_handle;
+typedef K::FT										RealValue;
+
+typedef Delaunay::Finite_vertices_iterator    		Vertex_iterator;
+typedef Delaunay::Finite_edges_iterator        		Edge_iterator;
+typedef Delaunay::Finite_faces_iterator        		Face_iterator;
+
+
+class AlphaSimplex2D: public SimplexWithVertices<Vertex_handle>
+{
+	public:
+		typedef 	std::set<AlphaSimplex2D>							SimplexSet;
+		typedef		SimplexWithVertices<Vertex_handle>					Parent;
+		typedef		Parent::VertexContainer								VertexSet;
+		typedef		std::list<AlphaSimplex2D>							Cycle;
+
+    public:
+									AlphaSimplex2D()					{}
+									AlphaSimplex2D(const Parent& p): 
+											Parent(p) 					{}
+									AlphaSimplex2D(const AlphaSimplex2D& s): 
+											Parent(s) 					{ attached_ = s.attached_; alpha_ = s.alpha_; }
+	    							AlphaSimplex2D(const ::Vertex& v);
+		
+								    AlphaSimplex2D(const Edge& e);
+								    AlphaSimplex2D(const Edge& e, const SimplexSet& simplices);
+		
+									AlphaSimplex2D(const Face& c);
+	    
+		RealType					value() const						{ return CGAL::to_double(alpha_); }
+		RealValue					alpha() const						{ return alpha_; }
+		bool						attached() const					{ return attached_; }
+		Cycle						boundary() const;
+
+		// Ordering
+		struct AlphaOrder
+		{ bool operator()(const AlphaSimplex2D& first, const AlphaSimplex2D& second) const; };
+		
+		std::ostream& 				operator<<(std::ostream& out) const;
+		
+	private:
+		RealValue 					alpha_;
+		bool 						attached_;
+};
+
+typedef 			std::vector<AlphaSimplex2D>								AlphaSimplex2DVector;
+void 				fill_alpha_order(const Delaunay& Dt, 
+									 AlphaSimplex2DVector& alpha_order);
+
+std::ostream& 		operator<<(std::ostream& out, const AlphaSimplex2D& s)	{ return s.operator<<(out); }
+
+#include "alphashapes2d.hpp"
+
+#endif // __ALPHASHAPES2D_H__
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/alphashapes2d.hpp	Tue Feb 19 06:56:49 2008 -0500
@@ -0,0 +1,121 @@
+AlphaSimplex2D::	    
+AlphaSimplex2D(const ::Vertex& v): alpha_(0), attached_(false)
+{
+	for (int i = 0; i < 3; ++i)
+		if (v.face()->vertex(i)->point() == v.point())
+			Parent::add(v.face()->vertex(i));
+}
+
+AlphaSimplex2D::	    
+AlphaSimplex2D(const Edge& e): attached_(false)
+{
+    Face_handle f = e.first;
+	for (int i = 0; i < 3; ++i)
+		if (i != e.second)
+			Parent::add(f->vertex(i));
+}
+
+AlphaSimplex2D::	    
+AlphaSimplex2D(const Edge& e, const SimplexSet& simplices): attached_(false)
+{
+    Face_handle f = e.first;
+	for (int i = 0; i < 3; ++i)
+		if (i != e.second)
+			Parent::add(f->vertex(i));
+
+	Face_handle o = f->neighbor(e.second);
+	int oi = o->index(f);
+
+	VertexSet::const_iterator v = Parent::vertices().begin();
+	const Point& p1 = (*v++)->point();
+	const Point& p2 = (*v)->point();
+	
+	attached_ = false;
+	if (CGAL::side_of_bounded_circle(p1, p2, 
+									 f->vertex(e.second)->point()) == CGAL::ON_BOUNDED_SIDE)
+		attached_ = true;
+	else if (CGAL::side_of_bounded_circle(p1, p2,
+										  o->vertex(oi)->point()) == CGAL::ON_BOUNDED_SIDE)
+			attached_ = true;
+	else
+		alpha_ = squared_radius(p1, p2);
+
+	if (attached_)
+	{
+		SimplexSet::const_iterator f_iter = simplices.find(AlphaSimplex2D(*f));
+		SimplexSet::const_iterator o_iter = simplices.find(AlphaSimplex2D(*o));
+		if (f_iter == simplices.end())			// f is infinite
+			alpha_ = o_iter->alpha();
+		else if (o_iter == simplices.end())		// o is infinite
+			alpha_ = f_iter->alpha();
+		else
+			alpha_ = std::min(f_iter->alpha(), o_iter->alpha());
+	}
+}
+
+AlphaSimplex2D::	    
+AlphaSimplex2D(const Face& f): attached_(false)
+{
+	for (int i = 0; i < 3; ++i)
+		Parent::add(f.vertex(i));
+	VertexSet::const_iterator v = Parent::vertices().begin();
+	Point p1 = (*v++)->point();
+	Point p2 = (*v++)->point();
+	Point p3 = (*v)->point();
+	alpha_ = CGAL::squared_radius(p1, p2, p3);
+}
+
+AlphaSimplex2D::Cycle
+AlphaSimplex2D::boundary() const
+{
+	Cycle bdry;
+	Parent::Cycle pbdry = Parent::boundary();
+	for (Parent::Cycle::const_iterator cur = pbdry.begin(); cur != pbdry.end(); ++cur)
+		bdry.push_back(*cur);
+	return bdry;
+}
+
+
+bool 
+AlphaSimplex2D::AlphaOrder::
+operator()(const AlphaSimplex2D& first, const AlphaSimplex2D& second) const
+{
+	if (first.alpha() == second.alpha())
+		return (first.dimension() < second.dimension());
+	else
+		return (first.alpha() < second.alpha()); 
+}
+
+std::ostream& 
+AlphaSimplex2D::
+operator<<(std::ostream& out) const
+{
+	for (VertexSet::const_iterator cur = Parent::vertices().begin(); 
+								   cur != Parent::vertices().end(); ++cur)
+		out << **cur << ", ";
+	out << "value = " << value();
+
+	return out;
+}
+
+
+void fill_alpha_order(const Delaunay& Dt, AlphaSimplex2DVector& alpha_order)
+{
+	// Compute all simplices with their alpha values and attachment information
+	AlphaSimplex2D::SimplexSet simplices;
+	for(Face_iterator cur = Dt.finite_faces_begin(); cur != Dt.finite_faces_end(); ++cur)
+		simplices.insert(AlphaSimplex2D(*cur));
+	rInfo("Faces inserted");
+	for(Edge_iterator cur = Dt.finite_edges_begin(); cur != Dt.finite_edges_end(); ++cur)
+		simplices.insert(AlphaSimplex2D(*cur, simplices));
+	rInfo("Edges inserted");
+	for(Vertex_iterator cur = Dt.finite_vertices_begin(); cur != Dt.finite_vertices_end(); ++cur)
+		simplices.insert(AlphaSimplex2D(*cur));
+	rInfo("Vertices inserted");
+    
+	// Sort simplices by their alpha values
+	alpha_order.resize(simplices.size());
+	std::copy(simplices.begin(), simplices.end(), alpha_order.begin());
+	std::sort(alpha_order.begin(), alpha_order.end(), AlphaSimplex2D::AlphaOrder());
+}
+
--- a/examples/alphashapes/alphashapes3d.cpp	Tue Feb 19 06:51:23 2008 -0500
+++ b/examples/alphashapes/alphashapes3d.cpp	Tue Feb 19 06:56:49 2008 -0500
@@ -1,5 +1,4 @@
-#include <utilities/sys.h>
-#include <utilities/debug.h>
+#include <utilities/log.h>
 
 #include "alphashapes3d.h"
 #include <topology/filtration.h>
@@ -11,6 +10,18 @@
 
 int main(int argc, char** argv) 
 {
+#ifdef LOGGING
+	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") );
+#endif
+
+	SetFrequency(GetCounter("filtration/pair"), 10000);
+	SetTrigger(GetCounter("filtration/pair"), GetCounter(""));
+
 	// Read in the point set and compute its Delaunay triangulation
 	std::istream& in = std::cin;
 	double x,y,z;
@@ -21,18 +32,27 @@
 		Point p(x,y,z);
 		Dt.insert(p);
 	}
-	std::cout << "Delaunay triangulation computed" << std::endl;
+	rInfo("Delaunay triangulation computed");
    
 	AlphaSimplex3DVector alpha_ordering;
 	fill_alpha_order(Dt, alpha_ordering);
-	std::cout << "Simplices: " << alpha_ordering.size() << std::endl;
+	rInfo("Simplices: %d", alpha_ordering.size());
 
 	// 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();
-	af.pair_simplices(af.begin(), af.end());
-	std::cout << "Simplices paired" << std::endl;
+	rInfo("Filled simplex-index map");
+	af.pair_simplices(af.begin(), af.end(), false);
+	rInfo("Simplices paired");
+
+	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;
+
 }
 
--- a/examples/alphashapes/alphashapes3d.hpp	Tue Feb 19 06:51:23 2008 -0500
+++ b/examples/alphashapes/alphashapes3d.hpp	Tue Feb 19 06:56:49 2008 -0500
@@ -163,16 +163,16 @@
 	AlphaSimplex3D::SimplexSet simplices;
 	for(Cell_iterator cur = Dt.finite_cells_begin(); cur != Dt.finite_cells_end(); ++cur)
 		simplices.insert(AlphaSimplex3D(*cur));
-	std::cout << "Cells inserted" << std::endl;
+	rInfo("Cells inserted");
 	for(Facet_iterator cur = Dt.finite_facets_begin(); cur != Dt.finite_facets_end(); ++cur)
 		simplices.insert(AlphaSimplex3D(*cur, simplices));
-	std::cout << "Facets inserted" << std::endl;
+	rInfo("Facets inserted");
 	for(Edge_iterator cur = Dt.finite_edges_begin(); cur != Dt.finite_edges_end(); ++cur)
 		simplices.insert(AlphaSimplex3D(*cur, simplices, Dt.incident_facets(*cur)));
-	std::cout << "Edges inserted" << std::endl;
+	rInfo("Edges inserted");
 	for(Vertex_iterator cur = Dt.finite_vertices_begin(); cur != Dt.finite_vertices_end(); ++cur)
 		simplices.insert(AlphaSimplex3D(*cur));
-	std::cout << "Vertices inserted" << std::endl;
+	rInfo("Vertices inserted");
     
 	// Sort simplices by their alpha values
 	alpha_order.resize(simplices.size());
--- a/examples/ar-vineyard/CMakeLists.txt	Tue Feb 19 06:51:23 2008 -0500
+++ b/examples/ar-vineyard/CMakeLists.txt	Tue Feb 19 06:56:49 2008 -0500
@@ -1,6 +1,7 @@
 set							(targets						
 							 ar-vineyard)
 							 
+add_definitions				(${cgal_cxxflags})
 foreach 					(t ${targets})
 	add_executable			(${t} ${t}.cpp ${external_sources})
 	target_link_libraries	(${t} ${libraries} ${cgal_libraries} ${synaps_libraries})
--- a/examples/ar-vineyard/ar-vineyard.cpp	Tue Feb 19 06:51:23 2008 -0500
+++ b/examples/ar-vineyard/ar-vineyard.cpp	Tue Feb 19 06:56:49 2008 -0500
@@ -1,5 +1,4 @@
-#include <utilities/sys.h>
-#include <utilities/debug.h>
+#include <utilities/log.h>
 
 #include "ar-vineyard.h"
 
@@ -9,14 +8,15 @@
 
 int main(int argc, char** argv) 
 {
-#ifdef CWDEBUG
-	Debug(dc::filtration.off());
-	Debug(dc::cycle.off());
-	Debug(dc::vineyard.off());
-	Debug(dc::transpositions.off());
-	Debug(dc::lsfiltration.off());
+#ifdef LOGGING
+	rlog::RLogInit(argc, argv);
 
-	dionysus::debug::init();
+	stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
+	//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/filtration") );
+	//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/cycle") );
+	//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/vineyard") );
+	//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/filtration/transpositions") );
+	//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/lowerstar") );
 #endif
 
 	// Read command-line arguments
--- a/examples/ar-vineyard/ar-vineyard.h	Tue Feb 19 06:51:23 2008 -0500
+++ b/examples/ar-vineyard/ar-vineyard.h	Tue Feb 19 06:56:49 2008 -0500
@@ -6,9 +6,6 @@
 #ifndef __AR_VINEYARD_H__
 #define __AR_VINEYARD_H__
 
-#include "utilities/sys.h"
-#include "utilities/debug.h"
-
 #include "topology/conesimplex.h"
 #include "topology/filtration.h"
 #include "geometry/kinetic-sort.h"
--- a/examples/ar-vineyard/ar-vineyard.hpp	Tue Feb 19 06:51:23 2008 -0500
+++ b/examples/ar-vineyard/ar-vineyard.hpp	Tue Feb 19 06:56:49 2008 -0500
@@ -1,3 +1,5 @@
+#include <utilities/log.h>
+
 /* Implementation */
 	
 void
--- a/examples/cech-complex/CMakeLists.txt	Tue Feb 19 06:51:23 2008 -0500
+++ b/examples/cech-complex/CMakeLists.txt	Tue Feb 19 06:56:49 2008 -0500
@@ -2,6 +2,6 @@
 							 cech-complex)
 							 
 foreach 					(t ${targets})
-	add_executable			(${t} ${t}.cpp ${external_sources})
-	target_link_libraries	(${t} ${libraries} ${cgal_libraries})
+	add_executable			(${t} ${t}.cpp)
+	target_link_libraries	(${t} ${libraries})
 endforeach 					(t ${targets})
--- a/examples/cech-complex/cech-complex.cpp	Tue Feb 19 06:51:23 2008 -0500
+++ b/examples/cech-complex/cech-complex.cpp	Tue Feb 19 06:56:49 2008 -0500
@@ -1,5 +1,4 @@
-#include <utilities/sys.h>
-#include <utilities/debug.h>
+#include <utilities/log.h>
 
 #include <topology/simplex.h>
 #include <topology/filtration.h>
@@ -74,6 +73,9 @@
 
 int main(int argc, char** argv) 
 {
+	SetFrequency(GetCounter("filtration/pair"), 10000);
+	SetTrigger(GetCounter("filtration/pair"), GetCounter(""));
+
 	// Read in the point set and compute its Delaunay triangulation
 	std::istream& in = std::cin;
 	int ambient_d, homology_d;	
@@ -115,7 +117,7 @@
 
 	// Compute persistence
 	cf.fill_simplex_index_map();
-	cf.pair_simplices(cf.begin(), cf.end());
+	cf.pair_simplices(cf.begin(), cf.end(), false);
 	std::cout << "Simplices paired" << std::endl;
 
 	for (CechFiltration::Index i = cf.begin(); i != cf.end(); ++i)
--- a/examples/grid/CMakeLists.txt	Tue Feb 19 06:51:23 2008 -0500
+++ b/examples/grid/CMakeLists.txt	Tue Feb 19 06:56:49 2008 -0500
@@ -3,7 +3,8 @@
 							 pdbdistance-vineyard
 							 combustion-vineyard)
 							 
+add_definitions				(${cgal_cxxflags})
 foreach 					(t ${targets})
-	add_executable			(${t} ${t}.cpp ${external_sources})
-	target_link_libraries	(${t} ${libraries} ${cgal_libraries})
+	add_executable			(${t} ${t}.cpp)
+	target_link_libraries	(${t} ${libraries} ${cgal_libraries} ${dsrpdb_LIBRARY})
 endforeach 					(t ${targets})
--- a/examples/grid/combustion-vineyard.cpp	Tue Feb 19 06:51:23 2008 -0500
+++ b/examples/grid/combustion-vineyard.cpp	Tue Feb 19 06:56:49 2008 -0500
@@ -3,8 +3,7 @@
  * Department of Computer Science, Duke University, 2005
  */
 
-#include <utilities/sys.h>
-#include <utilities/debug.h>
+#include <utilities/log.h>
 
 #include <iostream>
 #include <fstream>
--- a/examples/grid/grid2Dvineyard.h	Tue Feb 19 06:51:23 2008 -0500
+++ b/examples/grid/grid2Dvineyard.h	Tue Feb 19 06:56:49 2008 -0500
@@ -6,9 +6,6 @@
 #ifndef __GRID2DVINEYARD_H__
 #define __GRID2DVINEYARD_H__
 
-#include "utilities/sys.h"
-#include "utilities/debug.h"
-
 #include "grid2D.h"
 #include "topology/lowerstarfiltration.h"
 
--- a/examples/grid/grid2Dvineyard.hpp	Tue Feb 19 06:51:23 2008 -0500
+++ b/examples/grid/grid2Dvineyard.hpp	Tue Feb 19 06:56:49 2008 -0500
@@ -1,3 +1,5 @@
+#include <utilities/log.h>
+
 /* Implementation */
 	
 Grid2DVineyard::
@@ -108,25 +110,25 @@
 			VertexIndex vv(&vertices_[grid()->seq(x,y+1)]);
 			VertexIndex vd(&vertices_[grid()->seq(x+1,y+1)]);
 
-			Simplex sh(v);
+			Simplex sh(2, v);
 			sh.add(vh);	filtration_->append(sh);		// Horizontal edge
 			sh.add(vd);	filtration_->append(sh);		// "Horizontal" triangle
 			
-			Simplex sv(v);
+			Simplex sv(2, v);
 			sv.add(vv);	filtration_->append(sv);		// Vertical edge
 			sv.add(vd);	filtration_->append(sv);		// "Vertical" triangle
 			
-			Simplex sd(v);
+			Simplex sd(2, v);
 			sd.add(vd); filtration_->append(sd);		// Diagonal edge
 
 			if (y == grid()->ysize() - 2)
 			{
-				Simplex s(vv); 
+				Simplex s(1, vv); 
 				s.add(vd); filtration_->append(s);		// Top edge
 			}
 			if (x == grid()->xsize() - 2)
 			{
-				Simplex s(vh); 
+				Simplex s(1, vh); 
 				s.add(vd); filtration_->append(s);		// Right edge
 			}
 		}
--- a/examples/grid/pdbdistance-vineyard.cpp	Tue Feb 19 06:51:23 2008 -0500
+++ b/examples/grid/pdbdistance-vineyard.cpp	Tue Feb 19 06:56:49 2008 -0500
@@ -1,6 +1,5 @@
 //#include <boost/archive/binary_oarchive.hpp>
-#include "utilities/sys.h"
-#include "utilities/debug.h"
+#include "utilities/log.h"
 
 #include "pdbdistance.h"
 #include "grid2Dvineyard.h"
@@ -18,14 +17,13 @@
 
 int main(int argc, char** argv)
 {
-#ifdef CWDEBUG
-	Debug(dc::filtration.off());
-	Debug(dc::cycle.off());
-	Debug(dc::vineyard.off());
-	Debug(dc::transpositions.off());
-	Debug(dc::lsfiltration.off());
+#ifdef LOGGING
+	rlog::RLogInit(argc, argv);
 
-	dionysus::debug::init();
+	stdoutLog.subscribeTo( RLOG_CHANNEL("topology/filtration") );
+	stdoutLog.subscribeTo( RLOG_CHANNEL("topology/cycle") );
+	stdoutLog.subscribeTo( RLOG_CHANNEL("topology/vineyard") );
+	stdoutLog.subscribeTo( RLOG_CHANNEL("topology/lowerstar") );
 #endif
 
 	if (argc < 5)
--- a/examples/triangle/CMakeLists.txt	Tue Feb 19 06:51:23 2008 -0500
+++ b/examples/triangle/CMakeLists.txt	Tue Feb 19 06:56:49 2008 -0500
@@ -2,6 +2,6 @@
 							 triangle)
 							 
 foreach 					(t ${targets})
-	add_executable			(${t} ${t}.cpp ${external_sources})
+	add_executable			(${t} ${t}.cpp)
 	target_link_libraries	(${t} ${libraries})
 endforeach 					(t ${targets})
--- a/examples/triangle/triangle.cpp	Tue Feb 19 06:51:23 2008 -0500
+++ b/examples/triangle/triangle.cpp	Tue Feb 19 06:56:49 2008 -0500
@@ -1,3 +1,5 @@
+#include <utilities/log.h>
+
 #include "topology/filtration.h"
 #include "topology/simplex.h"
 #include <vector>
@@ -24,15 +26,14 @@
 	f.append(Simplex(bg,     bg + 3, 5));				// ABC
 }
 
-int main()
+int main(int argc, char** argv)
 {
-#ifdef CWDEBUG
-	dionysus::debug::init();
+#ifdef LOGGING
+	rlog::RLogInit(argc, argv);
 
-	Debug(dc::filtration.on());
-	Debug(dc::cycle.off());
-	Debug(dc::vineyard.on());
-	Debug(dc::transpositions.on());
+	stdoutLog.subscribeTo(RLOG_CHANNEL("topology/filtration"));
+	//stdoutLog.subscribeTo(RLOG_CHANNEL("topology/cycle"));
+	stdoutLog.subscribeTo(RLOG_CHANNEL("topology/vineyard"));
 #endif
 
 	Evaluator<Simplex> e;
--- a/include/topology/cycle.h	Tue Feb 19 06:51:23 2008 -0500
+++ b/include/topology/cycle.h	Tue Feb 19 06:56:49 2008 -0500
@@ -6,9 +6,6 @@
 #ifndef __CYCLE_H__
 #define __CYCLE_H__
 
-#include "utilities/sys.h"
-#include "utilities/debug.h"
-
 #include "utilities/types.h"
 #include "utilities/circular_list.h"
 #include <list>
--- a/include/topology/cycle.hpp	Tue Feb 19 06:51:23 2008 -0500
+++ b/include/topology/cycle.hpp	Tue Feb 19 06:56:49 2008 -0500
@@ -6,17 +6,29 @@
 #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(): sz(0)
+Cycle()
 {}
 
 template<class I, class OrderCmp, class ConsistencyCmp>
 Cycle<I,OrderCmp,ConsistencyCmp>::
-Cycle(const Cycle& c): CycleRepresentation(c), sz(c.sz)					
+Cycle(const Cycle& c): CycleRepresentation(c)
 {}
 
 template<class I, class OrderCmp, class ConsistencyCmp>
@@ -59,7 +71,6 @@
 swap(Cycle& c)
 {
 	CycleRepresentation::swap(c);
-	std::swap(sz, c.sz);
 }
 
 template<class I, class OrderCmp, class ConsistencyCmp>
@@ -80,7 +91,7 @@
 	AssertMsg(!empty(), "Cycle must not be empty for get_second()");
 	if (size() < 2)			return begin();					// Indicates that there is no second.
 
-	Dout(dc::cycle, "Looking for 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();
@@ -98,7 +109,7 @@
 		}
 	}
 	
-	Dout(dc::cycle, "Looked up: " << **second);
+	rLog(rlCycle, "Looked up: %s", tostring(**second).c_str());
 	return second;
 }
 
@@ -143,7 +154,7 @@
 	for (iterator cur = begin(); cur != end(); ++cur)
 		while ((cur != end()) && ((*cur == j) || (cmp(*cur, j) && cmp(i, *cur))))
 		{
-			Dout(dc::cycle, "Iteration of the erase while loop");
+			rLog(rlCycle, "Iteration of the erase while loop");
 			cur = erase(cur);
 		}
 }
@@ -173,7 +184,7 @@
 Cycle<I, OrderCmp, ConsistencyCmp>::
 add(const Self& c, const ConsistencyCmp& cmp)
 {
-	Dout(dc::cycle, "Adding cycles: " << *this << " + " << c);
+	rLog(rlCycle, "Adding cycles: %s + %s",  tostring(*this).c_str(), tostring(c).c_str());
 	
 	iterator 			cur1 = begin();
 	const_iterator 		cur2 = c.begin();
@@ -183,34 +194,37 @@
 		if (cur1 == end())
 		{
 			while (cur2 != c.end())
+			{
 				push_back(*cur2++);
-			Dout(dc::cycle, "After addition: " << *this);
+				Count(cCycleAddBasic);
+			}
+			rLog(rlCycle, "After addition: %s", tostring(*this).c_str());
 			return *this;
 		}
 
 		// mod 2
 		int res = cmp.compare(*cur1, *cur2);
-		Dout(dc::cycle, "Comparison result: " << res);
+		Count(cCycleAddComparison);
+		rLog(rlCycle, "Comparison result: %i", res);
 		if (res == 0)		// *cur1 == *cur2
 		{
-			Dout(dc::cycle, "Equality");
+			rLog(rlCycle, "Equality");
 			cur1 = erase(cur1);		// erase cur1 --- as a result cur1 will be pointing at old_cur1++
-			--sz;
 			++cur2;
 		} else if (res < 0)	// *cur1 < *cur2
 		{
-			Dout(dc::cycle, "Less than");
+			rLog(rlCycle, "Less than");
 			cur1++;
 		} else if (res > 0) // *cur1 > *cur2
 		{
-			Dout(dc::cycle, "Greater than");
+			rLog(rlCycle, "Greater than");
 			insert(cur1, *cur2);
 			++cur2;
-			++sz;
 		}
+		Count(cCycleAddBasic);
 	}
 
-	Dout(dc::cycle, "After addition: " << *this);
+	rLog(rlCycle, "After addition: %s", tostring(*this).c_str());
 	return *this;
 }
 
@@ -238,7 +252,6 @@
 serialize(Archive& ar, version_type )
 {
 	ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP(Parent);
-	ar & make_nvp("size", sz);;
 }
 
 
--- a/include/topology/filtration.h	Tue Feb 19 06:51:23 2008 -0500
+++ b/include/topology/filtration.h	Tue Feb 19 06:56:49 2008 -0500
@@ -6,8 +6,7 @@
 #ifndef __FILTRATION_H__
 #define __FILTRATION_H__
 
-#include "utilities/sys.h"
-#include "utilities/debug.h"
+#include "utilities/log.h"
 
 #include "filtrationcontainer.h"
 #include "filtrationsimplex.h"
@@ -62,8 +61,8 @@
 		/// \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);
-		void 							pair_simplices()							{ pair_simplices(begin(), end()); }
+		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 = true);
 		bool							is_paired() const;
 		Index							append(const Simplex& s);					///< Appends s to the filtration
@@ -94,8 +93,8 @@
 		
 		std::ostream& 					operator<<(std::ostream& out) const;
 
-	protected:
-		/// \name Comparator accessors (protected)
+	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; }
@@ -110,6 +109,7 @@
 		void							pairing_switch(Index i, Index j);
 		
 		bool 							paired;
+		bool							trails_stored;
 		SimplexMap						inverse_simplices;
 
 		Vineyard*						vineyard_;
--- a/include/topology/filtration.hpp	Tue Feb 19 06:51:23 2008 -0500
+++ b/include/topology/filtration.hpp	Tue Feb 19 06:56:49 2008 -0500
@@ -1,5 +1,5 @@
+#include "utilities/types.h"
 #include "utilities/counter.h"
-#include "utilities/types.h"
 #include <algorithm>
 
 #include <boost/utility.hpp>
@@ -11,6 +11,30 @@
 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*  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>::
@@ -20,44 +44,54 @@
 template<class S, class FS, class V>
 void 
 Filtration<S, FS, V>::
-pair_simplices(Index bg, Index end)
+pair_simplices(Index bg, Index end, bool store_trails)
 {
-	Dout(dc::filtration, "Entered: compute_pairing");
+	if (!is_paired())
+		trails_stored = store_trails;
+	else
+		trails_stored &= store_trails;
+
+	rLog(rlFiltration, "Entered: compute_pairing");
 	for (Index j = bg; j != end; ++j)
 	{
-		Dout(dc::filtration|flush_cf|continued_cf, *j << ": ");
+		rLog(rlFiltration, "Simplex %s", tostring(*j).c_str());
 		init_cycle_trail(j); 
 		Cycle& bdry = j->cycle();
-		Dout(dc::finish, bdry);
+		rLog(rlFiltration, "  has boundary: %s", tostring(bdry).c_str());
 		
-		CountNum("Boundaries", j->dimension());
-		Count("SimplexCount");
+		CountNum(cFiltrationPairBoundaries, j->dimension());
+		Count(cFiltrationPair);
 
 		while(!bdry.empty())
 		{
 			Index i = bdry.top(cycles_cmp);
-			Dout(dc::filtration, *i << ": " << *(i->pair()));
-			AssertMsg(!cycles_cmp(i, j), "Simplices in the cycle must precede current simplex: " << 
-										 "(" << *i << " in cycle of " << *j << ")");
+			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)
 			{
-				Dout(dc::filtration, "Pairing " << *i << " and " << *j << " with cycle " << j->cycle());
+				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("DepositedCycleLength", j->cycle().size());
+				CountNum(cFiltrationPairCycleLength, j->cycle().size());
+				CountBy(cFiltrationPairCycleLength, j->cycle().size());
 				break;
 			}
 
-			// continue searching --- change the Dout to the continued mode with newlines FIXME
-			Dout(dc::filtration, "  Adding: [" << bdry << "] + ");
-			Dout(dc::filtration, "          [" << i->pair()->cycle() << "]");
+			rLog(rlFiltration, "  Adding: [%s] + [%s]", 
+							   tostring(bdry).c_str(), tostring(i->pair()->cycle()).c_str());
 			bdry.add(i->pair()->cycle(), get_consistency_cmp());
-			i->pair()->trail().append(j, get_consistency_cmp());
-			Dout(dc::filtration, "After addition: " << bdry);
+			if (store_trails)	i->pair()->trail().append(j, get_consistency_cmp());
+			Count(cFiltrationPairTrailLength);
+			rLog(rlFiltration, "After addition: %s", tostring(bdry).c_str());
 		}
-		Dout(dc::filtration, "Finished with " << *j << ": " << *(j->pair()));
+		rLog(rlFiltration, "Finished with %s: %s", 
+						   tostring(*j).c_str(), tostring(*(j->pair())).c_str());
 	}
 	paired = true;
 }
@@ -78,6 +112,7 @@
 transpose(Index i, bool maintain_lazy)
 {
 	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++;
 	
@@ -167,28 +202,28 @@
 transpose_simplices(Index i, bool maintain_lazy)
 {
 	AssertMsg(is_paired(), "Pairing must be computed before transpositions");
-	Count("SimplexTransposition");
+	Count(cFiltrationTransposition);
 	
 	Index i_prev = i++;
 
 	if (i_prev->dimension() != i->dimension())
 	{
 		swap(i_prev, i);
-		Dout(dc::transpositions, "Different dimension");
-		Count("Case DiffDim");
+		rLog(rlFiltrationTranspositions, "Different dimension");
+		Count(cFiltrationTranspositionDiffDim);
 		return false;
 	}
 	
 	bool si = i_prev->sign(), sii = i->sign();
 	if (si && sii)
 	{
-		Dout(dc::transpositions, "Trail prev: " << i_prev->trail());
+		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())
 		{
-			Dout(dc::transpositions, "Case 1, U[i,i+1] = 1");
+			rLog(rlFiltrationTranspositions, "Case 1, U[i,i+1] = 1");
 			i_prev->trail().erase(i_in_i_prev);
 		}
 
@@ -199,9 +234,9 @@
 		if (l == i)
 		{
 			swap(i_prev, i);
-			Dout(dc::transpositions, "Case 1.2 --- unpaired");
-			Dout(dc::transpositions, *i_prev);
-			Count("Case 1.2");
+			rLog(rlFiltrationTranspositions, "Case 1.2 --- unpaired");
+			rLog(rlFiltrationTranspositions, "%s", tostring(*i_prev).c_str());
+			Count(cFiltrationTranspositionCase12);
 			return false;
 		} else if (k == i_prev)
 		{
@@ -209,23 +244,23 @@
 			{
 				// Case 1.2
 				swap(i_prev, i);
-				Dout(dc::transpositions, "Case 1.2 --- unpaired");
-				Dout(dc::transpositions, *i_prev);
-				Count("Case 1.2");
+				rLog(rlFiltrationTranspositions, "Case 1.2 --- unpaired");
+				rLog(rlFiltrationTranspositions, tostring(*i_prev).c_str());
+				Count(cFiltrationTranspositionCase12);
 				return false;
 			} else
 			{
 				// Case 1.1.2 --- special version (plain swap, but pairing switches)
 				swap(i_prev, i);
 				pairing_switch(i_prev, i);
-				Dout(dc::transpositions, "Case 1.1.2 --- unpaired");
-				Dout(dc::transpositions, *i_prev);
-				Count("Case 1.1.2");
+				rLog(rlFiltrationTranspositions, "Case 1.1.2 --- unpaired");
+				rLog(rlFiltrationTranspositions, tostring(*i_prev).c_str());
+				Count(cFiltrationTranspositionCase112);
 				return true;
 			}
 		}
 		
-		Dout(dc::transpositions, "l cycle: " << l->cycle());
+		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
@@ -239,8 +274,8 @@
 				}
 			}
 			swap(i_prev, i);
-			Dout(dc::transpositions, "Case 1.2");
-			Count("Case 1.2");
+			rLog(rlFiltrationTranspositions, "Case 1.2");
+			Count(cFiltrationTranspositionCase12);
 			return false;
 		} else
 		{
@@ -251,8 +286,8 @@
 				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
-				Dout(dc::transpositions, "Case 1.1.1");
-				Count("Case 1.1.1");
+				rLog(rlFiltrationTranspositions, "Case 1.1.1");
+				Count(cFiltrationTranspositionCase111);
 				return false;
 			} else
 			{
@@ -261,8 +296,8 @@
 				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);
-				Dout(dc::transpositions, "Case 1.1.2");
-				Count("Case 1.1.2");
+				rLog(rlFiltrationTranspositions, "Case 1.1.2");
+				Count(cFiltrationTranspositionCase112);
 				return true;
 			}
 		}
@@ -273,8 +308,8 @@
 		{
 			// Case 2.2
 			swap(i_prev, i);
-			Dout(dc::transpositions, "Case 2.2");
-			Count("Case 2.2");
+			rLog(rlFiltrationTranspositions, "Case 2.2");
+			Count(cFiltrationTranspositionCase22);
 			return false;
 		} else
 		{
@@ -290,14 +325,14 @@
 				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);
-				Dout(dc::transpositions, "Case 2.1.2");
-				Count("Case 2.1.2");
+				rLog(rlFiltrationTranspositions, "Case 2.1.2");
+				Count(cFiltrationTranspositionCase212);
 				return true;
 			} 
 			
 			// Case 2.1.1
-			Dout(dc::transpositions, "Case 2.1.1");
-			Count("Case 2.1.1");
+			rLog(rlFiltrationTranspositions, "Case 2.1.1");
+			Count(cFiltrationTranspositionCase211);
 			return false;
 		}
 	} else if (!si && sii)
@@ -307,8 +342,8 @@
 		{
 			// Case 3.2
 			swap(i_prev, i);
-			Dout(dc::transpositions, "Case 3.2");
-			Count("Case 3.2");
+			rLog(rlFiltrationTranspositions, "Case 3.2");
+			Count(cFiltrationTranspositionCase32);
 			return false;
 		} else
 		{
@@ -319,8 +354,8 @@
 			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);
-			Dout(dc::transpositions, "Case 3.1");
-			Count("Case 3.1");
+			rLog(rlFiltrationTranspositions, "Case 3.1");
+			Count(cFiltrationTranspositionCase31);
 			return true;
 		}
 	} else if (si && !sii)
@@ -329,12 +364,12 @@
 		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())
 		{
-			Dout(dc::transpositions, "Case 4, U[i,i+1] = 1");
+			rLog(rlFiltrationTranspositions, "Case 4, U[i,i+1] = 1");
 			i_prev->trail().erase(i_in_i_prev);
 		}
 		swap(i_prev, i);
-		Dout(dc::transpositions, "Case 4");
-		Count("Case 4");
+		rLog(rlFiltrationTranspositions, "Case 4");
+		Count(cFiltrationTranspositionCase4);
 		return false;
 	}
 	
@@ -352,7 +387,8 @@
 
 	for (typename Simplex::Cycle::const_iterator cur = bdry.begin(); cur != bdry.end(); ++cur)
 	{
-		Dout(dc::filtration, "Appending in init_cycle_trail(): " << *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());
 	}
@@ -400,7 +436,7 @@
 
 	SizeType sz = size();
 	ar << make_nvp("size", sz);
-	Dout(dc::filtration, "Size: " << sz);
+	rLog(rlFiltration, "Size: %i", sz);
 
 	/* Record integer indices */
 	IndexIntMap index_map; SizeType i = 0;
@@ -416,7 +452,7 @@
 		//FiltrationSimplexSerialization simplex = FiltrationSimplexSerialization(*cur, index_map);
 		//ar << make_nvp("FiltrationSimplex", simplex);	
 	}
-	Dout(dc::filtration, count << " simplices serialized");
+	rLog(rlFiltration, "%i simplices serialized", count);
 }
 
 template<class S, class FS, class V>
@@ -425,16 +461,16 @@
 Filtration<S, FS, V>::
 load(Archive& ar, version_type )
 {
-	Dout(dc::filtration, "Starting to read filtration");
+	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);
-	Dout(dc::filtration, "In Filtration: first block read");
+	rLog(rlFiltration, "In Filtration: first block read");
 
 	SizeType sz;
 	ar >> make_nvp("size", sz);
-	Dout(dc::filtration, "In Filtration: size read " << sz);
+	rLog(rlFiltration, "In Filtration: size read %i", sz);
 	
 	IndexVector index_vector(sz);
 	for (SizeType i = 0; i < sz; ++i)
@@ -449,10 +485,10 @@
 		//FiltrationSimplexSerialization simplex;
 		//ar >> make_nvp("FiltrationSimplex", simplex);
 		count++;
-		Dout(dc::filtration, "In Filtration: simplex read (" << count << ")");
+		rLog(rlFiltration, "In Filtration: simplex read (%i)", count);
 		//simplex.set_filtration_simplex(*index_vector[i], index_vector);
 	}
-	Dout(dc::filtration, "In Filtration: simplices read");
+	rLog(rlFiltration, "In Filtration: simplices read");
 }
 
 template<class S, class FS, class V>
--- a/include/topology/filtrationsimplex.h	Tue Feb 19 06:51:23 2008 -0500
+++ b/include/topology/filtrationsimplex.h	Tue Feb 19 06:56:49 2008 -0500
@@ -6,9 +6,6 @@
 #ifndef __FILTRATIONSIMPLEX_H__
 #define __FILTRATIONSIMPLEX_H__
 
-#include "utilities/sys.h"
-#include "utilities/debug.h"
-
 #include "filtrationcontainer.h"
 #include "vineyard.h"
 #include "utilities/types.h"
--- a/include/topology/lowerstarfiltration.h	Tue Feb 19 06:51:23 2008 -0500
+++ b/include/topology/lowerstarfiltration.h	Tue Feb 19 06:56:49 2008 -0500
@@ -6,8 +6,7 @@
 #ifndef __LOWERSTARFILTRATION_H__
 #define __LOWERSTARFILTRATION_H__
 
-#include "utilities/sys.h"
-#include "utilities/debug.h"
+#include "utilities/log.h"
 
 #include "filtration.h"
 #include "simplex.h"
--- a/include/topology/lowerstarfiltration.hpp	Tue Feb 19 06:51:23 2008 -0500
+++ b/include/topology/lowerstarfiltration.hpp	Tue Feb 19 06:56:49 2008 -0500
@@ -1,5 +1,16 @@
+#include "utilities/counter.h"
+
 /* Implementations */
 
+#ifdef LOGGING
+static rlog::RLogChannel* rlLowerStar = 				DEF_CHANNEL("lowerstar", rlog::Log_Debug);
+#endif // LOGGING
+
+#ifdef COUNTERS
+static Counter*  cLowerStarTransposition =		 		GetCounter("lowerstar");
+static Counter*  cLowerStarChangedAttachment =		 	GetCounter("lowerstar/ChangedAttachment");
+#endif // COUNTERS
+
 template<class VI, class Smplx, class FltrSmplx, class Vnrd>
 template<class VertexCmp>
 LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::
@@ -17,7 +28,7 @@
 
 	// Record vertex order
 	for(typename VertexIndexList::const_iterator cur = tmp_list.begin(); cur != tmp_list.end(); ++cur)
-		(*cur)->set_order(vertex_order.push_back(VertexDescriptor(*cur, Parent::append(Simplex(*cur)))));
+		(*cur)->set_order(vertex_order.push_back(VertexDescriptor(*cur, Parent::append(Simplex(0, *cur)))));
 }
 
 template<class VI, class Smplx, class FltrSmplx, class Vnrd>
@@ -59,22 +70,10 @@
 LowerStarFiltration<VI,Smplx,FltrSmplx,Vnrd>::
 transpose_vertices(const VertexOrderIndex& order)
 {
-	Count("VertexTransposition");
-
-#if COUNTERS
-	if ((counters.lookup("VertexTransposition") % 1000000) == 0)
-	{
-		Dout(dc::lsfiltration, "Vertex transpositions:  " << counters.lookup("VertexTransposition"));
-		Dout(dc::lsfiltration, "Simplex transpositions: " << counters.lookup("SimplexTransposition"));
-		Dout(dc::lsfiltration, "Attachment changed:     " << counters.lookup("ChangedAttachment"));
-		Dout(dc::lsfiltration, "Regular disconnected:   " << counters.lookup("RegularDisconnected"));
-		Dout(dc::lsfiltration, "Pairing Changed:        " << counters.lookup("ChangedPairing"));
-		Dout(dc::lsfiltration, "------------------------");
-	}
-#endif // COUNTERS
-	
-	Dout(dc::lsfiltration, "Transposing vertices (" << order->vertex_index << ", " 
-													<< boost::next(order)->vertex_index << ")");
+	Count(cLowerStarTransposition);
+	rLog(rlLowerStar, "Transposing vertices (%s, %s)", 
+					  tostring(order->vertex_index).c_str(),
+					  tostring(boost::next(order)->vertex_index).c_str());
 
 	Index i = order->simplex_index;
 	Index i_prev = boost::prior(i);
@@ -94,17 +93,17 @@
 		result |= transpose(i_next_prev);
 		i_next_prev = boost::prior(i_next);
 	}
-	Dout(dc::lsfiltration, "Done moving the vertex");
+	rLog(rlLowerStar, "Done moving the vertex");
 
 	// Second, move the simplices attached to it
-	Dout(dc::lsfiltration, "Moving attached simplices");
+	rLog(rlLowerStar, "Moving attached simplices");
 	while (j != Parent::end() && j->get_attachment() == v_i_next)
 	{
-		Dout(dc::lsfiltration, "  Considering " << *j);
+		rLog(rlLowerStar, "  Considering %s", tostring(*j).c_str());
 		if (nbghrs && j->contains(v_i))			// short circuit
 		{
-			Count("ChangedAttachment");
-			Dout(dc::lsfiltration, "  Attachment changed for " << *j);
+			Count(cLowerStarChangedAttachment);
+			rLog(rlLowerStar, "  Attachment changed for %s", tostring(*j).c_str());
 			j->set_attachment(v_i);
 			++j;
 			continue;
@@ -114,13 +113,14 @@
 		while ((--j_prev)->get_attachment() != v_i_next) 		// i.e., until we have reached v_i_next 
 															// (and the simplices that follow it) again
 		{
-			Dout(dc::lsfiltration, "    Moving: " << *j_prev << ", " << *boost::next(j_prev));
+			rLog(rlLowerStar, "    Moving: %s, %s", 
+							  tostring(*j_prev).c_str(), tostring(*boost::next(j_prev)).c_str());
 			AssertMsg(j_prev->get_attachment() == v_i, "Simplex preceding the one being moved must be attached to v_i");
 			result |= transpose(j_prev);
 			--j_prev;
 		}
 	}
-	Dout(dc::lsfiltration, "Done moving attached simplices");
+	rLog(rlLowerStar, "Done moving attached simplices");
 	vertex_order.swap(order, boost::next(order));
 	
 	return result;
@@ -133,14 +133,17 @@
 {
 	Index j = boost::next(i);
 	
-	Dout(dc::lsfiltration, "    Transposing (" << *i << ", " << *(i->pair()) << ") and (" 
-											   << *j << ", " << *(j->pair()) << ")");
+	rLog(rlLowerStar, "    Transposing (%s, %s) and (%s, %s)", 
+					  tostring(*i).c_str(), tostring(*(i->pair())).c_str(),
+					  tostring(*j).c_str(), tostring(*(j->pair())).c_str());
 
 	assert_pairing(i);
 	assert_pairing(j);
 
 	bool res = Parent::transpose(i, false);
-	Dout(dc::lsfiltration, "    " << *j << ": " << *(j->pair()) << ", " << *i << ": " << *(i->pair()));
+	rLog(rlLowerStar, "    %s: %s, %s: %s",
+					  tostring(*j).c_str(), tostring(*(j->pair())).c_str(),
+					  tostring(*i).c_str(), tostring(*(i->pair())).c_str());
 
 	assert_pairing(j);
 	assert_pairing(i);
@@ -159,10 +162,11 @@
 	{
 		if (i->pair() != i->cycle().top(Parent::get_cycles_cmp()))
 		{
-			Dout(dc::lsfiltration, "i (negative): " << *i);
-			Dout(dc::lsfiltration, "pair(i): " << *(i->pair()));
-			Dout(dc::lsfiltration, "i->cycle().top(): " << *(i->cycle().top(Parent::get_cycles_cmp())));
-			DoutFatal(dc::fatal, "Pairing not matching the matrix at " << *i);
+			rLog(rlLowerStar, "i (negative): %s", tostring(*i).c_str());
+			rLog(rlLowerStar, "pair(i): %s", tostring(*(i->pair())).c_str());
+			rLog(rlLowerStar, "i->cycle().top(): %s", 
+							  tostring(*(i->cycle().top(Parent::get_cycles_cmp()))).c_str());
+			AssertMsg(0, "Pairing not matching the matrix at %s", tostring(*i).c_str());
 		}
 	} else
 	{
@@ -170,11 +174,11 @@
 		{
 			if (i->pair()->cycle().top(Parent::get_cycles_cmp()) != i)
 			{
-				Dout(dc::lsfiltration, "i (positive): " << *i);
-				Dout(dc::lsfiltration, "pair(i): " << *(i->pair()));
-				Dout(dc::lsfiltration, "pair(i)->cycle(): " << i->pair()->cycle());
-				Dout(dc::lsfiltration, "pair->cycle().top(): " << *(i->pair()->cycle().top(Parent::get_cycles_cmp())));
-				DoutFatal(dc::fatal, "Pairing not matching the matrix at " << *(i->pair()));
+				rLog(rlLowerStar, "i (positive): %s", tostring(*i).c_str());
+				rLog(rlLowerStar, "pair(i): %s", tostring(*(i->pair())).c_str());
+				rLog(rlLowerStar, "pair(i)->cycle(): %s", tostring(i->pair()->cycle()).c_str());
+				rLog(rlLowerStar, "pair->cycle().top(): %s", tostring(*(i->pair()->cycle().top(Parent::get_cycles_cmp()))).c_str());
+				AssertMsg(0, "Pairing not matching the matrix at %s", tostring(*(i->pair())).c_str());
 			}
 		}
 	}
--- a/include/topology/simplex.h	Tue Feb 19 06:51:23 2008 -0500
+++ b/include/topology/simplex.h	Tue Feb 19 06:56:49 2008 -0500
@@ -1,15 +1,13 @@
 /*
  * Author: Dmitriy Morozov
- * Department of Computer Science, Duke University, 2005 -- 2006
+ * Department of Computer Science, Duke University, 2005 -- 2007
  */
 
 #ifndef __SIMPLEX_H__
 #define __SIMPLEX_H__
 
-#include "utilities/sys.h"
-#include "utilities/debug.h"
-
-#include <set>
+#include <vector>
+#include <algorithm>
 #include <list>
 #include <iostream>
 
@@ -30,7 +28,7 @@
 		typedef		V																Vertex;
 		typedef		SimplexWithVertices<Vertex>										Self;
 	
-		typedef		std::set<Vertex>												VertexContainer;
+		typedef		std::vector<Vertex>												VertexContainer;
 		typedef		std::list<Self>													Cycle;
 		
 		/// \name Constructors 
@@ -40,11 +38,13 @@
 			vertices_(s.vertices_)													{}
 		template<class Iterator>
 		SimplexWithVertices(Iterator bg, Iterator end):
-			vertices_(bg, end)														{}
+			vertices_(bg, end)														{ std::sort(vertices_.begin(), vertices_.end()); }
 		SimplexWithVertices(const VertexContainer& v):	
-			vertices_(v)															{}
-		SimplexWithVertices(Vertex v):	
-			vertices_()																{ vertices_.insert(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)															{}
 		/// @}
 		
 		/// \name Core 
@@ -55,9 +55,9 @@
 		
 		/// \name Vertex manipulation
 		/// @{
-		bool					contains(const Vertex& v) const						{ return (vertices_.find(v) != vertices_.end()); }
+		bool					contains(const Vertex& v) const;
 		const VertexContainer&	vertices() const									{ return vertices_; }
-		void					add(const Vertex& v) 								{ vertices_.insert(v); }
+		void					add(const Vertex& v);
 		/// @}
 
 		/// \name Assignment and comparison
@@ -147,8 +147,8 @@
 			Parent(bg, end)															{}
 		SimplexWithAttachment(const Parent& s):
 			Parent(s)																{}
-		SimplexWithAttachment(VertexIndex vi):
-			Parent(vi), attachment(vi)												{}
+		SimplexWithAttachment(Dimension d, VertexIndex vi):
+			Parent(d, vi), attachment(vi)											{}
 		/// @}
 
 		void 					set_attachment(VertexIndex v)						{ attachment = v; }
--- a/include/topology/simplex.hpp	Tue Feb 19 06:51:23 2008 -0500
+++ b/include/topology/simplex.hpp	Tue Feb 19 06:56:49 2008 -0500
@@ -16,15 +16,30 @@
 
 	for (typename VertexContainer::const_iterator cur = vertices_.begin(); cur != vertices_.end(); ++cur)
 	{
-		bdry.push_back(*this);
+		bdry.push_back(Self(dimension() - 1));
 		Self& s = bdry.back();
-		s.vertices_.erase(*cur);
+		std::remove_copy(vertices_.begin(), vertices_.end(), s.vertices_.begin(), *cur);
 	}
 
 	return bdry;
 }
 
 template<class V>
+bool
+SimplexWithVertices<V>::
+contains(const Vertex& v) const
+{ 
+	typename VertexContainer::const_iterator location = std::lower_bound(vertices_.begin(), vertices_.end(), v); 
+	return ((location == vertices_.end()) || (*location == v)); 
+}
+		
+template<class V>
+void
+SimplexWithVertices<V>::
+add(const Vertex& v)
+{ vertices_.push_back(v); std::sort(vertices_.begin(), vertices_.end()); }
+
+template<class V>
 std::ostream&			
 SimplexWithVertices<V>::
 operator<<(std::ostream& out) const
--- a/include/topology/vineyard.hpp	Tue Feb 19 06:51:23 2008 -0500
+++ b/include/topology/vineyard.hpp	Tue Feb 19 06:56:49 2008 -0500
@@ -3,6 +3,12 @@
 #include <fstream>
 #include <sstream>
 
+#include "utilities/log.h"
+
+#ifdef LOGGING
+static rlog::RLogChannel* rlVineyard =			DEF_CHANNEL("topology/vineyard", rlog::Log_Debug);
+#endif // LOGGING
+
 template<class FS>
 void
 Vineyard<FS>::
@@ -50,7 +56,7 @@
 Vineyard<FS>::
 start_vine(Index i)
 {
-	Dout(dc::vineyard, "Starting new vine");
+	rLog(rlVineyard, "Starting new vine");
 	AssertMsg(i->sign(), "Can only start vines for positive simplices");
 		
 	Dimension dim = i->dimension();
@@ -65,7 +71,7 @@
 Vineyard<FS>::
 record_diagram(Index bg, Index end)
 {
-	Dout(dc::vineyard, "Entered: record_diagram()");
+	rLog(rlVineyard, "Entered: record_diagram()");
 	AssertMsg(evaluator != 0, "Cannot record diagram with a null evaluator");
 	
 	for (Index i = bg; i != end; ++i)
@@ -104,7 +110,7 @@
 Vineyard<FS>::
 record_knee(Index i)
 {
-	Dout(dc::vineyard, "Entered record_knee()");
+	rLog(rlVineyard, "Entered record_knee()");
 	AssertMsg(evaluator != 0, "Cannot record knee with a null evaluator");
 	AssertMsg(i->vine() != 0, "Cannot add a knee to a null vine");
 	AssertMsg(i->sign(), "record_knee() must be called on a positive simplex");
@@ -113,23 +119,23 @@
 		i->vine()->add(evaluator->value(*i), Infinity, evaluator->time());
 	else
 	{
-		Dout(dc::vineyard, "Creating knee");
+		rLog(rlVineyard, "Creating knee");
 		Knee k(evaluator->value(*i), evaluator->value(*(i->pair())), evaluator->time());
-		Dout(dc::vineyard, "Knee created: " << k);
+		rLog(rlVineyard, "Knee created: %s", tostring(k).c_str());
 
 		if (!k.is_diagonal() || i->vine()->empty())			// non-diagonal k, or empty vine
 		{
-			Dout(dc::vineyard, "Extending a vine");
+			rLog(rlVineyard, "Extending a vine");
 			i->vine()->add(k);
 		}
 		else if (i->vine()->back().is_diagonal())			// last knee is diagonal
 		{
 			AssertMsg(i->vine()->size() == 1, "Only first knee may be diagonal for a live vine");
-			Dout(dc::vineyard, "Overwriting first diagonal knee");
+			rLog(rlVineyard, "Overwriting first diagonal knee");
 			i->vine()->back() = k;
 		} else												// finish this vine
 		{
-			Dout(dc::vineyard, "Finishing a vine");
+			rLog(rlVineyard, "Finishing a vine");
 			i->vine()->add(k);
 			start_vine(i);
 			i->vine()->add(k);
@@ -137,7 +143,7 @@
 	}
 	
 	i->vine()->back().set_cycle(resolve_cycle(i));
-	Dout(dc::vineyard, "Leaving record_knee()");
+	rLog(rlVineyard, "Leaving record_knee()");
 }
 
 template<class FS>
@@ -145,7 +151,7 @@
 Vineyard<FS>::
 resolve_cycle(Index i) const
 {
-	Dout(dc::vineyard, "Entered resolve_cycle");
+	rLog(rlVineyard, "Entered resolve_cycle");
 	const Cycle& cycle = i->cycle();
 	
 	// Resolve simplices
--- a/include/utilities/counter.h	Tue Feb 19 06:51:23 2008 -0500
+++ b/include/utilities/counter.h	Tue Feb 19 06:56:49 2008 -0500
@@ -1,68 +1,80 @@
 /*
  * Author: Dmitriy Morozov
- * Department of Computer Science, Duke University, 2005 -- 2006
+ * Department of Computer Science, Duke University, 2005 -- 2007
  */
 
 #ifndef __COUNTER_H__
 #define __COUNTER_H__
 
-/* This should be integrated with the logging facilities. Written in the same framework. */
+
+#ifndef COUNTERS
+	#define 	GetCounter(path) 		0
+	#define 	Count(x)
+	#define		CountNum(x,y)
+	#define		CountBy(x,y)
+	#define		SetFrequency(x, freq)
+	#define		SetTrigger(x, y)
+#else // COUNTERS
+	#define 	GetCounter(path) 			get_counter(path)
+	#define 	Count(x) 					do { x->count++; if ((x->count % x->frequency == 0)) x->trigger->print(); } while (0)
+	#define 	CountNum(x,y) 				do { x->subcount[y]++; } while (0)
+	#define 	CountBy(x,y) 				do { x->count += y; } while (0)
+	#define		SetFrequency(x, freq)		do { x->frequency = freq; } while (0)
+	#define		SetTrigger(x, y)			do { x->trigger = y; } while(0)
+#endif // COUNTERS
+
 
 #include <map>
 #include <string>
 #include <iostream>
-
-#ifndef COUNTERS
-#define Count(x)
-#define CountNum(x,y)
-#else // COUNTERS
-#define Count(x) counters.inc(x)
-#define CountNum(x,y) counters.inc(x,y)
-#endif
+#include <limits>
+#include <unistd.h>
 
-typedef 		unsigned long 		CounterType;
-typedef			std::string			KeyType;
-
-class CounterFactory
+class Counter
 {
-	private:
-		typedef				std::map<int, CounterType> 			CountersMap;
-		typedef				std::map<KeyType, CountersMap>		KeyMap;
-		KeyMap				ctrs;
-		
+	public:
+		typedef 				unsigned long 							CounterType;
+		typedef					std::map<std::string, Counter*>			SubCounterMap;
+		typedef					std::map<int, CounterType>				SubCountMap;
+
+	public:
+		CounterType				count;
+		CounterType				frequency;
+		SubCountMap				subcount;
+		Counter*				trigger;
+
 	public:
-		~CounterFactory()
-		{
-#ifdef COUNTERS
-			std::cout << "Counters:" << std::endl;
-			for (KeyMap::const_iterator cur = ctrs.begin(); cur != ctrs.end(); ++cur)
-			{
-				std::cout << cur->first << ": ";
-				if (cur->second.size() == 1)
-				{
-					std::cout << cur->second.begin()->second << std::endl;
-					continue;
-				}
-				std::cout << std::endl;
-				for (CountersMap::const_iterator ccur = cur->second.begin();
-														  ccur != cur->second.end();
-														  ++ccur)
-					std::cout << "  " << ccur->first << ": " << ccur->second << std::endl;
-			}
-#endif // COUNTERS
-		}
+								Counter(const std::string& full_name = "",
+										CounterType freq = std::numeric_limits<CounterType>::max());
+								~Counter();
+
+		Counter*				get_child(const std::string& path, std::string::size_type pos);
+		void					print();
 
-		void inc(const KeyType& key, int num = 0)
-		{
-			ctrs[key][num]++;
-		}
+	private:	
+		SubCounterMap			subcounters_;
+		std::string				full_name_;
+		
+		static const char*		start_color;
+		static const char*		finish_color;
+		static const char		green_color[];
+		static const char 		normal_color[];
+		static const char 		empty_string[];
+};
+const char Counter::green_color[] 		= "\033[32m";
+const char Counter::normal_color[] 		= "\033[0m";
+const char Counter::empty_string[] 		= "";
+const char* Counter::start_color 		= 0;
+const char* Counter::finish_color 		= 0;
 
-		CounterType lookup(const KeyType& key, int num = 0) const
-		{
-			return const_cast<KeyMap&>(ctrs)[key][num];
-		}
-};
+static		Counter				rootCounter;
 
-static CounterFactory counters;
+Counter*	get_counter(const char* path)
+{
+	return rootCounter.get_child(path, 0);
+}
+
+#include "counter.hpp"
+
 
 #endif // __COUNTER_H__
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/utilities/counter.hpp	Tue Feb 19 06:56:49 2008 -0500
@@ -0,0 +1,69 @@
+#include <ctime>
+#include <cstdio>
+
+Counter::
+Counter(const std::string& full_name,
+		CounterType freq):
+		count(0), frequency(freq), trigger(this), full_name_(full_name)
+{ 
+	if (isatty(STDOUT_FILENO)) 
+	{
+		start_color = green_color; 
+		finish_color = normal_color;
+	} 
+	else
+	{	start_color = finish_color = empty_string;	}
+}
+
+
+Counter*
+Counter::
+get_child(const std::string& path, std::string::size_type pos)
+{
+	if (pos >= path.size())
+		return this;
+
+	std::string::size_type slash_pos = path.find('/', pos);
+	if (slash_pos == std::string::npos)
+		slash_pos = path.size();
+
+	std::string child_name = path.substr(pos, slash_pos - pos);
+	SubCounterMap::iterator child = subcounters_.find(child_name);
+
+	if (child != subcounters_.end())
+		return child->second->get_child(path, slash_pos + 1);
+	else	
+		return (subcounters_[child_name] = new Counter(path.substr(0, slash_pos)))->get_child(path, slash_pos + 1);
+}
+
+Counter::
+~Counter()
+{ 
+	if (full_name_ == "" && (!subcounters_.empty() || count))
+		print(); 
+
+	for (SubCounterMap::iterator cur = subcounters_.begin(); cur != subcounters_.end(); ++cur)
+		delete cur->second;
+}
+
+inline
+void
+Counter::
+print()
+{
+	time_t rawtime; time(&rawtime);
+	struct tm* timeinfo = localtime(&rawtime);
+
+	printf("%s(%02i:%02i:%02i)%s ",
+		    start_color,	
+			timeinfo->tm_hour,
+			timeinfo->tm_min,
+			timeinfo->tm_sec,
+			finish_color);
+	std::cout << "Counter [" << full_name_ << "]: " << count << std::endl;
+	for (SubCountMap::const_iterator cur = subcount.begin(); cur != subcount.end(); ++cur)
+		std::cout << "    " << cur->first << ": " << cur->second << std::endl;
+	for (SubCounterMap::iterator cur = subcounters_.begin(); cur != subcounters_.end(); ++cur)
+		cur->second->print();
+}	
+
--- a/include/utilities/debug.h	Tue Feb 19 06:51:23 2008 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,90 +0,0 @@
-#ifndef DEBUG_H
-#define DEBUG_H
-
-#ifndef CWDEBUG
-
-#include <iostream>     // std::cerr
-#include <cstdlib>      // std::exit, EXIT_FAILURE
-#include <cassert>
-
-#define AllocTag1(p)
-#define AllocTag2(p, desc)
-#define AllocTag_dynamic_description(p, data)
-#define AllocTag(p, data)
-#define Debug(STATEMENT...)
-#define Dout(cntrl, data)
-#define DoutFatal(cntrl, data) LibcwDoutFatal(, , cntrl, data)
-#define ForAllDebugChannels(STATEMENT...)
-#define ForAllDebugObjects(STATEMENT...)
-#define LibcwDebug(dc_namespace, STATEMENT...)
-#define LibcwDout(dc_namespace, d, cntrl, data)
-#define LibcwDoutFatal(dc_namespace, d, cntrl, data) do { ::std::cerr << data << ::std::endl; ::std::exit(EXIT_FAILURE); } while(1)
-#define LibcwdForAllDebugChannels(dc_namespace, STATEMENT...)
-#define LibcwdForAllDebugObjects(dc_namespace, STATEMENT...)
-#define NEW(x) new x
-#define CWDEBUG_ALLOC 0
-#define CWDEBUG_MAGIC 0
-#define CWDEBUG_LOCATION 0
-#define CWDEBUG_LIBBFD 0
-#define CWDEBUG_DEBUG 0
-#define CWDEBUG_DEBUGOUTPUT 0
-#define CWDEBUG_DEBUGM 0
-#define CWDEBUG_DEBUGT 0
-#define CWDEBUG_MARKER 0
-#define AssertMsg(TEST,MSG)
-//#define AssertMsg(TEST,STRM,MSG)
-
-
-#else // CWDEBUG
-
-// This must be defined before <libcwd/debug.h> is included and must be the
-// name of the namespace containing your `dc' (Debug Channels) namespace
-// (see below).  You can use any namespace(s) you like, except existing
-// namespaces (like ::, ::std and ::libcwd).
-#define DEBUGCHANNELS ::dionysus::debug::channels
-#include <libcwd/debug.h>
-
-namespace dionysus
-{
-	namespace debug 
-	{
-		void init(void);		// Initialize debugging code from main().
-		void init_thread(void);	// Initialize debugging code from new threads.
-
-		namespace channels 	// This is the DEBUGCHANNELS namespace, see above.
-		{
-			namespace dc 		// 'dc' is defined inside DEBUGCHANNELS.
-			{
-				using namespace libcwd::channels::dc;
-				using libcwd::channel_ct;
-
-				// Add the declaration of new debug channels here
-				// and add their definition in a custom debug.cc file.
-				extern channel_ct filtration;
-				extern channel_ct transpositions;
-				extern channel_ct vineyard;
-				extern channel_ct cycle;
-				extern channel_ct lsfiltration;
-				extern channel_ct orderlist;
-			} // namespace dc
-		} // namespace DEBUGCHANNELS
-	}
-}
-
-
-#define AssertMsg(TEST,MSG)                                           \
-				 ( (TEST) ? (void)0                                   \
-						  : (std::cerr << __FILE__ ":" << __LINE__    \
-								  << ": Assertion failed " #TEST      \
-								  << " - " << MSG << std::endl,abort()))
-/*
-#define AssertMsg(TEST,STRM,MSG)                                      \
-				 ( (TEST) ? (void)0                                   \
-						  : (DoutFatal(STRM, __FILE__ "(" << __LINE__ \
-								  << "): Assertion failed " #TEST     \
-								  << MSG << std::endl)))
-*/
-								  
-#endif // CWDEBUG
-
-#endif // DEBUG_H
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/utilities/log.h	Tue Feb 19 06:56:49 2008 -0500
@@ -0,0 +1,63 @@
+#ifndef __LOG_H__
+#define __LOG_H__
+
+#if LOGGING
+
+#define RLOG_COMPONENT dionysus
+
+#include <rlog/rlog.h>
+#include <rlog/RLogChannel.h>
+#include <rlog/StdioNode.h>
+#include <sstream>
+
+template<class T>
+std::string tostring(const T& t) { std::ostringstream out; out << t; return out.str(); }
+
+#define AssertMsg(cond, message, ...)		do { if (!(cond)) { rError(message, ##__VA_ARGS__); rAssertSilent(cond); } } while (0)
+	
+#else // LOGGING
+
+#include <unistd.h>		// for STDOUT_FILENO and STDERR_FILENO
+
+#define rDebug(...)
+#define rInfo(...)
+#define rWarning(...)
+#define rError(...)
+#define rLog(...)
+
+#define rAssert(...)
+#define rAssertSilent(...)
+
+#define DEF_CHANNEL(...) 0
+#define RLOG_CHANNEL(...) 0
+
+#define AssertMsg(cond, ...)
+
+// To avoid undefined errors for RLogChannel, we create a dummy namespace
+namespace rlog
+{
+	typedef		void			RLogChannel;
+
+	class StdioNode
+	{
+		public:
+								StdioNode(int,int)					{}
+			void				subscribeTo(RLogChannel*)			{}
+
+			static const int 	OutputColor = 0;
+			static const int	OutputChannel = 0;
+	};
+}
+
+#endif // LOGGING
+
+static rlog::StdioNode stdoutLog(STDOUT_FILENO,
+								 rlog::StdioNode::OutputColor + 
+								 rlog::StdioNode::OutputChannel);
+
+static rlog::StdioNode stderrLog(STDERR_FILENO,
+								 rlog::StdioNode::OutputColor + 
+								 rlog::StdioNode::OutputChannel);
+
+
+#endif //__LOG_H__
--- a/include/utilities/orderlist.h	Tue Feb 19 06:51:23 2008 -0500
+++ b/include/utilities/orderlist.h	Tue Feb 19 06:56:49 2008 -0500
@@ -12,15 +12,13 @@
 #ifndef __ORDERLIST_H__
 #define __ORDERLIST_H__
 
-#include "sys.h"
-#include "debug.h"
+#include "log.h"
 
 #include <iterator>
 #include <iostream>
 #include <list>
 
 #include "types.h"
-//#include "counter.h"
 
 #include <boost/utility.hpp>
 #include <boost/iterator/iterator_adaptor.hpp>
--- a/include/utilities/orderlist.hpp	Tue Feb 19 06:51:23 2008 -0500
+++ b/include/utilities/orderlist.hpp	Tue Feb 19 06:56:49 2008 -0500
@@ -1,5 +1,9 @@
 /* Implementations */
 
+#ifdef LOGGING
+static rlog::RLogChannel* rlOrderList = 					DEF_CHANNEL("utilities/orderlist", rlog::Log_Debug);
+#endif // LOGGING
+
 template<class T>
 void 
 OrderList<T>::
@@ -70,16 +74,16 @@
 	} while (inv_density * num_elements >= maximum);
 	++num_elements;			// for the extra element inserted
 
-	Dout(dc::orderlist, num_elements << ", " << lower << ", " << upper);
-	Dout(dc::orderlist, "prev is at the end: " << (prev == Parent::end()));
-	Dout(dc::orderlist, "next is at the end: " << (next == Parent::end()));
+	rLog(rlOrderList, "%i, %i, %i", num_elements, lower, upper);
+	rLog(rlOrderList, "prev is at the end: %i", (prev == Parent::end()));
+	rLog(rlOrderList, "next is at the end: %i", (next == Parent::end()));
 	
 	// Reorder
 	AssertMsg((upper - lower + 1)/num_elements > 0, "Spacing between new tags must be non-zero");
 	for (unsigned int i = 0; i < num_elements; ++i)
 	{
 		(++prev)->tag = lower + i*((upper - lower + 1)/num_elements);
-		Dout(dc::orderlist, prev->tag);
+		rLog(rlOrderList, "%i", prev->tag);
 		AssertMsg(prev->tag != 0 || prev == Parent::begin(), "Cannot assign 0 tag except at the beginning of OrderList");
 	}
 
--- a/include/utilities/sys.h	Tue Feb 19 06:51:23 2008 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,26 +0,0 @@
-// sys.h
-//
-// This header file is included at the top of every source file,
-// before any other header file.
-// It is intended to add defines that are needed globally and
-// to work around Operating System dependend incompatibilities.
-
-// EXAMPLE: If you use autoconf you can add the following here.
-// #ifdef HAVE_CONFIG_H
-// #include "config.h"
-// #endif
-
-// EXAMPLE: You could add stuff like this here too
-// (Otherwise add -DCWDEBUG to your CFLAGS).
-// #if defined(WANTSDEBUGGING) && defined(HAVE_LIBCWD_BLAHBLAH)
-// #define CWDEBUG
-// #endif
-
-// The following is the libcwd related mandatory part.
-// It must be included before any system header file is included!
-#ifdef CWDEBUG
-#ifndef _GNU_SOURCE
-#define _GNU_SOURCE
-#endif
-#include <libcwd/sys.h>
-#endif
--- a/include/utilities/types.h	Tue Feb 19 06:51:23 2008 -0500
+++ b/include/utilities/types.h	Tue Feb 19 06:56:49 2008 -0500
@@ -5,7 +5,7 @@
 
 /* Types */
 typedef 	bool					Sign;
-typedef		short int				Dimension;
+typedef		unsigned short int		Dimension;
 const 		Sign	 				POS = true;
 const 		Sign					NEG = false;
 typedef		double					RealType;
--- a/src/debug.cpp	Tue Feb 19 06:51:23 2008 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,81 +0,0 @@
-#include "sys.h"
-#include "debug.h"
-
-#ifdef CWDEBUG
-
-namespace dionysus
-{
-	namespace debug 
-	{
-		namespace channels 		// DEBUGCHANNELS is set to this in debug.h
-		{
-			namespace dc 
-			{
-		        // Add new debug channels here.
-				channel_ct filtration("FILTRATION");
-				channel_ct transpositions("TRANSPOS");
-				channel_ct vineyard("VINEYARD");
-				channel_ct cycle("CYCLE");
-				channel_ct lsfiltration("LSFILTRATION");
-				channel_ct orderlist("ORDERLIST");
-			}
-    	} // namespace DEBUGCHANNELS
-
-		// Initialize debugging code from new threads.
-		void init_thread(void)
-		{
-			// Everything below needs to be repeated at the start of every
-			// thread function, because every thread starts in a completely
-			// reset state with all debug channels off etc.
-		
-			#if LIBCWD_THREAD_SAFE		// For the non-threaded case this is set by the rcfile.
-			// Turn on all debug channels by default.
-			ForAllDebugChannels(while(!debugChannel.is_on()) debugChannel.on());
-			// Turn off specific debug channels.
-			Debug(dc::bfd.off());
-			Debug(dc::malloc.off());
-			#endif
-		
-			// Turn on debug output.
-			// Only turn on debug output when the environment variable SUPPRESS_DEBUG_OUTPUT is not set.
-			Debug(if (getenv("SUPPRESS_DEBUG_OUTPUT") == NULL) libcw_do.on());
-			#if LIBCWD_THREAD_SAFE
-			Debug(libcw_do.set_ostream(&std::cout, &cout_mutex));
-		
-			// Set the thread id in the margin.
-			char margin[12];
-			sprintf(margin, "%-10lu ", pthread_self());
-			Debug(libcw_do.margin().assign(margin, 11));
-			#else
-			Debug(libcw_do.set_ostream(&std::cout));
-			#endif
-		
-			// Write a list of all existing debug channels to the default debug device.
-			Debug(list_channels_on(libcw_do));
-		}
-		
-		// Initialize debugging code from main().
-		void init(void)
-		{
-			// You want this, unless you mix streams output with C output.
-			// Read  http://gcc.gnu.org/onlinedocs/libstdc++/27_io/howto.html#8 for an explanation.
-			// We can't use it, because other code uses printf to write to the console.
-			//std::ios::sync_with_stdio(false);
-
-			// This will warn you when you are using header files that do not belong to the
-			// shared libcwd object that you linked with.
-			Debug(check_configuration());
-		
-			#if CWDEBUG_ALLOC
-			// Remove all current (pre- main) allocations from the Allocated Memory Overview.
-			libcwd::make_all_allocations_invisible_except(NULL);
-			#endif
-			
-			//Debug(read_rcfile());
-			
-			init_thread();
-		}
-	} // namespace debug
-} // namespace dionysus
-
-#endif // CWDEBUG
--- a/tests/utilities/CMakeLists.txt	Tue Feb 19 06:51:23 2008 -0500
+++ b/tests/utilities/CMakeLists.txt	Tue Feb 19 06:56:49 2008 -0500
@@ -1,8 +1,9 @@
 set							(targets						
 							 test-consistencylist
-							 test-orderlist)
+							 test-orderlist
+							 test-counters)
 
 foreach 					(t ${targets})
-	add_executable			(${t} ${t}.cpp ${external_sources})
+	add_executable			(${t} ${t}.cpp)
 	target_link_libraries	(${t} ${libraries})
 endforeach 					(t ${targets})
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/utilities/test-counters.cpp	Tue Feb 19 06:56:49 2008 -0500
@@ -0,0 +1,33 @@
+#include <utilities/counter.h>
+
+static Counter* cTestElaborate = GetCounter("test/elaborate");
+static Counter* cTestBasic = GetCounter("test/basic");
+static Counter* cTestBasicSub = GetCounter("test/basic/sub");
+
+int main()
+{
+	SetFrequency(cTestBasic, 2);
+
+	Count(cTestBasic);
+	Count(cTestBasicSub);
+	Count(cTestBasicSub);
+	Count(cTestBasicSub);
+	Count(cTestElaborate);
+	Count(cTestBasic);
+	Count(cTestElaborate);
+	Count(cTestBasic);
+	Count(cTestBasic);
+	CountNum(cTestBasic, 25);
+	CountNum(cTestBasic, 132);
+	CountNum(cTestBasic, 25);
+	CountNum(cTestBasic, 121);
+	CountNum(cTestBasic, 132);
+	CountNum(cTestBasic, 25);
+	
+	SetTrigger(cTestBasic, &rootCounter);
+	Count(cTestBasic);
+	Count(cTestBasic);
+	
+	SetFrequency(cTestElaborate, 3);
+	Count(cTestElaborate);
+}