Switched to internal KineticSort dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Sat, 06 Feb 2010 23:22:43 -0800
branchdev
changeset 199 2bde4c56101c
parent 198 e95766342e5f
child 200 73e8dce642be
Switched to internal KineticSort
examples/pl-functions/combustion-vineyard.cpp
examples/pl-functions/pl-vineyard.cpp
examples/pl-functions/test-grid2D-vineyard.cpp
include/geometry/kinetic-sort.h
include/geometry/kinetic-sort.hpp
include/geometry/simulator.h
include/geometry/simulator.hpp
include/topology/lsvineyard.h
include/topology/lsvineyard.hpp
--- a/examples/pl-functions/combustion-vineyard.cpp	Fri Feb 05 11:49:52 2010 -0800
+++ b/examples/pl-functions/combustion-vineyard.cpp	Sat Feb 06 23:22:43 2010 -0800
@@ -81,7 +81,7 @@
     std::cout << "Pairing computed" << std::endl;
     
     // Compute vineyard
-    v.compute_vineyard(g1, true);
+    v.compute_vineyard(g1);
     std::cout << "Vineyard computed" << std::endl;
 
     v.vineyard().save_edges("combustion");
--- a/examples/pl-functions/pl-vineyard.cpp	Fri Feb 05 11:49:52 2010 -0800
+++ b/examples/pl-functions/pl-vineyard.cpp	Sat Feb 06 23:22:43 2010 -0800
@@ -15,7 +15,7 @@
 namespace bl = boost::lambda;
 
 
-typedef     float                                       VertexValue;
+typedef     double                                      VertexValue;
 typedef     unsigned                                    Vertex;
 typedef     std::vector<VertexValue>                    VertexVector;
 struct SubscriptFunctor: public std::unary_function<Vertex, VertexValue>
@@ -74,7 +74,7 @@
     for (size_t i = 1; i < vertices.size(); ++i)
     {
         veval = VertexEvaluator(vertices[i]);
-        v.compute_vineyard(veval, explicit_events);
+        v.compute_vineyard(veval);
         std::cout << "Processed frame: " << i << std::endl;
     }
     std::cout << "Vineyard computed" << std::endl;
--- a/examples/pl-functions/test-grid2D-vineyard.cpp	Fri Feb 05 11:49:52 2010 -0800
+++ b/examples/pl-functions/test-grid2D-vineyard.cpp	Sat Feb 06 23:22:43 2010 -0800
@@ -42,7 +42,7 @@
         std::cout << "  " << v.filtration().simplex(cur) << std::endl;
 
     // Compute vineyard
-    v.compute_vineyard(g1, true);
+    v.compute_vineyard(g1);
     std::cout << "Vineyard computed" << std::endl;
     
     // Simplex order after
--- a/include/geometry/kinetic-sort.h	Fri Feb 05 11:49:52 2010 -0800
+++ b/include/geometry/kinetic-sort.h	Sat Feb 06 23:22:43 2010 -0800
@@ -53,7 +53,7 @@
 		/// \name Core Functionality
 		/// @{
 									KineticSort();
-									KineticSort(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator);
+									KineticSort(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator, const TrajectoryExtractor& te);
 		void						initialize(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator);
 
 		void						insert(iterator pos, ElementIterator f, ElementIterator l, Simulator* simulator);
@@ -68,6 +68,8 @@
 		iterator					begin() 									{ return list_.begin(); }
 		iterator					end() 										{ return list_.end(); }
 
+        const TrajectoryExtractor&  trajectory_extractor() const                { return te_; }
+
 	private:
 		class SwapEvent;
 		void						schedule_swaps(iterator b, iterator e, Simulator* s);
@@ -76,6 +78,7 @@
 	private:
 		NodeList					list_;
 		Swap						swap_;	
+        TrajectoryExtractor         te_;
 };
 
 #include "kinetic-sort.hpp"
