Merge in AR branch dev
authorDmitriy Morozov <morozov@cs.duke.edu>
Fri, 04 Apr 2008 21:14:44 -0400
branchdev
changeset 87 2c2e2f3b5d15
parent 54 8f88bd32153a (current diff)
parent 86 73a54447b54a (diff)
child 88 9e00cf44ea5d
Merge in AR branch
CMakeLists.txt
examples/ar-vineyard/CMakeLists.txt
examples/ar-vineyard/ar-vineyard.h
examples/ar-vineyard/ar-vineyard.hpp
include/geometry/euclidean.h
include/geometry/polynomial.h
include/geometry/polynomial.hpp
include/geometry/simulator.h
include/geometry/simulator.hpp
include/topology/filtration.hpp
include/topology/filtrationcontainer.h
include/topology/filtrationsimplex.h
include/topology/vineyard.h
include/topology/vineyard.hpp
include/utilities/eventqueue.h
tests/geometry/CMakeLists.txt
tests/geometry/Makefile
tests/geometry/polynomial.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/89ae955518665a61	Fri Apr 04 21:14:44 2008 -0400
@@ -0,0 +1,16 @@
+From artemis Tue Feb 26 23:24:44 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Tue, 26 Feb 2008 18:22:06 -0500
+State: new
+Subject: Get rid of intostring() and .c_str()
+Message-Id: <89ae955518665a61-0-artemis@metatron>
+
+Get rid of the need for intostring() (in addition to tostring()), and having to
+place .c_str() after tostring() in rLog calls.
+
+The former is necessary because of some problem with disambiguating which
+operator<<(ostream,T) to use when Event is being output, so intostring() calls
+T.operator<<(ostream) explicitly. This problem seems to exist only for Events.
+
+It should be possible to solve the latter by returning char* from tostring()
+rather than std::string.
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/b5d8a8403ae3a0d2	Fri Apr 04 21:14:44 2008 -0400
@@ -0,0 +1,18 @@
+From artemis Thu Feb 21 09:54:54 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Thu, 21 Feb 2008 04:54:00 -0500
+State: fixed
+Subject: Segfault in tests/geometry/test-kinetic-sort.cpp
+Message-Id: <b5d8a8403ae3a0d2-0-artemis@metatron>
+
+test-kinetic-sort segfaults when we subscribe to geometry/simulator RLog channel.
+
+From artemis Fri Feb 22 18:05:55 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Fri, 22 Feb 2008 13:05:55 -0500
+Subject: properties changes (state)
+Message-Id: <b5d8a8403ae3a0d2-a0288cca74095157-artemis@metatron>
+References: <b5d8a8403ae3a0d2-0-artemis@metatron>
+In-Reply-To: <b5d8a8403ae3a0d2-0-artemis@metatron>
+
+state=fixed
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/c420501cc5285bbc	Fri Apr 04 21:14:44 2008 -0400
@@ -0,0 +1,71 @@
+From artemis Tue Feb 26 10:25:07 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Tue, 26 Feb 2008 05:22:56 -0500
+State: fixed
+Subject: Non-optimized CGAL runtime error
+Message-Id: <c420501cc5285bbc-0-artemis@metatron>
+
+If the code is compiled with optimizations off or with debug on, CGAL gives a
+runtime error.
+
+From artemis Tue Feb 26 13:00:23 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Tue, 26 Feb 2008 08:00:23 -0500
+Subject: properties changes (state)
+Message-Id: <c420501cc5285bbc-e983173b6cd399a6-artemis@metatron>
+References: <c420501cc5285bbc-0-artemis@metatron>
+In-Reply-To: <c420501cc5285bbc-0-artemis@metatron>
+
+state=fixed
+
+From artemis Tue Feb 26 13:01:57 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Tue, 26 Feb 2008 08:00:50 -0500
+Subject: Fixed by adding CGAL_NO_ASSERTIONS
+Message-Id: <c420501cc5285bbc-d75d85d67d421f7e-artemis@metatron>
+References: <c420501cc5285bbc-0-artemis@metatron>
+In-Reply-To: <c420501cc5285bbc-0-artemis@metatron>
+
+I suspect the bug in CGAL (since everything works fine with CGAL_NO_ASSERTIONS
+set), so I disabled the offending assertion (and all the rest of them) by
+setting a CXX flag in CMakeLists.txt.
+
+From artemis Wed Mar 19 16:47:55 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Wed, 19 Mar 2008 12:47:24 -0400
+Subject: Fixed in 0a18d6902a55
+Message-Id: <c420501cc5285bbc-db67b15857938f7b-artemis@metatron>
+References: <c420501cc5285bbc-0-artemis@metatron>
+In-Reply-To: <c420501cc5285bbc-0-artemis@metatron>
+
+Fixed by dealing with infinite simplices in 0a18d6902a55
+
+From artemis Wed Mar 19 16:49:36 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Wed, 19 Mar 2008 12:49:01 -0400
+Subject: Reopened for ar-vineyard
+Message-Id: <c420501cc5285bbc-8f74ccde6b3f0bea-artemis@metatron>
+References: <c420501cc5285bbc-0-artemis@metatron>
+In-Reply-To: <c420501cc5285bbc-0-artemis@metatron>
+
+Fix for ar-vineyard like for alphashapes
+
+From artemis Wed Mar 19 16:49:47 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Wed, 19 Mar 2008 12:49:47 -0400
+Subject: properties changes (state)
+Message-Id: <c420501cc5285bbc-ef5f506b45aba335-artemis@metatron>
+References: <c420501cc5285bbc-0-artemis@metatron>
+In-Reply-To: <c420501cc5285bbc-0-artemis@metatron>
+
+state=open
+
+From artemis Wed Mar 19 17:17:47 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Wed, 19 Mar 2008 13:17:30 -0400
+Subject: Fixed for ar-vineyard
+Message-Id: <c420501cc5285bbc-feb676745fc16b8c-artemis@metatron>
+References: <c420501cc5285bbc-8f74ccde6b3f0bea-artemis@metatron>
+In-Reply-To: <c420501cc5285bbc-8f74ccde6b3f0bea-artemis@metatron>
+
+Fixed for ar-vineyard.
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/c664b2f69b5f6ea3	Fri Apr 04 21:14:44 2008 -0400
@@ -0,0 +1,29 @@
+From artemis Thu Feb 28 10:34:23 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Thu, 28 Feb 2008 05:33:34 -0500
+State: fixed
+Subject: Two simulators
+Message-Id: <c664b2f69b5f6ea3-0-artemis@metatron>
+
+Consider using two simulators instead of one to make sure that trajectory
+changes are processed before simplex swaps.
+
+From artemis Sat Mar  1 09:45:05 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Sat, 01 Mar 2008 04:45:05 -0500
+Subject: properties changes (state)
+Message-Id: <c664b2f69b5f6ea3-f8bf16bdff01098d-artemis@metatron>
+References: <c664b2f69b5f6ea3-0-artemis@metatron>
+In-Reply-To: <c664b2f69b5f6ea3-0-artemis@metatron>
+
+state=fixed
+
+From artemis Sat Mar  1 09:46:20 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Sat, 01 Mar 2008 04:45:29 -0500
+Subject: Fixed in f236c7d659d0
+Message-Id: <c664b2f69b5f6ea3-a1742f0eb7b5e1c0-artemis@metatron>
+References: <c664b2f69b5f6ea3-0-artemis@metatron>
+In-Reply-To: <c664b2f69b5f6ea3-0-artemis@metatron>
+
+Fixed in f236c7d659d0.
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/cd79223a108d3900	Fri Apr 04 21:14:44 2008 -0400
@@ -0,0 +1,9 @@
+From artemis Fri Feb 22 23:41:38 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Fri, 22 Feb 2008 18:39:31 -0500
+State: new
+Subject: Incorrect sign_at() in UPolynomial<double>
+Message-Id: <cd79223a108d3900-0-artemis@metatron>
+
+UPolynomial<double>::sign_at() reports -1 as the sign of -2*x + 4 at x = 2. 
+Should be 0.
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/d2ab07329c3588ca	Fri Apr 04 21:14:44 2008 -0400
@@ -0,0 +1,29 @@
+From artemis Wed Feb 27 21:31:41 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Wed, 27 Feb 2008 16:30:57 -0500
+State: fixed
+Subject: Segfault with topology/vineyard log enabled
+Message-Id: <d2ab07329c3588ca-0-artemis@metatron>
+
+The code segfaults in Vineyard::record_knee() if we subscribe to
+topology/vineyard log.
+
+From artemis Sat Mar  1 09:46:46 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Sat, 01 Mar 2008 04:46:46 -0500
+Subject: properties changes (state)
+Message-Id: <d2ab07329c3588ca-8c184b421f538482-artemis@metatron>
+References: <d2ab07329c3588ca-0-artemis@metatron>
+In-Reply-To: <d2ab07329c3588ca-0-artemis@metatron>
+
+state=fixed
+
+From artemis Sat Mar  1 09:47:07 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Sat, 01 Mar 2008 04:46:51 -0500
+Subject: Fixed in abba2950aced
+Message-Id: <d2ab07329c3588ca-fa65c35dfff6366f-artemis@metatron>
+References: <d2ab07329c3588ca-0-artemis@metatron>
+In-Reply-To: <d2ab07329c3588ca-0-artemis@metatron>
+
+Fixed in abba2950aced
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/ebda8db3f9908e33	Fri Apr 04 21:14:44 2008 -0400
@@ -0,0 +1,8 @@
+From artemis Mon Feb 25 16:30:10 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Mon, 25 Feb 2008 11:29:27 -0500
+State: new
+Subject: Efficient EventQueue
+Message-Id: <ebda8db3f9908e33-0-artemis@metatron>
+
+Change EventQueue to an efficient implementation, e.g., using a Fibonacci heap.
--- a/CMakeLists.txt	Fri Apr 04 16:17:28 2008 -0400
+++ b/CMakeLists.txt	Fri Apr 04 21:14:44 2008 -0400
@@ -5,6 +5,7 @@
 option						(debug				"Build Dionysus with debugging on" 		OFF)
 option						(optimize			"Build Dionysus with optimization"		ON)
 option						(use_dsrpdb			"Build examples that use DSR-PDB"		OFF)
+option						(use_synaps			"Build examples that use SYNAPS"		OFF)
 
 # Find everything that's always required
 find_package				(Boost REQUIRED)
@@ -43,6 +44,15 @@
 														${gmpxx_LIBRARY} 
 														${m_LIBRARY})
 
+# SYNAPS
+if                          (use_synaps)
+    add_definitions			(-DBOOST_UBLAS_TYPE_CHECK=0)
+    find_library			(synaps_LIBRARY				NAMES synaps)
+    set						(synaps_libraries			${synaps_LIBRARY}
+														${gmp_LIBRARY}
+														${gmpxx_LIBRARY})
+endif                       (use_synaps)
+
 # Debugging
 if							(debug)
 	if 						(optimize)
--- a/examples/ar-vineyard/CMakeLists.txt	Fri Apr 04 16:17:28 2008 -0400
+++ b/examples/ar-vineyard/CMakeLists.txt	Fri Apr 04 21:14:44 2008 -0400
@@ -1,8 +1,15 @@
 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)
-	target_link_libraries	(${t} ${libraries} ${cgal_libraries})
+	add_executable			(${t} ${t}.cpp ${external_sources})
+	target_link_libraries	(${t}   ${libraries} ${cgal_libraries} 
+                                    ${boost_signals_LIBRARY} 
+                                    ${boost_program_options_LIBRARY})
 endforeach 					(t)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/ar-vineyard/ar-function-kernel.h	Fri Apr 04 21:14:44 2008 -0400