--- a/include/geometry/kinetic-sort.hpp	Fri Feb 05 11:49:52 2010 -0800
+++ b/include/geometry/kinetic-sort.hpp	Sat Feb 06 23:22:43 2010 -0800
@@ -21,7 +21,8 @@
 
 template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
 KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
-KineticSort(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator)
+KineticSort(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator, const TrajectoryExtractor& te = TrajectoryExtractor()):
+    te_(te)
 {
 	initialize(b, e, swap, simulator);
 }
@@ -129,16 +130,14 @@
 	Time t = simulator->audit_time();
 	rLog(rlKineticSortAudit, "Auditing at %s", tostring(t).c_str());
 
-	TrajectoryExtractor	te;
-	
 	typename NodeList::const_iterator next = list_.begin();
 	typename NodeList::const_iterator cur = next++;
-	Function cur_trajectory = te(cur->element);
+	Function cur_trajectory = te_(cur->element);
 	while (next != list_.end())
 	{
 		rLog(rlKineticSortAudit, "  %s", intostring(**(cur->swap_event_key)).c_str());
 
-		Function next_trajectory = te(next->element);
+		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());
@@ -146,8 +145,8 @@
                                                          tostring(FunctionKernel::sign_at(next_trajectory - cur_trajectory, t)).c_str());
 		if (FunctionKernel::sign_at(next_trajectory - cur_trajectory, t) == -1)
 		{
-			rError("Audit failed at %s, %s", tostring(*cur->element).c_str(), 
-                                             tostring(*next->element).c_str());
+			// rError("Audit failed at %s, %s", tostring(*cur->element).c_str(), 
+            //                                  tostring(*next->element).c_str());
 			return false;
 		}
 
@@ -166,14 +165,12 @@
 {
 	typedef 		typename Simulator::Function		        Function;
 	
-	TrajectoryExtractor	te;
-	
 	iterator next = b; 
 	iterator cur = next++;
-	Function cur_trajectory = te(cur->element);
+	Function cur_trajectory = te_(cur->element);
 	while (next != e)
 	{
-		Function next_trajectory = te(next->element);
+		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));
@@ -197,11 +194,9 @@
 		return;
 	}
 
-	TrajectoryExtractor	te;
-	
 	iterator next = boost::next(i); 
-	Function i_trajectory = te(i->element);
-	Function 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;
@@ -247,6 +242,6 @@
 KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::SwapEvent::
 operator<<(std::ostream& out) const
 {
-	Parent::operator<<(out) << "SwapEvent at " << TrajectoryExtractor_()(position()->element);
+	Parent::operator<<(out) << "SwapEvent at " << sort_->trajectory_extractor()(position()->element);
 	return out;
 }
--- a/include/geometry/simulator.h	Fri Feb 05 11:49:52 2010 -0800
+++ b/include/geometry/simulator.h	Sat Feb 06 23:22:43 2010 -0800
@@ -1,9 +1,10 @@
 #ifndef __SIMULATOR_H__
 #define __SIMULATOR_H__
 
-#include "utilities/eventqueue.h"
+#include <utilities/eventqueue.h>
+#include <utilities/indirect.h>
 
-template<class Comparison>  						class IndirectComparison;
+#include <limits>
 
 /**
  * Simulator class. Keeps a queue of events. Infinity is reached if the Event 
@@ -17,50 +18,54 @@
 template<class FuncKernel_, template<class Event> class EventComparison_ = std::less>
 class Simulator
 {
-	public:
-		typedef						FuncKernel_								    FunctionKernel;
-		typedef						typename FunctionKernel::Function	        Function;
-		typedef						typename FunctionKernel::RootStack		    RootStack;
-		typedef						typename FunctionKernel::RootType			RootType;
-		typedef						RootType									Time;
+    public:
+        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> >			
-																				EventQueueS;
-		typedef						typename EventQueueS::iterator				Key;
-		typedef						typename EventQueueS::const_iterator		const_Key;
+        class Event;
+        typedef                     EventComparison_<Event>                     EventComparison;
+
+        class IndirectEventComparison;
+        typedef                     EventQueue<Event*, IndirectEventComparison> EventQueueS;
+        typedef                     typename EventQueueS::iterator              Key;
+        typedef                     typename EventQueueS::const_iterator        const_Key;
 
 
-									Simulator(Time start = FunctionKernel::root(0)):
-										current_(start)         				{}
+                                    Simulator(Time start = FunctionKernel::root(0)):
+                                        current_(start)                         {}
 
 
-		template<class Event_> 
-		Key							add(const Event_& e);
-		template<class Event_> 
-		Key							add(const Function& f, const Event_& e);
-		void						process();
-		void						update(Key k, const Function& f);
-		
-		void						remove(Key k)								{ queue_.remove(k); }
-		Key							null_key() 									{ return queue_.end(); }
+        template<class Event_> 
+        Key                         add(const Event_& e);
+        template<class Event_> 
+        Key                         add(const Function& f, const Event_& e);
+        void                        process();
+        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 queue_.empty() || (*queue_.top())->root_stack().empty(); }
+        Time                        current_time() const                        { return current_; }
+        Time                        audit_time() const;
+        bool                        reached_infinity() const                    { return queue_.empty() || (*queue_.top())->root_stack().empty(); }
+        Time                        next_event_time() const                     { return queue_.empty() ? std::numeric_limits<Time>::infinity():(*queue_.top())->root_stack().top(); } 
         
         Event*                      top() const                                 { return *(queue_.top()); }
         unsigned                    size() const                                { return queue_.size(); }
+        unsigned                    event_count() const                         { return count_; }
 
-		std::ostream&				operator<<(std::ostream& out) const;
+        std::ostream&               operator<<(std::ostream& out) const;
 
-	private:
-		void						update(Key i);
+    private:
+        void                        update(Key i);
 
-	private:
-		Time						current_;
-		EventQueueS					queue_;
+    private:
+        Time                        current_;
+        EventQueueS                 queue_;
+        unsigned                    count_;
 };
 
 
@@ -72,30 +77,30 @@
 template<class FuncKernel_, template<class Event> class EventComparison_>
 class Simulator<FuncKernel_, EventComparison_>::Event
 {
-	public:
-		typedef						FuncKernel_									FunctionKernel;
-		typedef						typename FunctionKernel::RootStack		    RootStack;
+    public:
+        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).
-		virtual	bool				process(Simulator* s) const					=0;
-		
-		RootStack&					root_stack()								{ return root_stack_; }
-		const RootStack&			root_stack() const							{ return root_stack_; }
+        /// Returns true if the event needs to remain in the Simulator 
+        /// (top of the root_stack() will be used for new time).
+        virtual bool                process(Simulator* s) const                 =0;
+        
+        RootStack&                  root_stack()                                { return root_stack_; }
+        const RootStack&            root_stack() const                          { return root_stack_; }
 
-		bool						operator<(const Event& e) const				
-		{ 
-			if (root_stack().empty())
-				return false;
-			else if (e.root_stack().empty())
-				return true;
-			else
-				return root_stack().top() < e.root_stack().top();
-		}
+        bool                        operator<(const Event& e) const             
+        { 
+            if (root_stack().empty())
+                return false;
+            else if (e.root_stack().empty())
+                return true;
+            else
+                return root_stack().top() < e.root_stack().top();
+        }
 
-		virtual std::ostream&		operator<<(std::ostream& out) const			
+        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();
@@ -103,28 +108,30 @@
             return out;
         }
 
-	private:
-		RootStack					root_stack_;
+    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*, 
-													  bool>
+template<class FuncKernel_, template<class Event> class EventComparison_>
+class Simulator<FuncKernel_, EventComparison_>::IndirectEventComparison: 
+    public std::binary_function<const typename EventComparison::first_argument_type*, 
+                                const typename EventComparison::second_argument_type*, 
+                                bool>
 {
-	public:
-		typedef						Comparison_											Comparison;
-		typedef						const typename Comparison::first_argument_type*		first_argument_type;
-		typedef						const typename Comparison::second_argument_type*	second_argument_type;
+    public:
+        typedef                     EventComparison                                     Comparison;
+        typedef                     const typename Comparison::first_argument_type*     first_argument_type;
+        typedef                     const typename Comparison::second_argument_type*    second_argument_type;
 
-		bool						operator()(first_argument_type e1, second_argument_type e2) const
-		{ return Comparison()(*e1, *e2); }
+        bool                        operator()(first_argument_type e1, second_argument_type e2) const
+        { return Comparison()(*e1, *e2); }
 };
 
+
 #include "simulator.hpp"
 
 #endif // __SIMULATOR_H__
--- a/include/geometry/simulator.hpp	Fri Feb 05 11:49:52 2010 -0800
+++ b/include/geometry/simulator.hpp	Sat Feb 06 23:22:43 2010 -0800
@@ -88,10 +88,12 @@
 	current_ = e->root_stack().top(); e->root_stack().pop();
 	
     // Get the top element out of the queue, put it back depending on what process() says
-    EventQueueS tmp; tmp.prepend(top, queue_);
+    EventQueueS             tmp; tmp.prepend(top, queue_);
 
 	if (e->process(this))				{ queue_.prepend(top, tmp); update(top); }
 	else								{ delete e; }
+
+    ++count_;
 }
 
 template<class FuncKernel_, template<class Event> class EventComparison_>
--- a/include/topology/lsvineyard.h	Fri Feb 05 11:49:52 2010 -0800
+++ b/include/topology/lsvineyard.h	Sat Feb 06 23:22:43 2010 -0800
@@ -1,6 +1,6 @@
 /**
  * Author: Dmitriy Morozov
- * Department of Computer Science, Duke University, 2005 -- 2009
+ * Department of Computer Science, Duke University, 2005 -- 2010
  */
 
 #ifndef __LSVINEYARD_H__