@@ -0,0 +1,80 @@
+#ifndef __AR_FUNCTION_KERNEL_H__
+#define __AR_FUNCTION_KERNEL_H__
+
+
+#include <stack>
+#include <iostream>
+#include "ar-simplex3d.h"
+
+
+
+/**
+ * Represents function suitable for the FunctionKernel. Albeit very restrictive 
+ * (it only supports subtraction), although that should suffice for KineticSort.
+ */
+class ARFunction
+{
+    public:
+        /// Represents the three forms of the involved functions. See Appendix B of LHI paper.
+        enum                        FunctionForm                                        { none, rho, lambda, phi };
+
+    public:
+                                    ARFunction(FunctionForm f, ARSimplex3D* s): 
+                                        form_(f), form2_(none),
+                                        simplex_(s), simplex2_(0)                       {}
+                                    ARFunction(const ARFunction& other):
+                                        form_(other.form_), form2_(other.form2_),
+                                        simplex_(other.simplex_), 
+                                        simplex2_(other.simplex2_)                      {}
+
+        ARFunction&                 operator-=(const ARFunction& other)                 { form2_ = other.form_; simplex2_ = other.simplex_; return *this; }
+        ARFunction                  operator-(const ARFunction& other)                  { return (ARFunction(*this) -= other); }
+
+        FunctionForm                form() const                                        { return form_; }
+        FunctionForm                form2() const                                       { return form2_; }
+        ARSimplex3D*                simplex() const                                     { return simplex_; }
+        ARSimplex3D*                simplex2() const                                    { return simplex2_; }
+
+        std::ostream&               operator<<(std::ostream& out) const                 { return (out << "(" << form_ << ", " << simplex_ << "), (" 
+                                                                                                             << form2_ << ", " << simplex2_ << ")"); }
+
+    private:
+        FunctionForm                form_, form2_;      // the function is form_ - form2_
+        ARSimplex3D                 *simplex_, *simplex2_;
+};
+
+std::ostream&
+operator<<(std::ostream& out, const ARFunction& f)
+{ return f.operator<<(out); }
+
+/**
+ * Function kernel specialized at solving the kinds of functions involved in 
+ * ARVineyard construction. We cannot use a polynomial kernel (which includes 
+ * rational functions) that we already have because not all involved functions
+ * are rational.
+ */
+class ARFunctionKernel
+{
+    public:
+        typedef                     double                                              FieldType;
+        typedef                     FieldType                                           RootType;
+        typedef                     std::stack<RootType>                                RootStack;
+        typedef                     ARSimplex3D::RealValue                              SimplexFieldType;
+
+        typedef                     ARFunction                                          Function;
+        typedef                     Function::FunctionForm                              FunctionForm;
+    
+    public:
+        static void                 solve(const Function& f, RootStack& stack);
+        static RootType             root(RootType r)                                    { return r; }
+        static RootType             root(SimplexFieldType r)                            { return CGAL::to_double(r); }
+        static int                  sign_at(const Function& f, RootType r);
+        static RootType             between(RootType r1, RootType r2)                   { return (r1 + r2)/2; }
+        static int                  sign_at_negative_infinity(const Function& f);
+
+        static FieldType            value_at(const Function& f, RootType v);
+};
+
+#include "ar-function-kernel.hpp"
+
+#endif
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/ar-vineyard/ar-function-kernel.hpp	Fri Apr 04 21:14:44 2008 -0400
@@ -0,0 +1,144 @@
+#include <utilities/log.h>
+#include <cmath>
+
+#if LOGGING
+static rlog::RLogChannel* rlARFunctionKernel =                      DEF_CHANNEL("ar/function-kernel", rlog::Log_Debug);
+static rlog::RLogChannel* rlARFunctionKernelValue =                 DEF_CHANNEL("ar/function-kernel/value", rlog::Log_Debug);
+#endif
+
+
+/* For now only rho and phi are implemented */
+void
+ARFunctionKernel::
+solve(const Function& f, RootStack& stack)
+{
+    AssertMsg(stack.empty(), "Stack must be empty before solving");
+    AssertMsg((f.form() != Function::none) && (f.form2() != Function::none), "Solving possible only for differences");
+
+    FunctionForm f1 = f.form(), f2 = f.form2();
+    const ARSimplex3D   *s1 = f.simplex(), 
+                        *s2 = f.simplex2();
+    if (f1 < f2)    { std::swap(f1,f2); std::swap(s1,s2); }         // for simplicity enforce function order 
+                                                                    // to handle fewer cases explicitly
+    AssertMsg(f1 != Function::lambda && f2 != Function::lambda, "Lambda not implemented yet");
+
+    rLog(rlARFunctionKernel,    "Solve: function1 = (%i, %p), function2 = (%i, %p)",
+                                f1, s1, f2, s2);
+    
+    //if (f1 == Function::phi && f2 == Function::phi)         return;
+    //if (f1 == Function::rho && f2 == Function::rho)         return;
+
+    if (f1 == Function::phi && f2 == Function::rho)
+    {
+        SimplexFieldType r = s2->alpha() - s1->phi_const();
+        rLog(rlARFunctionKernel, "  phi = rho => r^2 = %s (%lf)", tostring(r).c_str(), root(r));
+        stack.push(root(r));
+    }
+
+    if (f1 == Function::phi && f2 == Function::lambda)
+    {
+        rLog(rlARFunctionKernel, "  phi = lambda");
+        SimplexFieldType r2 = (s2->rho() + s2->v() - s2->s() - s1->phi_const());
+        r2 *= r2;
+        r2 /= 4*s2->v();
+        r2 += s2->s();
+        if (r2 >= s2->s() + s2->v())
+            stack.push(root(r2));
+
+        SimplexFieldType r1 = s2->rho() - s1->phi_const();
+        if (r1 <= s2->s() + s2->v())
+            stack.push(root(r1));
+    }
+    
+    // FIXME: this is far from complete!
+    if (f1 == Function::lambda && f2 == Function::lambda)
+    {
+        rLog(rlARFunctionKernel, "  lambda = lambda");
+        if ((s1->s() + s1->v() < s2->s() + s2->v()))                // let f1 be the function with larger break point
+        {   std::swap(f1,f2); std::swap(s1,s2); }
+
+        if (s1->rho() > s2->rho())
+        {
+            RootType r = root(s2->s() + s2->v() + s1->rho() - s2->rho()) + 2*sqrt(root(s2->v()*(s1->rho() - s2->rho())));
+            if (r < root(s1->s() + s1->v()) && r > root(s2->s() + s2->v()))
+                stack.push(r);
+        }
+    }
+
+    if (f1 == Function::lambda && f2 == Function::rho)
+    {
+        rLog(rlARFunctionKernel, "  lambda = rho");
+        // perhaps no solutions instead of an assertion is the right way to deal with this
+        AssertMsg(s2->alpha() > s1->rho(), "Rho_0^2 must be greater than Rho^2");
+
+        RootType r = sqrt(root(s2->v()*(s2->alpha() - s1->rho())));         // damn square roots
+        r *= 2;
+        r += root(s1->s() + s1->v() + s2->alpha() - s1->rho());
+    }
+    rLog(rlARFunctionKernel, "  Stack size: %i", stack.size());
+    if (stack.size() > 0) rLog(rlARFunctionKernel, "  Top: %lf", stack.top());
+}
+
+int
+ARFunctionKernel::
+sign_at(const Function& f, RootType r)
+{
+    FieldType v = value_at(f,r);
+    if (v > 0)  return true;
+    else        return false;
+}
+
+int
+ARFunctionKernel::
+sign_at_negative_infinity(const Function& f)
+{
+    FunctionForm f1 = f.form(), f2 = f.form2();
+    const ARSimplex3D   *s1 = f.simplex(), 
+                        *s2 = f.simplex2();
+    int multiplier = 1;
+    if (f1 < f2) { std::swap(f1, f2); std::swap(s1, s2); multiplier = -1; }
+    
+    AssertMsg(f1 != Function::lambda && f2 != Function::lambda, "Lambda not implemented yet");
+
+    if (f1 == Function::phi && f2 == Function::phi)
+    {
+        if (s1->phi_const() == s2->phi_const()) return  0;
+        if (s1->phi_const() > s2->phi_const())  return  1;        // multiplier must be 1
+        else                                    return -1;
+    }
+    
+    if (f1 == Function::phi && f2 == Function::rho)
+        return -multiplier;
+
+    if (f1 == Function::rho && f2 == Function::rho)
+    {
+        if (s1->alpha() == s2->alpha())         return  0;
+        if (s1->alpha() > s2->alpha())          return  1;        // multiplier must be 1
+        else                                    return -1;
+    }
+
+    AssertMsg(false, "The case analysis should be exhaustive in sign at -infinity");
+    return false;
+}
+        
+ARFunctionKernel::FieldType            
+ARFunctionKernel::
+value_at(const Function& f, RootType v)
+{
+    FunctionForm f1 = f.form(), f2 = f.form2();
+    ARSimplex3D     *s1 = f.simplex(), 
+                    *s2 = f.simplex2();
+
+    AssertMsg(f2 == Function::none && s2 == 0, "Value_at knows only about functions themselves, not their differences");
+    AssertMsg(f1 != Function::lambda, "Lambda not implemented yet");
+    rLog(rlARFunctionKernelValue,    "Value_at: function = (%i, %p)", f1, s1);
+
+    if (f1 == Function::phi)
+        return v + root(s1->phi_const());
+
+    if (f1 == Function::rho)
+        return root(s1->alpha());
+    
+    AssertMsg(false, "The case analysis should be exhaustive in value_at");
+    return 0;
+}
--- a/examples/ar-vineyard/ar-simplex3d.h	Fri Apr 04 16:17:28 2008 -0400
+++ b/examples/ar-vineyard/ar-simplex3d.h	Fri Apr 04 21:14:44 2008 -0400
@@ -28,7 +28,6 @@
 typedef Delaunay::Facet								Facet;
 typedef Delaunay::Cell								Cell;
 typedef Delaunay::Cell_handle						Cell_handle;
-typedef K::FT										RealValue;
 
 typedef Delaunay::Finite_vertices_iterator    		Vertex_iterator;
 typedef Delaunay::Finite_edges_iterator        		Edge_iterator;
@@ -40,6 +39,7 @@
 class ARSimplex3D: public SimplexWithVertices<Vertex_handle>
 {
 	public:
+        typedef     K::FT										        RealValue;
 		typedef 	std::map<ARSimplex3D, RealValue>					SimplexPhiMap;
 		typedef		SimplexWithVertices<Vertex_handle>					Parent;
 		typedef		Parent::VertexContainer								VertexSet;
@@ -56,10 +56,11 @@
 		
 								    ARSimplex3D(const Edge& e);
 								    ARSimplex3D(const Edge& e, const Point& z, SimplexPhiMap& simplices, 
-												Facet_circulator facet_bg);
+                                                const Delaunay& Dt, Facet_circulator facet_bg);
 		
 								    ARSimplex3D(const Facet& f);
-								    ARSimplex3D(const Facet& f, const Point& z, const SimplexPhiMap& simplices);
+								    ARSimplex3D(const Facet& f, const Point& z, const SimplexPhiMap& simplices,
+                                                const Delaunay& Dt);
 	    
 									ARSimplex3D(const Cell& c, const Point& z);
 	    
@@ -86,9 +87,11 @@
 		RealValue					rho_;		// rho_ 	is the squared radius of the smallest circumsphere
 		RealValue					s_;			// s_ 		is the squared distance from z to the affine hull of the simplex
 		RealValue					v_;			// v_ 		is the squared distance from z to the affine hull of the dual Voronoi cell
-		RealValue					phi_const_;	// see LHI paper, Appendix 1
+		RealValue					phi_const_;	// see LHI paper, Appendices A and B
 		bool 						attached_;
 
+        // in paper's notation: s_ = v^2; v_ = d^2
+
 };
 
 typedef 			std::vector<ARSimplex3D>								ARSimplex3DVector;
--- a/examples/ar-vineyard/ar-simplex3d.hpp	Fri Apr 04 16:17:28 2008 -0400
+++ b/examples/ar-vineyard/ar-simplex3d.hpp	Fri Apr 04 21:14:44 2008 -0400
@@ -1,3 +1,9 @@
+#include <utilities/log.h>
+
+#if LOGGING
+static rlog::RLogChannel* rlARSimplex3D =                       DEF_CHANNEL("ar/simplex3d", rlog::Log_Debug);
+#endif
+
 ARSimplex3D::	    
 ARSimplex3D(const ARSimplex3D& s): Parent(s)
 { 
@@ -35,19 +41,15 @@
 }
 
 ARSimplex3D::	    
-ARSimplex3D(const Edge& e, const Point& z, SimplexPhiMap& simplices, Facet_circulator facet_bg)
+ARSimplex3D(const Edge& e, const Point& z, SimplexPhiMap& simplices, const Delaunay& Dt, Facet_circulator facet_bg)
 {
     Cell_handle c = e.first;
 	Parent::add(c->vertex(e.second));
 	Parent::add(c->vertex(e.third));
 
 	Facet_circulator cur = facet_bg;
+	while (Dt.is_infinite(*cur))    ++cur; 
 	SimplexPhiMap::const_iterator cur_iter = simplices.find(ARSimplex3D(*cur));
-	while (cur_iter == simplices.end())
-	{
-		++cur; 
-		cur_iter = simplices.find(ARSimplex3D(*cur));
-	}
 	RealValue min = cur_iter->first.alpha();
 	RealValue phi_const_min = cur_iter->first.phi_const();
 	
@@ -63,10 +65,8 @@
 		int i0 = (*cur).first->index(*v1);
 		int i1 = (*cur).first->index(*v2);
 		int i = 6 - i0 - i1 - (*cur).second;
-		Point p3 = (*cur).first->vertex(i)->point();
 
-		cur_iter = simplices.find(ARSimplex3D(*cur));
-		if (cur_iter == simplices.end())			// cur is infinite
+		if (Dt.is_infinite(cur->first->vertex(i)))
 		{
 			++cur; continue;
 			// FIXME: what do we do with infinite cofaces (i.e., what
@@ -74,9 +74,11 @@
 			// infinite?
 		}
 		
+		Point p3 = (*cur).first->vertex(i)->point();
 		if (CGAL::side_of_bounded_sphere(p1, p2, p3) == CGAL::ON_BOUNDED_SIDE)
 			attached_ = true;
 		
+	    SimplexPhiMap::const_iterator cur_iter = simplices.find(ARSimplex3D(*cur));
 		RealValue 								val 				= cur_iter->first.alpha();
 		if (val < min) 							min 				= val;
 		RealValue 								phi_const_val 		= cur_iter->first.phi_const();
@@ -99,7 +101,7 @@
 
 	
 	s_ = CGAL::squared_distance(z, K::Segment_3(p1,p2).supporting_line());
-	Point origin;
+	Point origin(0,0,0);
 	Point cc = origin + ((p1 - origin) + (p2 - origin))/2;		// CGAL is funny
 	v_ = CGAL::squared_distance(z, cc) - s_;
 }
@@ -114,7 +116,7 @@
 }
 
 ARSimplex3D::	    