@@ -15,9 +15,9 @@
 
 #include <utilities/indirect.h>
 
-#include <CGAL/Kinetic/Inexact_simulation_traits.h>
-#include <CGAL/Kinetic/Sort.h>
-#include <CGAL/Kinetic/Sort_visitor_base.h>
+#include <geometry/simulator.h>
+#include <geometry/kinetic-sort.h>
+#include <geometry/linear-kernel.h>
 
 #include <boost/tuple/tuple.hpp>
 namespace b  = boost;
@@ -37,37 +37,14 @@
         typedef                     Filtration_                                         LSFiltration;
         typedef                     typename LSFiltration::Index                        LSFIndex;
 
-        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;
-        class                       KineticVertexType
-        {
-            public:
-                                        KineticVertexType(const Vertex& v):
-                                            vertex_(v)                                              {}
-
-                Key                     kinetic_key() const                                         { return key_; }
-                void                    set_kinetic_key(Key k)                                      { key_ = k; }
-
-                Vertex                  vertex() const                                              { return vertex_; }
-                void                    set_vertex(Vertex v)                                        { vertex_ = v; }
-
-                LSFIndex                simplex_index() const                                       { return simplex_index_; }
-                void                    set_simplex_index(LSFIndex i)                               { simplex_index_ = i; }
-                
-            private:
-                Key                     key_;
-                Vertex                  vertex_;
-                LSFIndex                simplex_index_;
-        };
+        typedef                     LinearKernel<VertexValue>                           KineticKernel;
+        typedef                     Simulator<KineticKernel>                            KineticSimulator;
+        class                       KineticVertexType;
         class                       KineticVertexComparison;
+        class                       TrajectoryExtractor;
         typedef                     typename OrderContainer<KineticVertexType>::Container        
                                                                                         VertexContainer;
         typedef                     typename VertexContainer::iterator                  VertexIndex;
-        typedef                     std::map<Key, VertexIndex>                          KeyVertexMap;
 
         struct                      AttachmentData: public VineData                     
         {   
@@ -103,7 +80,7 @@
                                                const VertexEvaluator& veval = VertexEvaluator());
                                     ~LSVineyard();
 
-        void                        compute_vineyard(const VertexEvaluator& veval, bool explicit_events = false);
+        void                        compute_vineyard(const VertexEvaluator& veval);
         bool                        transpose_vertices(VertexIndex vi);
         
         const LSFiltration&         filtration() const                                  { return filtration_; }
@@ -123,7 +100,7 @@
 
     public:
         // For Kinetic Sort
-        void                        swap(Key a, Key b);
+        void                        swap(VertexIndex a, KineticSimulator* simulator);
         
     private:
         void                        change_evaluator(Evaluator* eval);
@@ -151,8 +128,6 @@
         Vnrd                        vineyard_;
         Evaluator*                  evaluator_;
         unsigned                    time_count_;
-
-        KeyVertexMap                kinetic_map_;
                 
 #if 0
     private:
@@ -172,34 +147,48 @@
 };
 
 //BOOST_CLASS_EXPORT(LSVineyard)
-        
-// template<class V, class VE, class S, class C>
-// class LSVineyard<V,VE,S,C>::KineticVertexType
-// {
-//     public:
-//                                 KineticVertexType(const Vertex& v):
-//                                     vertex_(v)                                              {}
-
-//         Key                     kinetic_key() const                                         { return key_; }
-//         void                    set_kinetic_key(Key k)                                      { key_ = k; }
-
-//         Vertex                  vertex() const                                              { return vertex_; }
-//         void                    set_vertex(Vertex v)                                        { vertex_ = v; }
-
-//         LSFIndex                simplex_index() const                                       { return simplex_index_; }
-//         void                    set_simplex_index(iterator i)                               { simplex_index_ = i; }
-        
-//     private:
-//         Key                     key_;
-//         Vertex                  vertex_;
-//         LSFIndex                simplex_index_;
-// };
 
 template<class V, class VE, class S, class C>
 std::ostream& 
 operator<<(std::ostream& out, const typename LSVineyard<V,VE,S,C>::VertexIndex& vi)    
 { return out << vi->vertex(); }
+
+
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::KineticVertexType
+{
+    public:
+                                KineticVertexType(const Vertex& v):
+                                    vertex_(v)                                              {}
+
+        Vertex                  vertex() const                                              { return vertex_; }
+        void                    set_vertex(Vertex v)                                        { vertex_ = v; }
+
+        LSFIndex                simplex_index() const                                       { return simplex_index_; }
+        void                    set_simplex_index(LSFIndex i)                               { simplex_index_ = i; }
         
+    private:
+        Vertex                  vertex_;
+        LSFIndex                simplex_index_;
+};
+
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::TrajectoryExtractor: public std::unary_function<VertexIndex, typename KineticSimulator::Function>
+{
+    public:
+        typedef                 typename KineticSimulator::Function                         Function;
+
+                                TrajectoryExtractor(const VertexEvaluator& veval0, 
+                                                    const VertexEvaluator& veval1):
+                                    veval0_(veval0), veval1_(veval1)                        {}
+                                
+
+        Function                operator()(VertexIndex i) const                             { VertexValue v0 = veval0_(i->vertex()), v1 = veval1_(i->vertex()); return Function(v0, v1 - v0); }
+    
+    private:
+        const VertexEvaluator&  veval0_, veval1_;
+};
+
 template<class V, class VE, class S, class C>
 class LSVineyard<V,VE,S,C>::KineticVertexComparison: public std::binary_function<const KineticVertexType&, const KineticVertexType&, bool>
 {
@@ -245,20 +234,6 @@
 };
 
 template<class V, class VE, class S, class C>