-ARSimplex3D(const Facet& f, const Point& z, const SimplexPhiMap& simplices)
+ARSimplex3D(const Facet& f, const Point& z, const SimplexPhiMap& simplices, const Delaunay& Dt)
 {
     Cell_handle c = f.first;
 	for (int i = 0; i < 4; ++i)
@@ -133,29 +135,33 @@
 	rho_ = squared_radius(p1, p2, p3);
 	
 	attached_ = false;
-	if (CGAL::side_of_bounded_sphere(p1, p2, p3,
+	if (!Dt.is_infinite(c->vertex(f.second)) &&
+        CGAL::side_of_bounded_sphere(p1, p2, p3,
 									 c->vertex(f.second)->point()) == CGAL::ON_BOUNDED_SIDE)
 		attached_ = true;
-	else if (CGAL::side_of_bounded_sphere(p1, p2, p3,
+	else if (!Dt.is_infinite(o->vertex(oi)) &&
+             CGAL::side_of_bounded_sphere(p1, p2, p3,
 										  o->vertex(oi)->point()) == CGAL::ON_BOUNDED_SIDE)
 		attached_ = true;
 	else
 		alpha_ = rho_;
 	
-	SimplexPhiMap::const_iterator c_iter = simplices.find(ARSimplex3D(*c,z));
-	SimplexPhiMap::const_iterator o_iter = simplices.find(ARSimplex3D(*o,z));
-	if (c_iter == simplices.end())			// c is infinite
+	if (Dt.is_infinite(c))
 	{
+	    SimplexPhiMap::const_iterator o_iter = simplices.find(ARSimplex3D(*o,z));
 		if (attached_) alpha_ = o_iter->first.alpha();
 		phi_const_ = o_iter->first.phi_const();					// FIXME: it's probably the other way around
 	}
-	else if (o_iter == simplices.end())		// o is infinite
+	else if (Dt.is_infinite(o))
 	{
+	    SimplexPhiMap::const_iterator c_iter = simplices.find(ARSimplex3D(*c,z));
 		if (attached_) alpha_ = c_iter->first.alpha();
 		phi_const_ = c_iter->first.phi_const();					// FIXME: it's probably the other way around
 	}
 	else
 	{
+	    SimplexPhiMap::const_iterator o_iter = simplices.find(ARSimplex3D(*o,z));
+	    SimplexPhiMap::const_iterator c_iter = simplices.find(ARSimplex3D(*c,z));
 		if (attached_) alpha_ = std::min(c_iter->first.alpha(), o_iter->first.alpha());
 		phi_const_ = std::min(c_iter->first.phi_const(), o_iter->first.phi_const());
 	}
@@ -208,8 +214,9 @@
 ARSimplex3D::
 operator<<(std::ostream& out) const
 {
+    out << this << ": ";
 	for (VertexSet::const_iterator cur = Parent::vertices().begin(); cur != Parent::vertices().end(); ++cur)
-		out << **cur << ", ";
+		out << &(**cur) << ", ";
 	out << "value = " << value();
 
 	return out;
@@ -226,19 +233,19 @@
 	ARSimplex3D::SimplexPhiMap simplices;
 	for(Cell_iterator cur = Dt.finite_cells_begin(); cur != Dt.finite_cells_end(); ++cur)
 		update_simplex_phi_map(ARSimplex3D(*cur, z), simplices);
-	std::cout << "Cells inserted" << std::endl;
+	rLog(rlARSimplex3D, "Cells inserted");
 	for(Vertex_iterator cur = Dt.finite_vertices_begin(); cur != Dt.finite_vertices_end(); ++cur)
 		simplices[ARSimplex3D(*cur, z)] = 0;			// only one tetrahedron can have non-negative phi_const value 
 														// (namely the one containing z); all other simplices will have a 
 														// negative phi_const value, so 0 is safe
-	std::cout << "Vertices inserted" << std::endl;
+	rLog(rlARSimplex3D, "Vertices inserted");
 
 	for(Facet_iterator cur = Dt.finite_facets_begin(); cur != Dt.finite_facets_end(); ++cur)
-		update_simplex_phi_map(ARSimplex3D(*cur, z, simplices), simplices);
-	std::cout << "Facets inserted" << std::endl;
+		update_simplex_phi_map(ARSimplex3D(*cur, z, simplices, Dt), simplices);
+	rLog(rlARSimplex3D, "Facets inserted");
 	for(Edge_iterator cur = Dt.finite_edges_begin(); cur != Dt.finite_edges_end(); ++cur)
-		update_simplex_phi_map(ARSimplex3D(*cur, z, simplices, Dt.incident_facets(*cur)), simplices);
-	std::cout << "Edges inserted" << std::endl;
+		update_simplex_phi_map(ARSimplex3D(*cur, z, simplices, Dt, Dt.incident_facets(*cur)), simplices);
+	rLog(rlARSimplex3D, "Edges inserted");
     
 	// Sort simplices by their alpha values
 	alpha_order.resize(simplices.size()); ARSimplex3DVector::iterator out = alpha_order.begin();
--- a/examples/ar-vineyard/ar-vineyard.cpp	Fri Apr 04 16:17:28 2008 -0400
+++ b/examples/ar-vineyard/ar-vineyard.cpp	Fri Apr 04 21:14:44 2008 -0400
@@ -4,37 +4,77 @@
 
 #include <iostream>
 #include <fstream>
+#include <string>
+#include <vector>
+
+#include <boost/program_options.hpp>
+namespace po = boost::program_options;
 
 
 int main(int argc, char** argv) 
 {
 #ifdef LOGGING
 	rlog::RLogInit(argc, argv);
+	stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+#endif
 
-	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") );
+	std::string infilename;
+	double zx,zy,zz,r;
+	std::string outfilename;
+
+    po::options_description hidden("Hidden options");
+    hidden.add_options()
+        ("input-file",  po::value<std::string>(&infilename), "Points file")
+        ("x",  po::value<double>(&zx), "x")
+        ("y",  po::value<double>(&zy), "y")
+        ("z",  po::value<double>(&zz), "z")
+        ("r",  po::value<double>(&r), "r")
+        ("output-file",  po::value<std::string>(&outfilename), "Vineyard edges output");
+
+    std::vector<std::string> log_channels;
+    po::options_description visible("Allowed options");
+    visible.add_options()
+        ("help,h",      "produce help message");
+#if LOGGING
+    visible.add_options()
+        ("log,l",       po::value< std::vector<std::string> >(&log_channels),
+                        "log channels to turn on");
+#endif
+    
+    po::positional_options_description p;
+    p.add("input-file", 1).add("x", 1).add("y", 1).add("z", 1).add("r", 1).add("output-file", 1);
+    
+    po::options_description all; all.add(visible).add(hidden);
+
+    po::variables_map vm;
+    po::store(po::command_line_parser(argc, argv).
+                  options(all).positional(p).run(), vm);
+    po::notify(vm);
+
+#if LOGGING
+    for (std::vector<std::string>::const_iterator cur = log_channels.begin(); cur != log_channels.end(); ++cur)
+        stdoutLog.subscribeTo( RLOG_CHANNEL(cur->c_str()) );
+    /* Interesting channels
+     * "ar/vineyard", "ar/function-kernel/value", "geometry/simulator",
+     * "topology/filtration", "topology/cycle", "topology/vineyard",
+     * "topology/filtration/transpositions", "topology/lowerstar"
+     */
 #endif
 
 	// Read command-line arguments
-	if (argc < 6)
+	if (vm.count("help") || !vm.count("input-file") || !vm.count("output-file") 
+                         || !vm.count("x") || !vm.count("y") || !vm.count("z") || !vm.count("r"))
 	{
-		std::cout << "Usage: ar-vineyard POINTS X Y Z OUTFILENAME" << std::endl;
+		std::cout << "Usage: ar-vineyard [OPTIONS] POINTS X Y Z R OUTFILENAME" << std::endl;
 		std::cout << "  POINTS       - filename containing points" << std::endl;
 		std::cout << "  X,Y,Z        - center-point z at which to compute the vineyard" << std::endl;
+		std::cout << "  R            - maximum radius" << std::endl;
 		std::cout << "  OUTFILENAME  - filename for the resulting vineyard" << std::endl;
+        std::cout << visible << std::endl;
 		std::cout << std::endl;
 		std::cout << "Computes an (alpha,r)-vineyard of the given pointset around the given point." << std::endl;
-		exit(0);
+		return 1;
 	}
-	std::string infilename = argv[1];
-	double zx,zy,zz; std::istringstream(argv[2]) >> zx;
-	std::istringstream(argv[3]) >> zy; std::istringstream(argv[4]) >> zz;
-	std::string outfilename = argv[5];
-	
 	
 	// Read in the point set and compute its Delaunay triangulation
 	std::ifstream in(infilename.c_str());
@@ -52,7 +92,7 @@
 	arv.compute_pairing();
 
 	// Compute vineyard
-	arv.compute_vineyard(true);
+	arv.compute_vineyard(r);
 	std::cout << "Vineyard computed" << std::endl;
 	arv.vineyard()->save_edges(outfilename);
 }
--- a/examples/ar-vineyard/ar-vineyard.h	Fri Apr 04 16:17:28 2008 -0400
+++ b/examples/ar-vineyard/ar-vineyard.h	Fri Apr 04 21:14:44 2008 -0400
@@ -6,70 +6,109 @@
 #ifndef __AR_VINEYARD_H__
 #define __AR_VINEYARD_H__
 
+#include <boost/signals.hpp>
+#include <boost/bind.hpp>
+#include <list>
+#include <vector>
+
 #include "topology/conesimplex.h"
 #include "topology/filtration.h"
+#include "geometry/kinetic-sort.h"
+#include "geometry/simulator.h"
 
-#include <CGAL/Kinetic/Inexact_simulation_traits.h>
-#include <CGAL/Kinetic/Event_base.h>
-#include <CGAL/Kinetic/Sort.h>
-#include <CGAL/Kinetic/Sort_visitor_base.h>
-
-#include <list>
 #include "ar-simplex3d.h"
-
-#include <vector>
+#include "ar-function-kernel.h"
 
 
-class ARVineyardBase
-{
-	public:
-		/// \name CGAL Kinetic Sort types
-		/// @{
-		class						SortVisitor;
-		typedef 					CGAL::Kinetic::Inexact_simulation_traits 					Traits;
-		typedef						CGAL::Kinetic::Sort<Traits, SortVisitor>					Sort;
-		typedef 					Traits::Simulator 											Simulator;
-		typedef 					Traits::Active_points_1_table								ActivePointsTable;
-		typedef 					ActivePointsTable::Key										Key;
-		
-		typedef 					Traits::Kinetic_kernel::
-											Function_kernel::Construct_function 				CF; 
-		typedef 					Traits::Kinetic_kernel::Motion_function 					F; 
-		/// @}
-		
-		class						ARConeSimplex;
-		class						MembershipFunctionChangeEvent;
-};
-
-class ARVineyardBase::ARConeSimplex: public ConeSimplex<ARSimplex3D>
+template <class Simulator_>
+class ARConeSimplex3D: public ConeSimplex<ARSimplex3D>
 {
 	public:
 		typedef						ConeSimplex<ARSimplex3D>									Parent;
 		typedef						ARSimplex3D													ARSimplex3D;
+		
+		/// \name Simulator types
+		/// @{
+        typedef                     Simulator_                                                  Simulator;
+        typedef                     typename Simulator::FunctionKernel                          FunctionKernel;
+        typedef                     typename FunctionKernel::Function                           Function;
+        /// @}
+		
+		/// \name ThresholdSort types
+		/// @{
+		typedef 					std::list<Function>										    ThresholdList;
+        typedef                     typename ThresholdList::iterator                            ThresholdListIterator;
 
-									ARConeSimplex(const ARSimplex3D& s, bool coned = false): 
-										Parent(s, coned)										{}
+		struct 						ThresholdTrajectoryExtractor
+		{	Function                operator()(ThresholdListIterator i) const		            { return *i; } };
+
+		typedef						KineticSort<ThresholdListIterator, 
+                                                ThresholdTrajectoryExtractor, Simulator>		ThresholdSort;
+		/// @}
 
-		Key							kinetic_key() const											{ return key_; }
-		void						set_kinetic_key(Key k)										{ key_ = k; }
+        typedef                     boost::signal<void (Simulator*)>                            NewMaxSignal;
+    
+    public:
+									ARConeSimplex3D(const ARSimplex3D& s, bool coned = false);
+									ARConeSimplex3D(const Parent& p): Parent(p)                 {}      // crucial for boundary() to work correctly
+									ARConeSimplex3D(const ARConeSimplex3D& other):              // need explicit copy-constructor because of the signal
+                                        Parent(other, other.coned()), 
+                                        thresholds_(other.thresholds_)                          {}
+
+		const ThresholdList&        thresholds() const											{ return thresholds_; }
+
+        NewMaxSignal&               new_max_signal()                                            { return new_max_signal_; }
+        const Function&             max_threshold() const                                       { return thresholds_.back(); }
+		void						schedule_thresholds(Simulator* simulator);
+
+        // need explicit operator= because of the signal
+        ARConeSimplex3D&            operator=(const ARConeSimplex3D& other)                     { Parent::operator=(other); thresholds_ = other.thresholds_; return *this; }
+
 								
 	private:
-		Key							key_;
+		ThresholdList				thresholds_;
+		ThresholdSort				thresholds_sort_;
+        NewMaxSignal                new_max_signal_;
+
+		void						swap_thresholds(ThresholdListIterator i, Simulator* simulator);
 };
 
-
-class ARVineyard: public ARVineyardBase
+/**
+ * Encapsulated filtration, and provides compute_vineyard() functionality.
+ */
+class ARVineyard
 {
 	public:
 		typedef						ARVineyard													Self;
+	
+        /// \name FunctionKernel and Simulator types
+        /// @{
+        typedef                     ARFunctionKernel                                            FunctionKernel;
+        typedef                     FunctionKernel::Function                                    Function;
+        typedef                     Simulator<FunctionKernel>                                   SimulatorAR;
+        /// @}
+
+        /// \name Filtration types
+        /// @{    
+        typedef                     ARConeSimplex3D<SimulatorAR>                                ConeSimplex3D;
+		typedef						Filtration<ConeSimplex3D>									FiltrationAR;
+		typedef						FiltrationAR::Simplex										Simplex;
+		typedef						FiltrationAR::Index											Index;
+		typedef						FiltrationAR::Vineyard										Vineyard;
+		typedef						Vineyard::Evaluator											Evaluator;
+        /// @}
 		
-		typedef						Filtration<ARConeSimplex>									ARFiltration;	
-		typedef						ARFiltration::Simplex										Simplex;
-		typedef						ARFiltration::Index											Index;
-		typedef						ARFiltration::Vineyard										Vineyard;
-		typedef						Vineyard::Evaluator											Evaluator;
-		typedef						std::map<Key, Index>										KeyIndexMap;
+        /// \name SimplexSort types
+        /// @{
+        struct 						SimplexTrajectoryExtractor
+		{	Function				operator()(Index i) const									{ return i->max_threshold(); } };
+
+		typedef						KineticSort<Index, SimplexTrajectoryExtractor, SimulatorAR> SimplexSort;
+		typedef						SimplexSort::iterator										SimplexSortIterator;
 		
+        class                       ThresholdChangeSlot;              // used to notify of change in max threshold
+		/// @}
+
 		typedef						std::list<Point>											PointList;
 
 		class						StaticEvaluator;
@@ -80,26 +119,23 @@
 									~ARVineyard();
 
 		void						compute_pairing();
-		void						compute_vineyard(bool explicit_events = false);
+		void						compute_vineyard(double max_radius);
 		
-		const ARFiltration*			filtration() const											{ return filtration_; }
+		const FiltrationAR*			filtration() const											{ return filtration_; }
 		const Vineyard*				vineyard() const											{ return vineyard_; }
 
 	public:
-		// For Kinetic Sort
-		void 						swap(Key a, Key b);
+		void 						swap(Index i, SimulatorAR* simulator);						///< For kinetic sort
 	
 	private:
 		void 						add_simplices();
 		void						change_evaluator(Evaluator* eval);
 
 	private:
-		ARFiltration*				filtration_;
+		FiltrationAR*				filtration_;
 		Vineyard*					vineyard_;
 		Evaluator*					evaluator_;
 
-		KeyIndexMap					kinetic_map_;
-
 		Point						z_;
 		Delaunay					dt_;
 				
@@ -120,68 +156,58 @@
 
 //BOOST_CLASS_EXPORT(ARVineyard)
 
+#ifdef COUNTERS
+static Counter*  cARVineyardTrajectoryKnee =		 GetCounter("ar/vineyard/trajectoryknee");
+#endif
 
-class ARVineyardBase::MembershipFunctionChangeEvent: public CGAL::Kinetic::Event_base<int*>
-{
-	public:
-									MembershipFunctionChangeEvent(Key k, F function, 
-																  ActivePointsTable::Handle apt):
-										key_(k), function_(function), apt_(apt)					{}
-		
-		void						process() const;
-		std::ostream&				operator<<(std::ostream& out) const;
-		std::ostream&				write(std::ostream& out) const								{ return this->operator<<(out); }
-
-	private:
-		Key							key_;
-		F							function_;
-		ActivePointsTable::Handle	apt_;
+class ARVineyard::ThresholdChangeSlot
+{   
+    public:
+                                ThresholdChangeSlot(SimplexSortIterator iter, SimplexSort* sort, 
+                                                    Vineyard* vineyard, SimulatorAR* sort_simulator):
+                                    iter_(iter), sort_(sort), vineyard_(vineyard), sort_simulator_(sort_simulator)      { iter_->element->new_max_signal().connect(*this); }
+        void                    operator()(SimulatorAR* simulator)                                
+        { 
+            Count(cARVineyardTrajectoryKnee); 
+            sort_->update_trajectory(iter_, sort_simulator_); 
+            if (iter_->element->sign()) 
+                vineyard_->record_knee(iter_->element);
+            else
+                vineyard_->record_knee(iter_->element->pair());
+        }
+    
+    private:
+        SimplexSortIterator     iter_;
+        SimplexSort*            sort_;              // could make both of these static
+        Vineyard*               vineyard_;          // currently inefficient since there is
+                                                    // only one SimplexSort and one Vineyard, 
+                                                    // but each is stored in every slot
+        SimulatorAR*              sort_simulator_;
 };
 
-std::ostream& operator<<(std::ostream& out, const ARVineyardBase::MembershipFunctionChangeEvent& e)
-{ return e.operator<<(out); }
-
 class ARVineyard::StaticEvaluator: public Evaluator
 {
 	public:
-									StaticEvaluator(RealType t): time_(t)						{}
+									StaticEvaluator()                   						{}
 
-		virtual RealType			time() const												{ return time_; }
+		virtual RealType			time() const												{ return 0; }
 		virtual RealType			value(const Simplex& s) const								{ return s.value(); }
-
-	private:
-		RealType					time_;
 };
 
 class ARVineyard::KineticEvaluator: public Evaluator
 {
 	public:
-									KineticEvaluator(Simulator::Handle sp, 
-													 ActivePointsTable::Handle apt, 
-													 RealType time_offset): 
-										sp_(sp), apt_(apt)										{}
+									KineticEvaluator(SimulatorAR* simplex_sort_simulator,
+                                                     SimulatorAR* trajectory_sort_simulator): 
+                                        simplex_sort_simulator_(simplex_sort_simulator),
+                                        trajectory_sort_simulator_(trajectory_sort_simulator)   {}
 
-		virtual RealType			time() const												{ return CGAL::to_double(get_time()); }
-		virtual RealType			value(const Simplex& s)	const								{ return CGAL::to_double(apt_->at(s.kinetic_key()).x()(get_time())); }
+		virtual RealType			time() const												{ return std::max(simplex_sort_simulator_->current_time(), trajectory_sort_simulator_->current_time()); }
+		virtual RealType			value(const Simplex& s)	const								{ return FunctionKernel::value_at(s.max_threshold(), time()); }
 
 	private:
-		Simulator::Time				get_time() const											{ return sp_->current_time(); }
-		
-		Simulator::Handle			sp_;
-		ActivePointsTable::Handle 	apt_;
-};
-
-
-class ARVineyardBase::SortVisitor: public CGAL::Kinetic::Sort_visitor_base
-{
-	public:
-									SortVisitor(ARVineyard* arv): arv_(arv)						{}
-
-		template<class Vertex_handle>
-		void						before_swap(Vertex_handle a, Vertex_handle b) const;
-
-	private:
-		ARVineyard*					arv_;
+		SimulatorAR*                simplex_sort_simulator_;
+        SimulatorAR*                trajectory_sort_simulator_;
 };
 
 
--- a/examples/ar-vineyard/ar-vineyard.hpp	Fri Apr 04 16:17:28 2008 -0400
+++ b/examples/ar-vineyard/ar-vineyard.hpp	Fri Apr 04 21:14:44 2008 -0400
@@ -1,26 +1,68 @@
 #include <utilities/log.h>
 
 /* Implementation */
+
+#ifdef LOGGING
+static rlog::RLogChannel* rlARVineyard =                        DEF_CHANNEL("ar/vineyard", rlog::Log_Debug);
+static rlog::RLogChannel* rlARVineyardComputing =               DEF_CHANNEL("ar/vineyard/computing", rlog::Log_Debug);
+static rlog::RLogChannel* rlARVineyardSwap =                    DEF_CHANNEL("ar/vineyard/swap", rlog::Log_Debug);
+static rlog::RLogChannel* rlARVineyardThresholdSwap =           DEF_CHANNEL("ar/vineyard/threshold/swap", rlog::Log_Debug);
+#endif
+
+
+template <class Simulator_>
+ARConeSimplex3D<Simulator_>::
+ARConeSimplex3D(const ARSimplex3D& s, bool coned): Parent(s, coned)
+{}
 	
+template <class Simulator_>
+void
+ARConeSimplex3D<Simulator_>::
+swap_thresholds(ThresholdListIterator i, Simulator* simulator)
+{
+    rLog(rlARVineyardThresholdSwap, "Transposing %s and %s", tostring(*i).c_str(),
+                                                             tostring(*boost::next(i)).c_str());
+	typename ThresholdList::iterator n = boost::next(i);
+	thresholds_.splice(i, thresholds_, n);
+	if (boost::next(i) == thresholds_.end())
+        new_max_signal_(simulator);
+}
+
+template <class Simulator_>
+void
+ARConeSimplex3D<Simulator_>::
+schedule_thresholds(Simulator* simulator)
+{
+    if (!coned()) thresholds_.push_back(Function(Function::rho, this));
+    else        
+    { 
+        thresholds_.push_back(Function(Function::phi, this)); 
+        thresholds_.push_back(Function(Function::rho, this)); 
+	    thresholds_sort_.initialize(thresholds_.begin(), thresholds_.end(), 
+		    						boost::bind(&ARConeSimplex3D::swap_thresholds, this, _1, _2), simulator);
+    }
+}
+
+
 ARVineyard::
 ARVineyard(const PointList& points, const Point& z): z_(z)
 {
 	for (PointList::const_iterator cur = points.begin(); cur != points.end(); ++cur)
 		dt_.insert(*cur);
-	std::cout << "Delaunay triangulation computed" << std::endl;
+	rLog(rlARVineyard, "Delaunay triangulation computed");
 
 	ARSimplex3DVector alpha_ordering;
 	fill_alpha_order(dt_, z_, alpha_ordering);
-	std::cout << "Delaunay simplices: " << alpha_ordering.size() << std::endl;
+	rLog(rlARVineyard, "Delaunay simplices: %i", alpha_ordering.size());
 		
-	evaluator_ = new StaticEvaluator(0);
+	evaluator_ = new StaticEvaluator;
 	vineyard_ = new Vineyard(evaluator_);
 
-	filtration_ = new ARFiltration(vineyard_);
+	filtration_ = new FiltrationAR(vineyard_);
 	for (ARSimplex3DVector::const_iterator cur = alpha_ordering.begin(); cur != alpha_ordering.end(); ++cur)
 	{
-		filtration_->append(*cur);						// Delaunay simplex
-		filtration_->append(ARConeSimplex(*cur, true));	// Coned off delaunay simplex
+		filtration_->append(ConeSimplex3D(*cur));         // Delaunay simplex
+		filtration_->append(ConeSimplex3D(*cur, true));   // Coned off delaunay simplex
 	}
 }
 
@@ -39,98 +81,64 @@
 	filtration_->fill_simplex_index_map();
 	filtration_->pair_simplices(filtration_->begin(), filtration_->end());
 	vineyard_->start_vines(filtration_->begin(), filtration_->end());
-	std::cout << "Simplices paired" << std::endl;
+	rLog(rlARVineyard, "Simplices paired");
 }
 
 void					
 ARVineyard::
-compute_vineyard(bool explicit_events)
+compute_vineyard(double max_radius)
 {
 	AssertMsg(filtration_->is_paired(), "Simplices must be paired for a vineyard to be computed");
 	
-	typedef Traits::Kinetic_kernel::Point_1 								Point_1;
-	typedef Simulator::Time													Time;
+	SimulatorAR simplex_sort_simulator, trajectory_sort_simulator;
+	SimplexSort	simplex_sort;
 	
-	Traits tr(0,1);
-	Simulator::Handle sp = tr.simulator_handle();
-	ActivePointsTable::Handle apt = tr.active_points_1_table_handle();
-	Sort sort(tr, SortVisitor(this));
-	
-	// Setup the kinetic sort and membership changes
-	std::cout << "Setting up the kinetic sort and membership events" << std::endl;
-	CF cf; 
-	kinetic_map_.clear();
+	// Schedule thresholds
 	for (Index cur = filtration_->begin(); cur != filtration_->end(); ++cur)
-	{
-		F x = cf(F::NT(CGAL::to_double(cur->alpha())));
-		Point_1 p(x);
-		cur->set_kinetic_key(apt->insert(p));
-		kinetic_map_[cur->kinetic_key()] = cur;
-			
-		if (!cur->coned()) continue;						// non-coned simplices stay put, so we are done
+        cur->schedule_thresholds(&trajectory_sort_simulator);
 
-		Time lambda_alpha = CGAL::to_double((cur->alpha() - cur->rho()));	// when lambda becomes greater than alpha
-		lambda_alpha += 2*CGAL::sqrt(CGAL::to_double(cur->s()*lambda_alpha));
-		lambda_alpha += CGAL::to_double(cur->s() + cur->v());
-
-		Time phi_alpha = CGAL::to_double(cur->alpha() - cur->phi_const());
-
-		Time phi_lambda = CGAL::to_double(cur->rho() + cur->s() - cur->v() - cur->phi_const());
-		phi_lambda *= phi_lambda;
-		phi_lambda /= CGAL::to_double(4*cur->s());
-		phi_lambda += CGAL::to_double(cur->v());
+	// Once thresholds are scheduled, we can initialize the simplex_sort
+	simplex_sort.initialize(filtration_->begin(), filtration_->end(), 
+							boost::bind(&ARVineyard::swap, this, _1, _2), &simplex_sort_simulator);
+    rLog(rlARVineyardComputing, "SimplexSort initialized");
 
-		Time sv = CGAL::to_double(cur->s() + cur->v());
-		
-		if (true || phi_lambda < sv || phi_lambda < phi_alpha)		// FIXME: remove true
-		{
-			sp->new_event(Time(phi_alpha), 
-						  MembershipFunctionChangeEvent(cur->kinetic_key(),
-														cf(F::NT(CGAL::to_double(cur->phi_const())), 1),
-														apt));		// \phi^2 = r^2 + \phi_c^2
-			std::cout << "Scheduled" << std::endl;
-		} else
-			std::cout << "Not scheduled" << std::endl;
-		
-
-		//sp->new_event(Time(...), MembershipFunctionChangeEvent(cur->kinetic_key()));
-		
-		std::cout << *cur << std::endl;
-		std::cout << "lambda_alpha: " 		<< lambda_alpha << std::endl;
-		std::cout << "phi_alpha: " 			<< phi_alpha << std::endl;
-		std::cout << "phi_lambda: " 		<< phi_lambda << std::endl;
-		std::cout << "s^2 + v^2: " 			<< sv << std::endl;
-		std::cout << std::endl;
-	}
+    // Connect signals and slots
+    std::vector<ThresholdChangeSlot> slots; 
+    slots.reserve(filtration_->size());
+    for (SimplexSortIterator cur = simplex_sort.begin(); cur != simplex_sort.end(); ++cur)
+        slots.push_back(ThresholdChangeSlot(cur, &simplex_sort, vineyard_, &simplex_sort_simulator));
+    rLog(rlARVineyardComputing, "Signals and slots connected");
+    rLog(rlARVineyardComputing, "SimplexSort size: %i", simplex_sort_simulator.size());
+    rLog(rlARVineyardComputing, "TrajectorySort size: %i", trajectory_sort_simulator.size());
 	
-	// Process all the events (compute the vineyard in the process)
-	// FIXME: the time should not be 1, but something like twice the radius of
-	// the pointset as seen from z
-	change_evaluator(new KineticEvaluator(sp, apt, 0));
-	if (explicit_events)
-	{
-		while (sp->next_event_time() < 1)
-		{
-			std::cout << "Next event time: " << sp->next_event_time() << std::endl;
-			sp->set_current_event_number(sp->current_event_number() + 1);
-			std::cout << "Processed event" << std::endl;
-		}
-	} else
-		sp->set_current_time(1.0);
-	std::cout << "Processed " << sp->current_event_number() << " events" << std::endl;
+    // Simulate
+	change_evaluator(new KineticEvaluator(&simplex_sort_simulator, &trajectory_sort_simulator));
+    while ((simplex_sort_simulator.current_time() < max_radius || trajectory_sort_simulator.current_time() < max_radius) && !(simplex_sort_simulator.reached_infinity() && trajectory_sort_simulator.reached_infinity()))
+    {
+        if (*(simplex_sort_simulator.top()) < *(trajectory_sort_simulator.top()))
+        {
+            rLog(rlARVineyardComputing, "Current time before: %lf (simplex sort)", simplex_sort_simulator.current_time());
+            rLog(rlARVineyardComputing, "Top event: %s (simplex sort)", intostring(*(simplex_sort_simulator.top())).c_str());
+            simplex_sort_simulator.process();
+        } else
+        {
+            rLog(rlARVineyardComputing, "Current time before: %lf (trajectory sort)", trajectory_sort_simulator.current_time());
+            rLog(rlARVineyardComputing, "Top event: %s (trajectory sort)", intostring(*(trajectory_sort_simulator.top())).c_str());
+            trajectory_sort_simulator.process();
+        }
+    }
 	
-	//change_evaluator(new StaticEvaluator(1));
 	vineyard_->record_diagram(filtration_->begin(), filtration_->end());
 }
 		
 void 					
 ARVineyard::
-swap(Key a, Key b)
+swap(Index i, SimulatorAR* simulator)
 {
-	Index ao = kinetic_map_[a], bo = kinetic_map_[b];
-	AssertMsg(filtration_->get_trails_cmp()(ao, bo), "In swap(a,b), a must precede b");
-	filtration_->transpose(ao);
-	AssertMsg(filtration_->get_trails_cmp()(bo, ao), "In swap(a,b), b must precede a after the transposition");
+    rLog(rlARVineyardSwap, "Transposing %p and %p:", &(*i), &(*boost::next(i)));
+    rLog(rlARVineyardSwap, "  %s and", tostring(*i).c_str());
+    rLog(rlARVineyardSwap, "  %s", tostring(*boost::next(i)).c_str());
+	filtration_->transpose(i);
 }
 
 void
@@ -143,30 +151,3 @@
 	evaluator_ = eval;
 	vineyard_->set_evaluator(evaluator_);
 }
-
-void 
-ARVineyardBase::MembershipFunctionChangeEvent::
-process() const
-{
-	apt_->set(key_, function_);
-	std::cout << "Updated for phi's dominance" << std::endl;
-}
-
-
-template<class Vertex_handle>
-void
-ARVineyardBase::SortVisitor::
-before_swap(Vertex_handle a, Vertex_handle b) const
-{ 
-	std::cout << "Swapping elements" << *a << " and " << *b << std::endl;
-	arv_->swap(*a,*b); 
-}
-
-
-std::ostream&
-ARVineyardBase::MembershipFunctionChangeEvent::
-operator<<(std::ostream& out) const
-{
-	return out << "Membership change" << std::endl;
-}
-
--- a/include/dionysus.h	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/dionysus.h	Fri Apr 04 21:14:44 2008 -0400
@@ -1,6 +1,6 @@
 /*
  * Author: Dmitriy Morozov
- * Department of Computer Science, Duke University, 2006
+ * Department of Computer Science, Duke University, 2006--2008
  *
  * For now file exists only to store the main page of the documentation
  */
@@ -9,3 +9,16 @@
  * Detailed description of Dionysus with references to all the classes, 
  * and code samples goes here.
  */
+
+/**
+ * \defgroup topology Topology Classes
+ */
+
+/**
+ * \defgroup geometry Geometry Classes
+ */
+
+/**
+ * \defgroup kinetic Kinetic Data Structures Classes
+ * \ingroup  geometry
+ */
--- a/include/geometry/euclidean.h	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/geometry/euclidean.h	Fri Apr 04 21:14:44 2008 -0400
@@ -14,15 +14,19 @@
 #include "number-traits.h"
 
 
+/**
+ * Geometric Kernel. Defines operations on geometric primitives.
+ * \ingroup geometry
+ */
 template<class NumberType_ = double>
 class Kernel
 {
 	public:
 		typedef 					unsigned int								DimensionType;
 		typedef						NumberType_									NumberType;
-		typedef						LinearAlgebra<NumberType>					LinearAlgebra;
-		typedef						typename LinearAlgebra::MatrixType			MatrixType;
-		typedef						typename LinearAlgebra::VectorType			VectorType;
+		typedef						LinearAlgebra<NumberType>					LinearAlgebraK;
+		typedef						typename LinearAlgebraK::MatrixType			MatrixType;
+		typedef						typename LinearAlgebraK::VectorType			VectorType;
 
 		class						Point;
 		class						Sphere;
@@ -56,7 +60,10 @@
 };
 
 
-/* Point */
+/** 
+ * Point class.
+ * \ingroup geometry 
+ */
 template<class NumberType_>
 class Kernel<NumberType_>::Point: public VectorType
 {
@@ -77,7 +84,10 @@
 };
 
 
-/* Sphere */
+/** 
+ * Sphere class.
+ * \ingroup geometry
+ */
 template<class NumberType_>
 class Kernel<NumberType_>::Sphere
 {
--- a/include/geometry/kinetic-sort.h	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/geometry/kinetic-sort.h	Fri Apr 04 21:14:44 2008 -0400
@@ -10,33 +10,37 @@
  * Maintains elements of the given data structure in the sorted order assuming the elements follow 
  * trajectories given by TrajectoryExtractor_.
  *
- *  \arg SortDS_              should be forward and backward iterable, swaps are handles via SwapCallback
- *  \arg TrajectoryExtractor_ applied to the iterator into SortDS_ should return a rational 
- *                            function describing the 
+ *  \arg ElementIterator_     iterator over the underlying data structure that's kept in sorted order
+ *  \arg TrajectoryExtractor_ applied to the iterator into SortDS_ should return a function 
+ *                            (of type Simulator_::FunctionKernel::Function) describing the trajectory of the element
  *  \arg Simulator_           the Simulator type, e.g. Simulator. Note that KineticSort does not store 
  *                            a pointer to the Simulator (so a pointer is passed in each relevant operation)
+ *  \arg Swap_                is called with an ElementIterator_ when a swap needs to be performed
+ *
+ *  \ingroup kinetic
  */
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+template<class ElementIterator_, class TrajectoryExtractor_, 
+		 class Simulator_, class Swap_ = boost::function<void(ElementIterator_ pos, Simulator_* simulator)> >
 class KineticSort
 {
 	public:
 		typedef						Simulator_									Simulator;
-		typedef						typename Simulator::PolynomialKernel		PolynomialKernel;
-		typedef						SortDS_										SortDS;
+		typedef						typename Simulator::FunctionKernel		    FunctionKernel;
+		typedef						ElementIterator_							ElementIterator;
+		typedef						Swap_										Swap;
 		typedef						TrajectoryExtractor_						TrajectoryExtractor;
 		
 		typedef						typename Simulator::Key						SimulatorKey;
-		typedef						typename SortDS::iterator					SortDSIterator;
 
 									
 	private:
 		/* Implementation */
 		struct Node
 		{
-			SortDSIterator			element;
+			ElementIterator			element;
 			SimulatorKey			swap_event_key;		
 
-									Node(SortDSIterator e, SimulatorKey k): 
+									Node(ElementIterator e, SimulatorKey k): 
 										element(e), swap_event_key(k)			{}
 		};
 
@@ -44,16 +48,15 @@
 
 	public:
 		typedef						typename NodeList::iterator					iterator;
-		typedef						boost::function<void(SortDS*, SortDSIterator pos)>	
-																				SwapCallback;
 
 
 		/// \name Core Functionality
 		/// @{
-									KineticSort(SortDS* sort, Simulator* simulator, SwapCallback swap_callback);
+									KineticSort();
+									KineticSort(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator);
+		void						initialize(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator);
 
-		template<class InputIterator>
-		void						insert(iterator pos, InputIterator f, InputIterator l, Simulator* simulator);
+		void						insert(iterator pos, ElementIterator f, ElementIterator l, Simulator* simulator);
 		void						erase(iterator pos, Simulator* simulator);
 		void						update_trajectory(iterator pos, Simulator* simulator);
 
@@ -62,6 +65,9 @@
 		bool						audit(Simulator* simulator) const;
 		/// @}
 
+		iterator					begin() 									{ return list_.begin(); }
+		iterator					end() 										{ return list_.end(); }
+
 	private:
 		class SwapEvent;
 		void						schedule_swaps(iterator b, iterator e, Simulator* s);
@@ -69,8 +75,7 @@
 
 	private:
 		NodeList					list_;
-		SortDS*						sort_;			
-		SwapCallback				swap_callback_;	
+		Swap						swap_;	
 };
 
 #include "kinetic-sort.hpp"
--- a/include/geometry/kinetic-sort.hpp	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/geometry/kinetic-sort.hpp	Fri Apr 04 21:14:44 2008 -0400
@@ -1,49 +1,78 @@
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
-KineticSort(SortDS* sort, Simulator* simulator, SwapCallback swap_callback):
-	sort_(sort), swap_callback_(swap_callback)
+#include "utilities/log.h"
+#include "utilities/counter.h"
+
+#ifdef LOGGING
+static rlog::RLogChannel* rlKineticSort =           DEF_CHANNEL("geometry/kinetic-sort", rlog::Log_Debug);
+static rlog::RLogChannel* rlKineticSortAudit =      DEF_CHANNEL("geometry/kinetic-sort/audit", rlog::Log_Debug);
+static rlog::RLogChannel* rlKineticSortSchedule =   DEF_CHANNEL("geometry/kinetic-sort/schedule", rlog::Log_Debug);
+static rlog::RLogChannel* rlKineticSortProcess =    DEF_CHANNEL("geometry/kinetic-sort/process", rlog::Log_Debug);
+#endif // LOGGING
+
+#ifdef COUNTERS
+static Counter*  cKineticSort =                     GetCounter("kinetic-sort");
+static Counter*  cKineticSortSwap =                 GetCounter("kinetic-sort/swap");
+#endif // COUNTERS
+
+
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
+KineticSort()
+{}	
+
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
+KineticSort(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator)
 {
-	for (SortDSIterator cur = sort->begin(); cur != sort->end(); ++cur)
+	initialize(b, e, swap, simulator);
+}
+
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
+void
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
+initialize(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator)
+{
+	swap_ = swap;
+	for (ElementIterator cur = b; cur != e; ++cur)
 		list_.push_back(Node(cur, simulator->null_key()));
 	schedule_swaps(list_.begin(), list_.end(), simulator);
 }
 
-
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
-template<class InputIterator>
+/// Adds elements in the range [f,l) of the underlying data structure to the management by the KineticSort.
+/// pos must be a valid iterator, i.e., it cannot be end().
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
 void
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
-insert(iterator pos, InputIterator f, InputIterator l, Simulator* simulator)
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
+insert(iterator pos, ElementIterator f, ElementIterator l, Simulator* simulator)
 {
 	iterator previous = pos; --previous;
 	if (previous != list_.end()) simulator->remove(previous->swap_event_key);
 
-	sort_->insert(pos->element, f, l);
+	ElementIterator cur = boost::next(previous)->element;
+	while(cur != pos->element)
+		list_.insert(pos->element, Node(cur++, simulator->null_key()));
 
-	SortDSIterator cur = boost::next(previous)->element;
-	while(cur != pos->element)
-		list_.insert(pos->element, Node(cur++));
 	if (previous != list_.end()) 
 		schedule_swaps(previous, pos, simulator);
 	else
 		schedule_swaps(list_.begin(), pos, simulator);
 }
 
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+/// Removes pos from the KineticSort structure. Assumes that the element of the underlying data structure 
+/// that *pos refers to has already been removed.
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
 void
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
 erase(iterator pos, Simulator* simulator)
 {
 	simulator->remove(pos->swap_event_key);
-	sort_->erase(pos->element);
 	iterator prev = pos; --prev;
 	list_.erase(pos);
 	schedule_swaps(prev, simulator);
 }
 
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
 void
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
 update_trajectory(iterator pos, Simulator* simulator)
 {
 	iterator prev = boost::prior(pos);
@@ -61,15 +90,16 @@
 }
 
 
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
 void						
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
 swap(iterator pos, Simulator* simulator)
 {
-	swap_callback_(sort_, pos->element);
-	
-	// TODO: add assertion that boost::next(pos) != list_.end()
-	
+	// TODO: AssertMsg(boost::next(pos) != list_.end(), "Cannot swap the last element");
+
+    Count(cKineticSortSwap);
+	swap_(pos->element, simulator);
+
 	// Remove events
 	iterator prev = boost::prior(pos);
 	if (prev != list_.end())
@@ -88,60 +118,63 @@
 	//audit(simulator);
 }
 
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
 bool
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
 audit(Simulator* simulator) const
 {
-	typedef 		typename Simulator::RationalFunction		RationalFunction;
+	typedef 		typename Simulator::Function		        Function;
 	typedef 		typename Simulator::Time					Time;
 	
 	Time t = simulator->audit_time();
-	std::cout << "Auditing at " << t << std::endl;
+	rLog(rlKineticSortAudit, "Auditing at %s", tostring(t).c_str());
 
 	TrajectoryExtractor	te;
 	
 	typename NodeList::const_iterator next = list_.begin();
 	typename NodeList::const_iterator cur = next++;
-	RationalFunction cur_trajectory = te(cur->element);
+	Function cur_trajectory = te(cur->element);
 	while (next != list_.end())
 	{
-		(*(cur->swap_event_key))->print(std::cout << "  ") << std::endl;
+		rLog(rlKineticSortAudit, "  %s", intostring(**(cur->swap_event_key)).c_str());
 
-		RationalFunction next_trajectory = te(next->element);
-		std::cout << "  Auditing:   " << cur_trajectory << ", " << next_trajectory << std::endl;
-		std::cout << "  Difference: " << next_trajectory - cur_trajectory << std::endl;
-		std::cout << "  Sign at:    " << t << ", " << PolynomialKernel::sign_at(next_trajectory - cur_trajectory, t) << std::endl;
-		if (PolynomialKernel::sign_at(next_trajectory - cur_trajectory, t) == -1)
+		Function next_trajectory = te(next->element);
+		rLog(rlKineticSortAudit, "  Auditing:   %s, %s", tostring(cur_trajectory).c_str(),
+                                                         tostring(next_trajectory).c_str());
+		rLog(rlKineticSortAudit, "  Difference: %s", tostring(next_trajectory - cur_trajectory).c_str());
+		rLog(rlKineticSortAudit, "  Sign at:    %s, %s", tostring(t).c_str(),
+                                                         tostring(FunctionKernel::sign_at(next_trajectory - cur_trajectory, t)).c_str());
+		if (FunctionKernel::sign_at(next_trajectory - cur_trajectory, t) == -1)
 		{
-			std::cout << "Audit failed at " << *cur->element << ", " << *next->element << std::endl;
+			rError("Audit failed at %s, %s", tostring(*cur->element).c_str(), 
+                                             tostring(*next->element).c_str());
 			return false;
 		}
 
 		cur_trajectory = next_trajectory;
 		cur = next++;
 	}
-	if (cur != list_.end()) (*(cur->swap_event_key))->print(std::cout << "  ") << std::endl;
+	if (cur != list_.end()) rLog(rlKineticSortAudit, "  %s", intostring(**(cur->swap_event_key)).c_str());
 	return true;
 }
 
 		
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
 void						
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
 schedule_swaps(iterator b, iterator e, Simulator* simulator)
 {
-	typedef 		typename Simulator::RationalFunction		RationalFunction;
+	typedef 		typename Simulator::Function		        Function;
 	
 	TrajectoryExtractor	te;
 	
 	iterator next = b; 
 	iterator cur = next++;
-	RationalFunction cur_trajectory = te(cur->element);
+	Function cur_trajectory = te(cur->element);
 	while (next != e)
 	{
-		RationalFunction next_trajectory = te(next->element);
-		std::cout << "Next trajectory: " << next_trajectory << std::endl;
+		Function next_trajectory = te(next->element);
+		rLog(rlKineticSortSchedule, "Next trajectory: %s", tostring(next_trajectory).c_str());
 		// TODO: add assertion that (next_trajectory - cur_trajectory)(s->curren_time()) > 0
 		cur->swap_event_key = simulator->add(next_trajectory - cur_trajectory, SwapEvent(this, cur));
 		cur = next++;
@@ -150,12 +183,12 @@
 	if (cur != e) schedule_swaps(cur, simulator);
 }
 
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
 void						
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
 schedule_swaps(iterator i, Simulator* simulator)
 {
-	typedef 		typename Simulator::RationalFunction		RationalFunction;
+	typedef 		typename Simulator::Function		        Function;
 	
 	if (i == list_.end()) return;
 	if (boost::next(i) == list_.end())
@@ -167,8 +200,8 @@
 	TrajectoryExtractor	te;
 	
 	iterator next = boost::next(i); 
-	RationalFunction i_trajectory = te(i->element);
-	RationalFunction next_trajectory = te(next->element);
+	Function i_trajectory = te(i->element);
+	Function next_trajectory = te(next->element);
 	
 	//std::cout << "Updating swaps for: " << i_trajectory << ", " << next_trajectory << std::endl;
 	//std::cout << "Difference:         " << next_trajectory - i_trajectory << std::endl;
@@ -178,8 +211,8 @@
 }
 
 /* SwapEvent */
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
-class KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::SwapEvent: public Simulator::Event
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
+class KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::SwapEvent: public Simulator::Event
 {
 	public:
 		typedef						typename Simulator::Event					Parent;
@@ -192,29 +225,28 @@
 		virtual bool				process(Simulator* s) const;
 		void						set_position(iterator i)					{ pos_ = i; }
 		iterator					position() const							{ return pos_; }
-		std::ostream&				print(std::ostream& out) const;
+		std::ostream&				operator<<(std::ostream& out) const;
 
 	private:
 		KineticSort*				sort_;
 		iterator					pos_;
 };
 
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
 bool
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::SwapEvent::
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::SwapEvent::
 process(Simulator* s) const
 { 
-	std::cout << "Swapping. Current time: " << s->current_time() << std::endl;
+	rLog(rlKineticSortProcess, "Swapping. Current time: %s", tostring(s->current_time()).c_str());
 	sort_->swap(pos_, s); 
 	return true; 
 }
 
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
 std::ostream&				
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::SwapEvent::
-print(std::ostream& out) const
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::SwapEvent::
+operator<<(std::ostream& out) const
 {
-	Parent::print(out) << ", SwapEvent at " << TrajectoryExtractor_()(position()->element);
+	Parent::operator<<(out) << "SwapEvent at " << TrajectoryExtractor_()(position()->element);
 	return out;
 }
-
--- a/include/geometry/polynomial.h	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/geometry/polynomial.h	Fri Apr 04 21:14:44 2008 -0400
@@ -8,6 +8,7 @@
 #include <synaps/usolve/bezier/SlvBzStd.h>
 //#include <synaps/usolve/Sturm.h>
 //#include <synaps/arithm/Infinity.h>
+#include <synaps/arithm/gmp.h>
 
 #include <stack>
 
@@ -20,30 +21,34 @@
 {
 	public:
 		typedef						typename SynapsTraits<T>::Polynomial				Polynomial;
-		typedef						RationalFunction<Polynomial>						RationalFunction;
+		typedef						RationalFunction<Polynomial>						RationalFunctionP;
+        typedef                     RationalFunctionP                                   Function;
 
 		typedef						typename SynapsTraits<T>::Solver					Solver;
 		typedef						typename SynapsTraits<T>::RootType					RootType;
 		typedef						std::stack<RootType>								RootStack;
 
-		static void					solve(const RationalFunction& rf, RootStack& stack);
+		static void					solve(const RationalFunctionP& rf, RootStack& stack);
 		static RootType				root(const T& r)									{ return SynapsTraits<T>::root(r); }
-		static int					sign_at(const RationalFunction& rf, const RootType& r);
+		static int					sign_at(const RationalFunctionP& rf, const RootType& r);
 		static RootType				between(const RootType& r1, const RootType& r2)		{ return SynapsTraits<T>::between(r1,r2); }
+		static int					sign_at_negative_infinity(const RationalFunctionP& rf);
 };
 
 template<class T>
 struct SynapsTraits				///< Suitable for double
 {
 		typedef						T													CoefficientType;
+		typedef						T													RootType;
 		typedef						SYNAPS::UPolDse<CoefficientType>					Polynomial;
-		typedef						SYNAPS::SlvBzStd<CoefficientType>					Solver;
-		typedef						T													RootType;
+		typedef						SYNAPS::SlvBzStd<RootType, CoefficientType>			Solver;
+		typedef						Polynomial											SolverPolynomial;
 
 		static RootType				root(CoefficientType r)								{ return r; }
 		static unsigned int			multiplicity(RootType r)							{ return 1; }
 		static int					sign_at(const Polynomial& p, RootType r)			{ return SYNAPS::UPOLDAR::sign_at(p, r); }
 		static RootType				between(RootType r1, RootType r2)					{ return (r1 + r2)/2; }
+		static SolverPolynomial		convert(const Polynomial& p)						{ return p; }
 };
 
 template<>
@@ -51,7 +56,8 @@
 {
 		typedef						ZZ													CoefficientType;
 		typedef						SYNAPS::UPolDse<CoefficientType>					Polynomial;
-		typedef						SYNAPS::Algebraic<CoefficientType>					Solver;
+		typedef						Polynomial											SolverPolynomial;
+		typedef						SYNAPS::Algebraic<ZZ>								Solver;
 		typedef						Solver::root_t										RootType;
 
 		static RootType				root(const CoefficientType& r)						{ CoefficientType p[2] = {-r, 1}; return SYNAPS::solve(Polynomial(2, p), Solver(), 0);}
@@ -59,6 +65,24 @@
 		static int					sign_at(const Polynomial& p, 
 											const RootType& r)							{ return SYNAPS::ALGEBRAIC::sign_at(p, r); }
 		//static RootType				between(const RootType& r1, const RootType& r2) 	{ RootType r = r1; r += r2; r /= root(2); return r; }
+		static SolverPolynomial		convert(const Polynomial& p)						{ return p; }
+};
+
+template<>
+struct SynapsTraits<QQ>
+{
+		typedef						QQ													CoefficientType;
+		typedef						SYNAPS::UPolDse<CoefficientType>					Polynomial;
+		typedef						SYNAPS::Algebraic<ZZ>								Solver;
+		typedef						SYNAPS::UPolDse<ZZ>									SolverPolynomial;
+		typedef						Solver::root_t										RootType;
+
+		static RootType				root(const CoefficientType& r);
+		static unsigned int			multiplicity(const RootType& r)						{ return r.multiplicity(); }
+		static int					sign_at(const Polynomial& p, 
+											const RootType& r)							{ return SYNAPS::ALGEBRAIC::sign_at(convert(p), r); }
+		//static RootType				between(const RootType& r1, const RootType& r2) 	{ RootType r = r1; r += r2; r /= root(2); return r; }
+		static SolverPolynomial		convert(const Polynomial& p);
 };
 
 #include "polynomial.hpp"
--- a/include/geometry/polynomial.hpp	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/geometry/polynomial.hpp	Fri Apr 04 21:14:44 2008 -0400
@@ -1,12 +1,16 @@
+#include <utilities/log.h>
+
 template<class T>
 void
 UPolynomial<T>::
-solve(const RationalFunction& rf, RootStack& stack)
+solve(const RationalFunctionP& rf, RootStack& stack)
 {
 	typedef SYNAPS::Seq<RootType>		RootSeq;
 
-	RootSeq seq_num = SYNAPS::solve(rf.numerator(), Solver());
-	RootSeq seq_den = SYNAPS::solve(rf.denominator(), Solver());
+    AssertMsg(stack.empty(), "Stack must be empty before solve");
+
+	RootSeq seq_num = SYNAPS::solve(SynapsTraits<T>::convert(rf.numerator()), Solver());
+	RootSeq seq_den = SYNAPS::solve(SynapsTraits<T>::convert(rf.denominator()), Solver());
 
 	// TODO: assert that all roots in seq_den have positive multiplicity
 	// TODO: deal with multiplicities for the numerator
@@ -28,7 +32,45 @@
 template<class T>
 int
 UPolynomial<T>::
-sign_at(const RationalFunction& rf, const RootType& r)
+sign_at(const RationalFunctionP& rf, const RootType& r)
 {
 	return SynapsTraits<T>::sign_at(rf.numerator(), r) * SynapsTraits<T>::sign_at(rf.denominator(), r);
 }
+
+template<class T>
+int
+UPolynomial<T>::
+sign_at_negative_infinity(const RationalFunctionP& rf)
+{
+	const Polynomial& num = rf.numerator();
+	const Polynomial& den = rf.denominator();
+	int ndegree = num.get_degree();
+	int ddegree = den.get_degree();
+    if (ndegree == -1) return 0;          // ndegree == -1 => num == 0, and 0 is 0 at -infinity
+	return !((((ndegree + 1) % 2 == 0) ^ (num[ndegree] > 0)) ^
+		     (((ddegree + 1) % 2 == 0) ^ (den[ddegree] > 0))) ? 1:-1;
+}
+
+SynapsTraits<QQ>::RootType
+SynapsTraits<QQ>::
+root(const CoefficientType& r)
+{
+	ZZ p[2] = { -SYNAPS::numerator(r), SYNAPS::denominator(r) };
+	return SYNAPS::solve(SolverPolynomial(2, p), Solver(), 0); 
+}
+
+SynapsTraits<QQ>::SolverPolynomial
+SynapsTraits<QQ>::
+convert(const Polynomial& p)
+{
+	SolverPolynomial result(1, p.size() - 1);
+	assert(result.size() == p.size());
+	ZZ denominator_product = 1;
+	for (unsigned int i = 0; i < p.size(); ++i)
+		denominator_product *= SYNAPS::denominator(p[i]);
+	for (unsigned int i = 0; i < p.size(); ++i)
+		result[i] = SYNAPS::numerator(p[i]) * denominator_product / 
+					SYNAPS::denominator(p[i]);
+	return result;
+}
+
--- a/include/geometry/simulator.h	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/geometry/simulator.h	Fri Apr 04 21:14:44 2008 -0400
@@ -2,69 +2,84 @@
 #define __SIMULATOR_H__
 
 #include "utilities/eventqueue.h"
-#include "polynomial.h"
 
 template<class Comparison>  						class IndirectComparison;
-template<class PolyKernel_, class Simulator_>		class Event;
 
-template<class PolyKernel_, template<class Event> class EventComparison_ = std::less>
+/**
+ * Simulator class. Keeps a queue of events. Infinity is reached if the Event 
+ * at the front of the queue has an empty root stack. Keeps track of current time, 
+ * Event addition, and processes events one by one. Degeneracies are handled by 
+ * assuming that the FunctionKernel::Function responsible for the event must be 
+ * positive before the Event occurs.
+ *
+ * \ingroup kinetic
+ */
+template<class FuncKernel_, template<class Event> class EventComparison_ = std::less>
 class Simulator
 {
 	public:
-		typedef						PolyKernel_									PolynomialKernel;
-		typedef						typename PolynomialKernel::Polynomial		Polynomial;
-		typedef						typename PolynomialKernel::RationalFunction	RationalFunction;
-		typedef						typename PolynomialKernel::RootStack		RootStack;
-		typedef						typename PolynomialKernel::RootType			RootType;
+		typedef						FuncKernel_								    FunctionKernel;
+		typedef						typename FunctionKernel::Function	        Function;
+		typedef						typename FunctionKernel::RootStack		    RootStack;
+		typedef						typename FunctionKernel::RootType			RootType;
 		typedef						RootType									Time;
 
 		class Event;
 		typedef						EventComparison_<Event>						EventComparison;
 		typedef						EventQueue<Event*, IndirectComparison<EventComparison> >			
-																				EventQueue;
-		typedef						typename EventQueue::iterator				Key;
-		typedef						typename EventQueue::const_iterator			const_Key;
+																				EventQueueS;
+		typedef						typename EventQueueS::iterator				Key;
+		typedef						typename EventQueueS::const_iterator		const_Key;
 
 
-									Simulator(Time start = PolynomialKernel::root(0)):
-										current_(start), 
-										reached_infinity_(false)				{}
+									Simulator(Time start = FunctionKernel::root(0)):
+										current_(start)         				{}
 
 
 		template<class Event_> 
 		Key							add(const Event_& e);
 		template<class Event_> 
-		Key							add(const RationalFunction& f, const Event_& e);
+		Key							add(const Function& f, const Event_& e);
 		void						process();
-		void						update(Key k, const RationalFunction& f);
+		void						update(Key k, const Function& f);
 		
 		void						remove(Key k)								{ queue_.remove(k); }
 		Key							null_key() 									{ return queue_.end(); }
 
 		Time						current_time() const						{ return current_; }
 		Time						audit_time() const;
-		bool						reached_infinity() const					{ return reached_infinity_; }
+		bool						reached_infinity() const					{ return queue_.empty() || (*queue_.top())->root_stack().empty(); }
+        
+        Event*                      top() const                                 { return *(queue_.top()); }
+        unsigned                    size() const                                { return queue_.size(); }
 
-		std::ostream&				print(std::ostream& out) const;
+		std::ostream&				operator<<(std::ostream& out) const;
 
 	private:
 		void						update(Key i);
 
 	private:
 		Time						current_;
-		EventQueue					queue_;
-		bool						reached_infinity_;
+		EventQueueS					queue_;
 };
 
-template<class PolyKernel_, template<class Event> class EventComparison_>
-class Simulator<PolyKernel_, EventComparison_>::Event
+
+/**
+ * Base class for events. Stores a root stack, subclasses need to define process(). 
+ * Event with an empty root stack compares greater than any other Event, 
+ * pushing those events to the end of the queue.
+ */
+template<class FuncKernel_, template<class Event> class EventComparison_>
+class Simulator<FuncKernel_, EventComparison_>::Event
 {
 	public:
-		typedef						PolyKernel_									PolynomialKernel;
-		typedef						typename PolynomialKernel::RootStack		RootStack;
+		typedef						FuncKernel_									FunctionKernel;
+		typedef						typename FunctionKernel::RootStack		    RootStack;
 
+        /// process() is called when the event is at the top of the queue 
+        /// in the simulator.
 		/// Returns true if the event needs to remain in the Simulator 
-		/// (top of the root_stack() will be used for new time)
+		/// (top of the root_stack() will be used for new time).
 		virtual	bool				process(Simulator* s) const					=0;
 		
 		RootStack&					root_stack()								{ return root_stack_; }
@@ -80,15 +95,22 @@
 				return root_stack().top() < e.root_stack().top();
 		}
 
-		int							sign_before() const							{ return root_stack().top().sign_low(); }
-		int							sign_after() const							{ return root_stack().top().sign_high(); }
-
-		virtual std::ostream&		print(std::ostream& out) const				{ return out << "Event with " << root_stack_.size() << " roots"; }
+		virtual std::ostream&		operator<<(std::ostream& out) const			
+        { 
+            out << "Event with " << root_stack_.size() << " roots"; 
+            if (!root_stack_.empty()) out << "; top root: " << root_stack_.top();
+            out << ", ";
+            return out;
+        }
 
 	private:
 		RootStack					root_stack_;
 };
 
+/**
+ * Compares elements pointed at by its arguments using the provided Comparison_ 
+ * (which must not take any arguments during construction).
+ */
 template<class Comparison_>
 class IndirectComparison: public std::binary_function<const typename Comparison_::first_argument_type*, 
 													  const typename Comparison_::second_argument_type*, 
--- a/include/geometry/simulator.hpp	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/geometry/simulator.hpp	Fri Apr 04 21:14:44 2008 -0400
@@ -1,82 +1,138 @@
-template<class PolyKernel_, template<class Event> class EventComparison_>
+#include "utilities/log.h"
+#include "utilities/counter.h"
+
+#ifdef LOGGING
+static rlog::RLogChannel* rlSimulator =             DEF_CHANNEL("geometry/simulator", rlog::Log_Debug);
+
+#endif // LOGGING
+
+#ifdef COUNTERS
+static Counter*  cSimulatorProcess =                GetCounter("simulator/process");
+#endif // COUNTERS
+
+
+template<class FuncKernel_, template<class Event> class EventComparison_>
 template<class Event_>
-typename Simulator<PolyKernel_, EventComparison_>::Key
-Simulator<PolyKernel_, EventComparison_>::
+typename Simulator<FuncKernel_, EventComparison_>::Key
+Simulator<FuncKernel_, EventComparison_>::
 add(const Event_& e)
 {
 	Event* ee = new Event_(e);
 	return queue_.push(ee);
 }
 
-template<class PolyKernel_, template<class Event> class EventComparison_>
+template<class FuncKernel_, template<class Event> class EventComparison_>
 template<class Event_>
-typename Simulator<PolyKernel_, EventComparison_>::Key
-Simulator<PolyKernel_, EventComparison_>::
-add(const RationalFunction& f, const Event_& e)
+typename Simulator<FuncKernel_, EventComparison_>::Key
+Simulator<FuncKernel_, EventComparison_>::
+add(const Function& f, const Event_& e)
 {
 	Event* ee = new Event_(e);
-	//std::cout << "Solving: " << f << std::endl;
-	PolynomialKernel::solve(f, ee->root_stack());
+	rLog(rlSimulator, "Solving: %s", tostring(f).c_str());
+	int sign = FunctionKernel::sign_at_negative_infinity(f);        // going to be sign after current time
+    rLog(rlSimulator, "Sign at -infinity: %i", sign);
+    if (sign != 0)
+    {
+    	FunctionKernel::solve(f, ee->root_stack());
+        rLog(rlSimulator, "Got solution with root stack size: %i", ee->root_stack().size());
+    }
+
 	while (!ee->root_stack().empty() && ee->root_stack().top() < current_time())
+	{
 		ee->root_stack().pop();
-	//std::cout << "Pushing: " << ee->root_stack().top() << std::endl;
+		sign *= -1;
+	}
+    if (sign == -1)
+    {
+        //AssertMsg(ee->root_stack().top() == current_time(), 
+        //          "If sign is negative, we must be in the degenerate case");
+        rLog(rlSimulator, "Popping the root because of negative sign (degeneracy)");
+        ee->root_stack().pop();
+    }
+
+	if (ee->root_stack().empty())
+        rLog(rlSimulator, "Pushing event with empty root stack");
+    else
+    {
+        rLog(rlSimulator, "Root stack size: %i", ee->root_stack().size());
+        rLog(rlSimulator, "Pushing: %s", tostring(ee->root_stack().top()).c_str());
+    }
 	return queue_.push(ee);
 }
 		
-template<class PolyKernel_, template<class Event> class EventComparison_>
+template<class FuncKernel_, template<class Event> class EventComparison_>
 void
-Simulator<PolyKernel_, EventComparison_>::
-update(Key k, const RationalFunction& f)
+Simulator<FuncKernel_, EventComparison_>::
+update(Key k, const Function& f)
 {
 	Event* ee = *k;
 	ee->root_stack() = RootStack();								// no clear() in std::stack
-	PolynomialKernel::solve(f, ee->root_stack());
+	FunctionKernel::solve(f, ee->root_stack());
 	while (!ee->root_stack().empty() && ee->root_stack().top() < current_time())
 		ee->root_stack().pop();
 	update(k);
 }
 
-template<class PolyKernel_, template<class Event> class EventComparison_>
+template<class FuncKernel_, template<class Event> class EventComparison_>
 void
-Simulator<PolyKernel_, EventComparison_>::
+Simulator<FuncKernel_, EventComparison_>::
 process()
 {
-	std::cout << "Queue size: " << queue_.size() << std::endl;
+    Count(cSimulatorProcess);
+    if (reached_infinity()) return;
+	rLog(rlSimulator, "Queue size: %i", queue_.size());
 	Key top = queue_.top();
 	Event* e = *top;
+    rLog(rlSimulator, "Processing event: %s", intostring(*e).c_str());
 	
-	if (e->root_stack().empty()) 		{ reached_infinity_ = true; return; }
-	else 								{ current_ = e->root_stack().top(); e->root_stack().pop();  }
+	current_ = e->root_stack().top(); e->root_stack().pop();
 	
-	if (e->process(this))				update(top);
-	else								{ queue_.pop(); delete e; }
+    // Get the top element out of the queue, put it back depending on what process() says
+    EventQueueS tmp; tmp.prepend(top, queue_);
+
+	if (e->process(this))				{ queue_.prepend(top, tmp); update(top); }
+	else								{ delete e; }
 }
 
-template<class PolyKernel_, template<class Event> class EventComparison_>
+template<class FuncKernel_, template<class Event> class EventComparison_>
 void
-Simulator<PolyKernel_, EventComparison_>::
+Simulator<FuncKernel_, EventComparison_>::
 update(Key i)
 {
 	queue_.update(i);
 }
 		
-template<class PolyKernel_, template<class Event> class EventComparison_>
-typename Simulator<PolyKernel_, EventComparison_>::Time
-Simulator<PolyKernel_, EventComparison_>::
+template<class FuncKernel_, template<class Event> class EventComparison_>
+typename Simulator<FuncKernel_, EventComparison_>::Time
+Simulator<FuncKernel_, EventComparison_>::
 audit_time() const
 {
 	const_Key top = queue_.top();
 	Event* e = *top;
 
 	if (e->root_stack().empty()) return current_ + 1;
-	else return PolynomialKernel::between(e->root_stack().top(), current_);
+	else return FunctionKernel::between(e->root_stack().top(), current_);
 }
 		
-template<class PolyKernel_, template<class Event> class EventComparison_>
+template<class FuncKernel_, template<class Event> class EventComparison_>
 std::ostream&
-Simulator<PolyKernel_, EventComparison_>::
-print(std::ostream& out) const
+Simulator<FuncKernel_, EventComparison_>::
+operator<<(std::ostream& out) const
 {
 	out << "Simulator: " << std::endl;
 	return queue_.print(out, "  ");
 }
+
+template<class FuncKernel_, template<class Event> class EventComparison_>
+std::ostream&
+operator<<(std::ostream& out, const Simulator<FuncKernel_, EventComparison_>& s)
+{
+    return s.operator<<(out);
+}
+
+template<class FuncKernel_, template<class Event> class EventComparison_>
+std::ostream&
+operator<<(std::ostream& out, const typename Simulator<FuncKernel_, EventComparison_>::Event& e)
+{
+    return e.operator<<(out);
+}
--- a/include/topology/conesimplex.h	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/conesimplex.h	Fri Apr 04 21:14:44 2008 -0400
@@ -9,6 +9,9 @@
 #include <list>
 #include <iostream>
 
+#include "utilities/types.h"
+
+
 template<class S>
 class ConeSimplex: public S
 {
@@ -18,12 +21,18 @@
 		typedef		std::list<Self>										Cycle;
 
     public:
+								ConeSimplex(const Self& s): 
+                                    Parent(s), coned_(s.coned_)         {}
 								ConeSimplex(const Parent& parent, 
 											bool coned = false):
 									Parent(parent), coned_(coned)		{}
 	    
 		Cycle					boundary() const;
 		bool					coned() const							{ return coned_; }
+        Dimension               dimension() const                       { return coned_ ? (Parent::dimension() + 1) : Parent::dimension(); }
+        
+        bool                    operator<(const Self& other) const      { if (coned_ ^ other.coned_) return !coned_; else return Parent::operator<(other); }
+        bool                    operator==(const Self& other) const     { return !(coned_ ^ other.coned_) && Parent::operator==(other); }
 
 		std::ostream& 			operator<<(std::ostream& out) const;
 		
--- a/include/topology/cycle.h	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/cycle.h	Fri Apr 04 21:14:44 2008 -0400
@@ -16,6 +16,8 @@
  * 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>
--- a/include/topology/cycle.hpp	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/cycle.hpp	Fri Apr 04 21:14:44 2008 -0400
@@ -166,7 +166,8 @@
 {
 	for (const_iterator cur = begin(); cur != end(); ++cur)
 	{
-		out << **cur << ", ";
+        if (cur != begin()) out << ", ";
+		out << **cur;
 	}
 	// out << "(last: " << *last << ")";  // For debugging only
 	return out;
--- a/include/topology/filtration.h	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/filtration.h	Fri Apr 04 21:14:44 2008 -0400
@@ -22,6 +22,8 @@
  * 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
@@ -37,6 +39,8 @@
 		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 
@@ -61,7 +65,7 @@
 		/// 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 = true);
+		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
--- a/include/topology/filtration.hpp	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/filtration.hpp	Fri Apr 04 21:14:44 2008 -0400
@@ -117,6 +117,15 @@
 	
 	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");
--- a/include/topology/filtrationcontainer.h	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/filtrationcontainer.h	Fri Apr 04 21:14:44 2008 -0400
@@ -13,6 +13,8 @@
  * FiltrationContainer class. Serves as a parent of Filtration that 
  * describes the container functionality. Used by FiltrationSimplex 
  * to get Cycle representation.
+ *
+ * \ingroup topology
  */
 template<class FltrSmplx>
 class FiltrationContainer: public ConsistencyList<FltrSmplx>
--- a/include/topology/filtrationsimplex.h	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/filtrationsimplex.h	Fri Apr 04 21:14:44 2008 -0400
@@ -21,6 +21,8 @@
 /**
  * Evaluator is a base class for the structures that are able to return a value
  * given a simplex.
+ *
+ * \ingroup topology
  */
 template<class Smplx>
 class Evaluator
@@ -37,6 +39,8 @@
 /**
  * FiltrationSimplex stores information needed for the RU-decomposition: 
  * cycle (column of R), trail (row of U), and pair.
+ *
+ * \ingroup topology
  */
 template<class Smplx>
 class FiltrationSimplex: public Smplx
--- a/include/topology/simplex.h	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/simplex.h	Fri Apr 04 21:14:44 2008 -0400
@@ -20,6 +20,8 @@
  * 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.
+ *
+ * \ingroup topology
  */
 template<class V>
 class SimplexWithVertices
@@ -83,6 +85,8 @@
 
 /**
  * SimplexWithValue explicitly adds a RealType value to the SimplexWithVertices.
+ *
+ * \ingroup topology
  */
 template<class Vert>
 class SimplexWithValue: public SimplexWithVertices<Vert>
@@ -130,6 +134,8 @@
 
 /**
  * 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>
--- a/include/topology/vineyard.h	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/vineyard.h	Fri Apr 04 21:14:44 2008 -0400
@@ -23,6 +23,8 @@
 /**
  * Vineyard class. Keeps track of vines and knees. switched() is the key function called
  * by filtration when pairing switches after a Filtration::transpose().
+ *
+ * \ingroup topology
  */
 template<class FltrSmplx>
 class Vineyard
@@ -32,8 +34,9 @@
 		typedef							typename FiltrationSimplex::Simplex			Simplex;
 		typedef							Vine<Simplex>								VineS;
 		typedef							Knee<Simplex>								KneeS;
-		typedef							std::list<VineS>							VineList;
-		typedef							std::vector<VineList>						VineListVector;
+		typedef							std::list<VineS>						    VineList;
+		typedef							std::list<VineList>						    VineListList;
+        typedef                         std::vector<typename VineListList::iterator> VineListVector;
 		typedef							typename FiltrationSimplex::Cycle			Cycle;
 
 		typedef							typename FiltrationSimplex::Index			Index;
@@ -45,6 +48,7 @@
 
 		void							start_vines(Index bg, Index end);
 		void							switched(Index i, Index j);
+		void							record_knee(Index i);
 		void							record_diagram(Index bg, Index end);
 
 		void							set_evaluator(Evaluator* eval)				{ evaluator = eval; }
@@ -55,16 +59,18 @@
 		typename KneeS::SimplexList  	resolve_cycle(Index i) const;
 
 	private:
-		void							record_knee(Index i);
 		void							start_vine(Index i);
 
 	private:
-		VineListVector					vines;
+		VineListList                    vines;            // stores vine lists
+		VineListVector                  vines_vector;     // stores pointers (iterators) to vine lists
 		Evaluator*						evaluator;
 };
 
 /**
  * Knee class stores the knee in R^3 as well as the cycle that is associated with the Simplex starting from the Knee.
+ *
+ * \ingroup topology
  */
 template<class S>
 class Knee
@@ -117,11 +123,17 @@
 		typedef					typename VineRepresentation::const_iterator		const_knee_iterator;
 		
 								Vine()											{}
-								Vine(const KneeS& k)							{ add(k); }
+								Vine(const Vine& other): 
+                                    VineRepresentation(other)	                {}
+								Vine(const VineRepresentation& other): 
+                                    VineRepresentation(other)	                {}
+								Vine(const KneeS& k)						    { add(k); }
 		
 		void 					add(RealType b, RealType d, RealType t)			{ push_back(KneeS(b,d,t)); }
 		void 					add(const KneeS& k)								{ push_back(k); }
 
+        std::ostream&           operator<<(std::ostream& out) const             { for (const_knee_iterator cur = begin(); cur != end(); ++cur) out << *cur; return out; }
+
 		using VineRepresentation::begin;
 		using VineRepresentation::end;
 		using VineRepresentation::front;
@@ -139,6 +151,9 @@
 		void 					serialize(Archive& ar, version_type );
 };
 
+template<class S>
+std::ostream& operator<<(std::ostream& out, const Vine<S>& v) 					{ return v.operator<<(out); }
+
 
 #include "vineyard.hpp"
 
--- a/include/topology/vineyard.hpp	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/vineyard.hpp	Fri Apr 04 21:14:44 2008 -0400
@@ -23,7 +23,8 @@
 		if (dim >= vines.size())
 		{
 			AssertMsg(dim == vines.size(), "New dimension has to be contiguous");
-			vines.push_back(std::list<VineS>());
+			vines.push_back(VineList());
+            vines_vector.push_back(boost::prior(vines.end()));
 		}
 
 		start_vine(cur);
@@ -60,8 +61,8 @@
 	AssertMsg(i->sign(), "Can only start vines for positive simplices");
 		
 	Dimension dim = i->dimension();
-	vines[dim].push_back(VineS());
-	i->set_vine(&vines[dim].back());
+	vines_vector[dim]->push_back(VineS());
+	i->set_vine(&vines_vector[dim]->back());
 	i->pair()->set_vine(i->vine());
 }
 	
@@ -88,12 +89,12 @@
 Vineyard<FS>::
 save_edges(const std::string& filename) const
 {
-	for (unsigned int i = 0; i < vines.size(); ++i)
+	for (unsigned int i = 0; i < vines_vector.size(); ++i)
 	{
 		std::ostringstream os; os << i;
 		std::string fn = filename + os.str() + ".edg";
 		std::ofstream out(fn.c_str());
-		for (typename VineList::const_iterator vi = vines[i].begin(); vi != vines[i].end(); ++vi)
+		for (typename VineList::const_iterator vi = vines_vector[i]->begin(); vi != vines_vector[i]->end(); ++vi)
 			for (typename VineS::const_iterator ki = vi->begin(), kiprev = ki++; ki != vi->end(); kiprev = ki++)
 			{
 				if (kiprev->is_infinite() || ki->is_infinite()) continue;
@@ -122,6 +123,7 @@
 		rLog(rlVineyard, "Creating knee");
 		KneeS k(evaluator->value(*i), evaluator->value(*(i->pair())), evaluator->time());
 		rLog(rlVineyard, "Knee created: %s", tostring(k).c_str());
+        rLog(rlVineyard, "Vine: %s", tostring(*(i->vine())).c_str());
 
 		if (!k.is_diagonal() || i->vine()->empty())			// non-diagonal k, or empty vine
 		{
--- a/include/utilities/consistencylist.h	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/utilities/consistencylist.h	Fri Apr 04 21:14:44 2008 -0400
@@ -144,7 +144,7 @@
 	T 				data;
 	OrderType		consistency;
 
-	std::ostream& 		operator<<(std::ostream& out) const					{ return out << data << ": " << consistency; }
+	std::ostream& 		operator<<(std::ostream& out) const					{ return out << consistency << ": " << data; }
 };
 
 template<class T>
--- a/include/utilities/eventqueue.h	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/utilities/eventqueue.h	Fri Apr 04 21:14:44 2008 -0400
@@ -8,6 +8,7 @@
 #include <iostream>
 #include <cassert>					// TODO: switch to internal debugging system
 #include <string>
+#include <algorithm>
 
 // TODO: change inefficient list-based implementation to something heap-based
 // Need a queue that supports deleting arbitrary items (given by iterator), 
@@ -34,6 +35,9 @@
 		void 					pop()						{ assert(!empty()); queue_.erase(queue_.begin()); }
 		void					remove(iterator i)			{ queue_.erase(i); }
 		void					update(iterator i);
+        void                    prepend(iterator i, 
+                                        EventQueue& other)  { queue_.splice(queue_.begin(), other.queue_, i); }
+        ///< intended for temporary storage of elements from other queues
 
 		iterator 				end()						{ return queue_.end(); }
 		const_iterator 			end() const					{ return queue_.end(); }
@@ -54,7 +58,8 @@
 {
 	QueueRepresentation tmp;
 	tmp.splice(tmp.end(), queue_, i);
-	iterator pos = std::find_if(queue_.begin(), queue_.end(), std::not1(std::bind2nd(EventComparison(), *i)));
+	//iterator pos = std::find_if(queue_.begin(), queue_.end(), std::not1(std::bind2nd(EventComparison(), *i)));
+	iterator pos = std::find_if(queue_.begin(), queue_.end(), std::bind1st(EventComparison(), *i));
 	queue_.splice(pos, tmp);
 }
 
@@ -64,7 +69,7 @@
 print(std::ostream& out, const std::string& prefix) const
 {
 	for (typename QueueRepresentation::const_iterator cur = queue_.begin(); cur != queue_.end(); ++cur)
-		(*cur)->print(out << prefix) << std::endl;
+		(*cur)->operator<<(out << prefix) << std::endl;
 	return out;
 }
 
--- a/include/utilities/log.h	Fri Apr 04 16:17:28 2008 -0400
+++ b/include/utilities/log.h	Fri Apr 04 21:14:44 2008 -0400
@@ -12,6 +12,8 @@
 
 template<class T>
 std::string tostring(const T& t) { std::ostringstream out; out << t; return out.str(); }
+template<class T>
+std::string intostring(const T& t) { std::ostringstream out; t.operator<<(out); return out.str(); }
 
 #define AssertMsg(cond, message, ...)		do { if (!(cond)) { rError(message, ##__VA_ARGS__); rAssertSilent(cond); } } while (0)
 	
--- a/tests/CMakeLists.txt	Fri Apr 04 16:17:28 2008 -0400
+++ b/tests/CMakeLists.txt	Fri Apr 04 21:14:44 2008 -0400
@@ -1,1 +1,2 @@
+add_subdirectory			(geometry)
 add_subdirectory			(utilities)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/geometry/CMakeLists.txt	Fri Apr 04 21:14:44 2008 -0400
@@ -0,0 +1,16 @@
+set							(targets
+							 euclidean
+							 test-eventqueue)
+
+if                          (use_synaps)
+    set                     (targets                    ${targets}
+                                                        polynomial
+                                                        test-kinetic-sort
+							                            test-linalg)
+endif                       (use_synaps)
+
+foreach 					(t ${targets})
+	add_executable			(${t} ${t}.cpp)
+	target_link_libraries	(${t} ${synaps_libraries} ${libraries})
+endforeach 					(t ${targets})
+
--- a/tests/geometry/Makefile	Fri Apr 04 16:17:28 2008 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,26 +0,0 @@
-SYNAPS_LIBS = -L/usr/lib -L/usr/lib -lsynaps -llapack -lblas -lgmpxx -lgmp -lmps -L/usr/lib/gcc/i686-pc-linux-gnu/4.1.2 -L/usr/lib/gcc/i686-pc-linux-gnu/4.1.2/../../.. -lgfortranbegin -lgfortran -lm -lgcc_s
-
-INCLUDE_PATH=-I../../include/geometry
-DEFINES=-DBOOST_UBLAS_TYPE_CHECK=0 -g
-LIBRARIES=-lsynaps -lgmp -lgmpxx #${SYNAPS_LIBS}
-#CPPFLAGS=-march=i686 -mtune=native -Wall
-#CPPFLAGS=-march=i686 -mtune=generic -Wall
-CPPFLAGS=-Wall
-GCC=g++
-
-all: test-linalg euclidean polynomial test-eventqueue test-kinetic-sort
-
-test-linalg: test-linalg.cpp
-	${GCC} test-linalg.cpp -o test-linalg ${INCLUDE_PATH} ${DEFINES} ${LIBRARIES} ${CPPFLAGS}
-
-euclidean: euclidean.cpp
-	${GCC} euclidean.cpp -o euclidean ${INCLUDE_PATH} ${DEFINES} ${CPPFLAGS}
-
-polynomial: polynomial.cpp
-	${GCC} polynomial.cpp -o polynomial ${INCLUDE_PATH} ${DEFINES} ${LIBRARIES} ${CPPFLAGS}
-
-test-eventqueue: test-eventqueue.cpp
-	${GCC} test-eventqueue.cpp -o test-eventqueue ${CPPFLAGS} ${DEFINES} ${INCLUDE_PATH}
-
-test-kinetic-sort: test-kinetic-sort.cpp
-	${GCC} test-kinetic-sort.cpp -o test-kinetic-sort ${INCLUDE_PATH} ${DEFINES} ${LIBRARIES} ${CPPFLAGS}
--- a/tests/geometry/euclidean.cpp	Fri Apr 04 16:17:28 2008 -0400
+++ b/tests/geometry/euclidean.cpp	Fri Apr 04 21:14:44 2008 -0400
@@ -1,4 +1,4 @@
-#include "euclidean.h"
+#include "geometry/euclidean.h"
 #include <vector>
 #include <iostream>
 #include <cmath>
--- a/tests/geometry/polynomial.cpp	Fri Apr 04 16:17:28 2008 -0400
+++ b/tests/geometry/polynomial.cpp	Fri Apr 04 21:14:44 2008 -0400
@@ -1,5 +1,5 @@
-#include <euclidean.h>
-#include <polynomial.h>
+#include <geometry/euclidean.h>
+#include <geometry/polynomial.h>
 
 #include <vector>
 #include <iostream>
@@ -7,7 +7,7 @@
 
 typedef UPolynomial<ZZ>						PolynomialKernel;
 typedef PolynomialKernel::Polynomial		Polynomial;
-typedef PolynomialKernel::RationalFunction	RationalF;
+typedef PolynomialKernel::Function	        RationalF;
 
 typedef Kernel<RationalF>					K;
 typedef K::Point							Point;
@@ -75,8 +75,6 @@
 		while (!roots.empty()) { std::cout << roots.top() << std::endl; roots.pop(); }
 	}
 
-	return 0;
-
 	// Edges
 	{
 		PointContainer vertices(2);
--- a/tests/geometry/test-eventqueue.cpp	Fri Apr 04 16:17:28 2008 -0400
+++ b/tests/geometry/test-eventqueue.cpp	Fri Apr 04 21:14:44 2008 -0400
@@ -1,4 +1,4 @@
-#include <eventqueue.h>
+#include <utilities/eventqueue.h>
 #include <functional>
 #include <iostream>
 
--- a/tests/geometry/test-kinetic-sort.cpp	Fri Apr 04 16:17:28 2008 -0400
+++ b/tests/geometry/test-kinetic-sort.cpp	Fri Apr 04 21:14:44 2008 -0400
@@ -1,25 +1,29 @@
-#include <polynomial.h>
-#include <simulator.h>
-#include <kinetic-sort.h>
+#include <geometry/polynomial.h>
+#include <geometry/simulator.h>
+#include <geometry/kinetic-sort.h>
 #include <iostream>
 
 #include <boost/utility.hpp>
+#include <boost/bind.hpp>
 
-//typedef		double							FieldType;
+#include <utilities/log.h>
+
+typedef		double							FieldType;
 //typedef		ZZ								FieldType;
-typedef		QQ								FieldType;
+//typedef		QQ								FieldType;
 typedef 	UPolynomial<FieldType>			PolyKernel;
 typedef		PolyKernel::Polynomial			Polynomial;
 typedef 	std::list<Polynomial>			SortDS;
+typedef 	SortDS::iterator				SortDSIterator;
 typedef		Simulator<PolyKernel>			SimulatorFT;
 
 class TrajectoryExtractor
 {
 	public:
-		Polynomial		operator()(SortDS::iterator i) const			{ return *i; }
+		Polynomial		operator()(SortDSIterator i) const			{ return *i; }
 };
 
-typedef		KineticSort<SortDS, TrajectoryExtractor, SimulatorFT>		KineticSortDS;
+typedef		KineticSort<SortDSIterator, TrajectoryExtractor, SimulatorFT>		KineticSortDS;
 
 struct EvaluatedComparison: public std::binary_function<const Polynomial&, const Polynomial&, bool>
 {
@@ -29,21 +33,33 @@
 	FieldType					vv;
 };
 
-void swap(SortDS* s, SortDS::iterator i)
+void swap(SortDS* s, SortDSIterator i, SimulatorFT* simulator)
 {
 	std::cout << "Swapping " << *i << " " << *boost::next(i) << std::endl;
 	s->splice(i, *s, boost::next(i));
 }
 
-int main()
+int main(int argc, char** argv)
 {
+#ifdef LOGGING
+    rlog::RLogInit(argc, argv);
+
+    stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+    stdoutLog.subscribeTo( RLOG_CHANNEL("geometry/simulator") );
+    stdoutLog.subscribeTo( RLOG_CHANNEL("geometry/kinetic-sort") );
+#endif
+
 	SimulatorFT		simulator;
 	SortDS 			list;
 
 	// Insert polynomials and sort the list for current time
 	list.push_back(Polynomial("x^3 - 3"));
 	list.push_back(Polynomial("x^2 - 2*x - 2"));
+	list.push_back(Polynomial("x^2 - 2*x - 2"));
 	list.push_back(Polynomial("2*x - 4"));
+	list.push_back(Polynomial("x"));
+	list.push_back(Polynomial("-x + 4"));
+	list.push_back(Polynomial("2"));
 	//list.sort(EvaluatedComparison(simulator.current_time()));
 	list.sort(EvaluatedComparison(0));
 
@@ -52,15 +68,23 @@
 		std::cout << *cur << std::endl;
 
 	// Setup kinetic sort
-	KineticSortDS	ks(&list, &simulator, swap);
+	KineticSortDS	ks(list.begin(), list.end(), boost::bind(swap, &list, _1, _2), &simulator);
 
-	while(!simulator.reached_infinity() && simulator.current_time() < 1)
+    std::cout << "Examining " << simulator;
+
+	while(!simulator.reached_infinity())
 	{
-		//std::cout << "Current time before: " << simulator.current_time() << std::endl;
-		if (!ks.audit(&simulator)) return 1;
-		//simulator.print(std::cout << "Auditing ");
+		std::cout << "Current time before: " << simulator.current_time() << std::endl;
+		//if (!ks.audit(&simulator)) return 1;
+		ks.audit(&simulator);
+        std::cout << "Examining " << simulator;
 		simulator.process();
-		//std::cout << "Current time after: " << simulator.current_time() << std::endl;
+		std::cout << "Current time after: " << simulator.current_time() << std::endl;
 	}
+	ks.audit(&simulator);
+    std::cout << "Examining " << simulator;
 
+    std::cout << "Done at " << simulator.current_time() << std::endl;
+    for (SortDS::const_iterator cur = list.begin(); cur != list.end(); ++cur)
+        std::cout << "  " << *cur << std::endl;    
 }
--- a/tests/geometry/test-linalg.cpp	Fri Apr 04 16:17:28 2008 -0400
+++ b/tests/geometry/test-linalg.cpp	Fri Apr 04 21:14:44 2008 -0400
@@ -1,9 +1,9 @@
-#include "linalg.h"
+#include "geometry/linalg.h"
 #include <iostream>
 
 #include <synaps/upol.h>
 #include <synaps/upol/gcd.h>
-#include "rational-function.h"
+#include "geometry/rational-function.h"
 
 
 typedef			UPolDse<double>					Polynomial;