-class LSVineyard<V,VE,S,C>::SortVisitor: public CGAL::Kinetic::Sort_visitor_base
-{
-    public:
-                                SortVisitor(LSVineyard& v): 
-                                    vineyard_(v)                                            {}
-
-        template<class Vertex_handle>
-        void                    pre_swap(Vertex_handle a, Vertex_handle b) const            { vineyard_.swap(*a,*b); }
-
-    private:
-        LSVineyard&             vineyard_;
-};
-
-template<class V, class VE, class S, class C>
 class LSVineyard<V,VE,S,C>::DimensionFromIterator: std::unary_function<iterator, Dimension>
 {
     public:
--- a/include/topology/lsvineyard.hpp	Fri Feb 05 11:49:52 2010 -0800
+++ b/include/topology/lsvineyard.hpp	Sat Feb 06 23:22:43 2010 -0800
@@ -111,50 +111,30 @@
 template<class V, class VE, class S, class F_>
 void                    
 LSVineyard<V,VE,S,F_>::
-compute_vineyard(const VertexEvaluator& veval, bool explicit_events)
+compute_vineyard(const VertexEvaluator& veval)
 {
-    typedef Traits::Kinetic_kernel::Point_1                                 Point;
-    typedef Traits::Kinetic_kernel::Function_kernel::Construct_function     CF; 
-    typedef Traits::Kinetic_kernel::Motion_function                         F; 
-    
-    Traits tr(0,1);
-    Simulator::Handle sp = tr.simulator_handle();
-    ActivePointsTable::Handle apt = tr.active_points_1_table_handle();
-    Sort sort(tr, SortVisitor(*this));
+    typedef     KineticSort<VertexIndex, TrajectoryExtractor, KineticSimulator>       KineticSortDS;
     
     // Setup the (linear) trajectories
     rLog(rlLSVineyard, "Setting up trajectories");
-    CF cf; 
-    kinetic_map_.clear();
+    KineticSimulator    simulator;
+    TrajectoryExtractor traj(veval_, veval);
     
-    // The reverse order takes advantage of how the kinetic sort arranges its elements 
-    // (inserts before the upper bound) to avoid problems with duplicates; it's unfortunate 
-    // that such dependence is necessary, it would be nice to eventually get rid of it.
-    for (VertexIndex cur = boost::prior(vertices_.end()); cur != boost::prior(vertices_.begin()); --cur)
-    {
-        VertexValue val0 = veval_(cur->vertex());
-        VertexValue val1 = veval(cur->vertex());
-        rLog(rlLSVineyardDebug, "Vertex %d: %f -> %f", cur->vertex(), val0, val1);
-        F x = cf(F::NT(val0), F::NT(val1 - val0));          // x = val0 + (val1 - val0)*t = (1-t)*val0 + t*val1
-        Point p(x);
-        vertices_.modify(cur,   b::bind(&KineticVertexType::set_kinetic_key, bl::_1, apt->insert(p)));    // cur->set_kinetic_key(apt->insert(p))
-        kinetic_map_[cur->kinetic_key()] = cur;
-        rLog(rlLSVineyardDebug, "Scheduling: %d %s", cur->vertex(), tostring(x).c_str());
-    }
+    KineticSortDS       sort(vertices_.begin(), vertices_.end(), 
+                             boost::bind(&LSVineyard::swap, this, bl::_1, bl::_2),
+                             &simulator, traj);
+    AssertMsg(sort.audit(&simulator), "Sort audit should succeed");
     
     // Process all the events (compute the vineyard in the process)
-    change_evaluator(new KineticEvaluator(*this, sp, apt, time_count_));
-    if (explicit_events)
+    change_evaluator(new KineticEvaluator(*this, simulator, time_count_, traj));
+    while (!simulator.reached_infinity() && simulator.next_event_time() < 1)
     {
-        while (sp->next_event_time() < 1)
-        {
-            rLog(rlLSVineyardDebug, "Next event time: %f", sp->next_event_time());
-            sp->set_current_event_number(sp->current_event_number() + 1);
-            rLog(rlLSVineyardDebug, "Processed event");
-        }
-    } else
-        sp->set_current_time(1.0);
-    rLog(rlLSVineyard, "Processed %d events", sp->current_event_number());
+        rLog(rlLSVineyardDebug, "Next event time: %f", simulator.next_event_time());
+        simulator.process();
+        rLog(rlLSVineyardDebug, "Processed event");
+    }
+    rLog(rlLSVineyard, "Processed %d events", simulator.event_count());
+    AssertMsg(sort.audit(&simulator), "Sort audit should succeed");
     
     veval_ = veval;
     change_evaluator(new StaticEvaluator(*this, ++time_count_));
@@ -164,14 +144,15 @@
 template<class V, class VE, class S, class F>
 void                    
 LSVineyard<V,VE,S,F>::
-swap(Key a, Key b)
+swap(VertexIndex a, KineticSimulator* simulator)
 {
+    VertexIndex b = boost::next(a);
     rLog(rlLSVineyardDebug, "Entered swap");
-    VertexIndex ao = kinetic_map_[a], bo = kinetic_map_[b];
-    rLog(rlLSVineyardDebug, "Vertices: %d %d compare %d", ao->vertex(), bo->vertex(), vcmp_(ao->vertex(), bo->vertex()));
-    AssertMsg(!vcmp_(bo->vertex(), ao->vertex()), "In swap(a,b), a must precede b");
-    transpose_vertices(ao);
-    // AssertMsg(vcmp_(bo->vertex(), ao->vertex()), "In swap(a,b), b must precede a after the transposition");
+    rLog(rlLSVineyardDebug, "Vertices: %d %d compare %d", a->vertex(), b->vertex(), vcmp_(a->vertex(), b->vertex()));
+    AssertMsg(!vcmp_(b->vertex(), a->vertex()), "In swap(a,b), a must precede b");
+    AssertMsg(a < b, "In swap(a,b), a must precede b");
+    transpose_vertices(a);
+    AssertMsg(b < a, "In swap(a,b), b must precede a after the transposition");
 }
 
 template<class V, class VE, class S, class F>
@@ -318,27 +299,30 @@
 class LSVineyard<V,VE,S,C>::KineticEvaluator: public Evaluator
 {
     public:
-                                KineticEvaluator(const LSVineyard& v, Simulator::Handle sp, ActivePointsTable::Handle apt, RealType time_offset): 
-                                    vineyard_(v), sp_(sp), apt_(apt), time_offset_(time_offset)           {}
+        typedef                 typename KineticSimulator::Time                             Time;
 
-        virtual RealType        time() const                                                { return time_offset_ + CGAL::to_double(get_time()); }
+                                KineticEvaluator(const LSVineyard& v, const KineticSimulator& sp, RealType time_offset, const TrajectoryExtractor& traj): 
+                                    vineyard_(v), sp_(sp), 
+                                    time_offset_(time_offset), traj_(traj)                  {}
+
+        virtual RealType        time() const                                                { return time_offset_ + get_time(); }
         virtual RealType        operator()(Index i) const                                   
         {
             rLog(rlLSVineyard, "%s (attached to %d): %s(%f) = %f", tostring(vineyard_.pfmap(i)).c_str(),
                                                                    i->attachment->vertex(),
-                                                                   tostring(apt_->at(i->attachment->kinetic_key()).x()).c_str(),
+                                                                   tostring(traj_(i->attachment)).c_str(),
                                                                    get_time(),
-                                                                   CGAL::to_double(apt_->at(i->attachment->kinetic_key()).x()(get_time())));
-            return CGAL::to_double(apt_->at(i->attachment->kinetic_key()).x()(get_time())); 
+                                                                   traj_(i->attachment)(get_time()));
+            return traj_(i->attachment)(get_time()); 
         }
         virtual Dimension       dimension(Index i) const                                    { return vineyard_.pfmap(i).dimension(); }
 
     private:
-        Simulator::Time         get_time() const                                            { return sp_->current_time(); }
+        Time                    get_time() const                                            { return sp_.current_time(); }
         
         const LSVineyard&           vineyard_;
-        Simulator::Handle           sp_;
-        ActivePointsTable::Handle   apt_;
+        const KineticSimulator&     sp_;
+        const TrajectoryExtractor&  traj_;
         RealType                    time_offset_;
 };