Switched to internal KineticSort dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Sat Feb 06 23:22:43 2010 -0800 (2010-02-06)
branchdev
changeset 1992bde4c56101c
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
     1.1 --- a/examples/pl-functions/combustion-vineyard.cpp	Fri Feb 05 11:49:52 2010 -0800
     1.2 +++ b/examples/pl-functions/combustion-vineyard.cpp	Sat Feb 06 23:22:43 2010 -0800
     1.3 @@ -81,7 +81,7 @@
     1.4      std::cout << "Pairing computed" << std::endl;
     1.5      
     1.6      // Compute vineyard
     1.7 -    v.compute_vineyard(g1, true);
     1.8 +    v.compute_vineyard(g1);
     1.9      std::cout << "Vineyard computed" << std::endl;
    1.10  
    1.11      v.vineyard().save_edges("combustion");
     2.1 --- a/examples/pl-functions/pl-vineyard.cpp	Fri Feb 05 11:49:52 2010 -0800
     2.2 +++ b/examples/pl-functions/pl-vineyard.cpp	Sat Feb 06 23:22:43 2010 -0800
     2.3 @@ -15,7 +15,7 @@
     2.4  namespace bl = boost::lambda;
     2.5  
     2.6  
     2.7 -typedef     float                                       VertexValue;
     2.8 +typedef     double                                      VertexValue;
     2.9  typedef     unsigned                                    Vertex;
    2.10  typedef     std::vector<VertexValue>                    VertexVector;
    2.11  struct SubscriptFunctor: public std::unary_function<Vertex, VertexValue>
    2.12 @@ -74,7 +74,7 @@
    2.13      for (size_t i = 1; i < vertices.size(); ++i)
    2.14      {
    2.15          veval = VertexEvaluator(vertices[i]);
    2.16 -        v.compute_vineyard(veval, explicit_events);
    2.17 +        v.compute_vineyard(veval);
    2.18          std::cout << "Processed frame: " << i << std::endl;
    2.19      }
    2.20      std::cout << "Vineyard computed" << std::endl;
     3.1 --- a/examples/pl-functions/test-grid2D-vineyard.cpp	Fri Feb 05 11:49:52 2010 -0800
     3.2 +++ b/examples/pl-functions/test-grid2D-vineyard.cpp	Sat Feb 06 23:22:43 2010 -0800
     3.3 @@ -42,7 +42,7 @@
     3.4          std::cout << "  " << v.filtration().simplex(cur) << std::endl;
     3.5  
     3.6      // Compute vineyard
     3.7 -    v.compute_vineyard(g1, true);
     3.8 +    v.compute_vineyard(g1);
     3.9      std::cout << "Vineyard computed" << std::endl;
    3.10      
    3.11      // Simplex order after
     4.1 --- a/include/geometry/kinetic-sort.h	Fri Feb 05 11:49:52 2010 -0800
     4.2 +++ b/include/geometry/kinetic-sort.h	Sat Feb 06 23:22:43 2010 -0800
     4.3 @@ -53,7 +53,7 @@
     4.4  		/// \name Core Functionality
     4.5  		/// @{
     4.6  									KineticSort();
     4.7 -									KineticSort(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator);
     4.8 +									KineticSort(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator, const TrajectoryExtractor& te);
     4.9  		void						initialize(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator);
    4.10  
    4.11  		void						insert(iterator pos, ElementIterator f, ElementIterator l, Simulator* simulator);
    4.12 @@ -68,6 +68,8 @@
    4.13  		iterator					begin() 									{ return list_.begin(); }
    4.14  		iterator					end() 										{ return list_.end(); }
    4.15  
    4.16 +        const TrajectoryExtractor&  trajectory_extractor() const                { return te_; }
    4.17 +
    4.18  	private:
    4.19  		class SwapEvent;
    4.20  		void						schedule_swaps(iterator b, iterator e, Simulator* s);
    4.21 @@ -76,6 +78,7 @@
    4.22  	private:
    4.23  		NodeList					list_;
    4.24  		Swap						swap_;	
    4.25 +        TrajectoryExtractor         te_;
    4.26  };
    4.27  
    4.28  #include "kinetic-sort.hpp"
     5.1 --- a/include/geometry/kinetic-sort.hpp	Fri Feb 05 11:49:52 2010 -0800
     5.2 +++ b/include/geometry/kinetic-sort.hpp	Sat Feb 06 23:22:43 2010 -0800
     5.3 @@ -21,7 +21,8 @@
     5.4  
     5.5  template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
     5.6  KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
     5.7 -KineticSort(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator)
     5.8 +KineticSort(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator, const TrajectoryExtractor& te = TrajectoryExtractor()):
     5.9 +    te_(te)
    5.10  {
    5.11  	initialize(b, e, swap, simulator);
    5.12  }
    5.13 @@ -129,16 +130,14 @@
    5.14  	Time t = simulator->audit_time();
    5.15  	rLog(rlKineticSortAudit, "Auditing at %s", tostring(t).c_str());
    5.16  
    5.17 -	TrajectoryExtractor	te;
    5.18 -	
    5.19  	typename NodeList::const_iterator next = list_.begin();
    5.20  	typename NodeList::const_iterator cur = next++;
    5.21 -	Function cur_trajectory = te(cur->element);
    5.22 +	Function cur_trajectory = te_(cur->element);
    5.23  	while (next != list_.end())
    5.24  	{
    5.25  		rLog(rlKineticSortAudit, "  %s", intostring(**(cur->swap_event_key)).c_str());
    5.26  
    5.27 -		Function next_trajectory = te(next->element);
    5.28 +		Function next_trajectory = te_(next->element);
    5.29  		rLog(rlKineticSortAudit, "  Auditing:   %s, %s", tostring(cur_trajectory).c_str(),
    5.30                                                           tostring(next_trajectory).c_str());
    5.31  		rLog(rlKineticSortAudit, "  Difference: %s", tostring(next_trajectory - cur_trajectory).c_str());
    5.32 @@ -146,8 +145,8 @@
    5.33                                                           tostring(FunctionKernel::sign_at(next_trajectory - cur_trajectory, t)).c_str());
    5.34  		if (FunctionKernel::sign_at(next_trajectory - cur_trajectory, t) == -1)
    5.35  		{
    5.36 -			rError("Audit failed at %s, %s", tostring(*cur->element).c_str(), 
    5.37 -                                             tostring(*next->element).c_str());
    5.38 +			// rError("Audit failed at %s, %s", tostring(*cur->element).c_str(), 
    5.39 +            //                                  tostring(*next->element).c_str());
    5.40  			return false;
    5.41  		}
    5.42  
    5.43 @@ -166,14 +165,12 @@
    5.44  {
    5.45  	typedef 		typename Simulator::Function		        Function;
    5.46  	
    5.47 -	TrajectoryExtractor	te;
    5.48 -	
    5.49  	iterator next = b; 
    5.50  	iterator cur = next++;
    5.51 -	Function cur_trajectory = te(cur->element);
    5.52 +	Function cur_trajectory = te_(cur->element);
    5.53  	while (next != e)
    5.54  	{
    5.55 -		Function next_trajectory = te(next->element);
    5.56 +		Function next_trajectory = te_(next->element);
    5.57  		rLog(rlKineticSortSchedule, "Next trajectory: %s", tostring(next_trajectory).c_str());
    5.58  		// TODO: add assertion that (next_trajectory - cur_trajectory)(s->curren_time()) > 0
    5.59  		cur->swap_event_key = simulator->add(next_trajectory - cur_trajectory, SwapEvent(this, cur));
    5.60 @@ -197,11 +194,9 @@
    5.61  		return;
    5.62  	}
    5.63  
    5.64 -	TrajectoryExtractor	te;
    5.65 -	
    5.66  	iterator next = boost::next(i); 
    5.67 -	Function i_trajectory = te(i->element);
    5.68 -	Function next_trajectory = te(next->element);
    5.69 +	Function i_trajectory = te_(i->element);
    5.70 +	Function next_trajectory = te_(next->element);
    5.71  	
    5.72  	//std::cout << "Updating swaps for: " << i_trajectory << ", " << next_trajectory << std::endl;
    5.73  	//std::cout << "Difference:         " << next_trajectory - i_trajectory << std::endl;
    5.74 @@ -247,6 +242,6 @@
    5.75  KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::SwapEvent::
    5.76  operator<<(std::ostream& out) const
    5.77  {
    5.78 -	Parent::operator<<(out) << "SwapEvent at " << TrajectoryExtractor_()(position()->element);
    5.79 +	Parent::operator<<(out) << "SwapEvent at " << sort_->trajectory_extractor()(position()->element);
    5.80  	return out;
    5.81  }
     6.1 --- a/include/geometry/simulator.h	Fri Feb 05 11:49:52 2010 -0800
     6.2 +++ b/include/geometry/simulator.h	Sat Feb 06 23:22:43 2010 -0800
     6.3 @@ -1,9 +1,10 @@
     6.4  #ifndef __SIMULATOR_H__
     6.5  #define __SIMULATOR_H__
     6.6  
     6.7 -#include "utilities/eventqueue.h"
     6.8 +#include <utilities/eventqueue.h>
     6.9 +#include <utilities/indirect.h>
    6.10  
    6.11 -template<class Comparison>  						class IndirectComparison;
    6.12 +#include <limits>
    6.13  
    6.14  /**
    6.15   * Simulator class. Keeps a queue of events. Infinity is reached if the Event 
    6.16 @@ -17,50 +18,54 @@
    6.17  template<class FuncKernel_, template<class Event> class EventComparison_ = std::less>
    6.18  class Simulator
    6.19  {
    6.20 -	public:
    6.21 -		typedef						FuncKernel_								    FunctionKernel;
    6.22 -		typedef						typename FunctionKernel::Function	        Function;
    6.23 -		typedef						typename FunctionKernel::RootStack		    RootStack;
    6.24 -		typedef						typename FunctionKernel::RootType			RootType;
    6.25 -		typedef						RootType									Time;
    6.26 +    public:
    6.27 +        typedef                     FuncKernel_                                 FunctionKernel;
    6.28 +        typedef                     typename FunctionKernel::Function           Function;
    6.29 +        typedef                     typename FunctionKernel::RootStack          RootStack;
    6.30 +        typedef                     typename FunctionKernel::RootType           RootType;
    6.31 +        typedef                     RootType                                    Time;
    6.32  
    6.33 -		class Event;
    6.34 -		typedef						EventComparison_<Event>						EventComparison;
    6.35 -		typedef						EventQueue<Event*, IndirectComparison<EventComparison> >			
    6.36 -																				EventQueueS;
    6.37 -		typedef						typename EventQueueS::iterator				Key;
    6.38 -		typedef						typename EventQueueS::const_iterator		const_Key;
    6.39 +        class Event;
    6.40 +        typedef                     EventComparison_<Event>                     EventComparison;
    6.41  
    6.42 +        class IndirectEventComparison;
    6.43 +        typedef                     EventQueue<Event*, IndirectEventComparison> EventQueueS;
    6.44 +        typedef                     typename EventQueueS::iterator              Key;
    6.45 +        typedef                     typename EventQueueS::const_iterator        const_Key;
    6.46  
    6.47 -									Simulator(Time start = FunctionKernel::root(0)):
    6.48 -										current_(start)         				{}
    6.49  
    6.50 +                                    Simulator(Time start = FunctionKernel::root(0)):
    6.51 +                                        current_(start)                         {}
    6.52  
    6.53 -		template<class Event_> 
    6.54 -		Key							add(const Event_& e);
    6.55 -		template<class Event_> 
    6.56 -		Key							add(const Function& f, const Event_& e);
    6.57 -		void						process();
    6.58 -		void						update(Key k, const Function& f);
    6.59 -		
    6.60 -		void						remove(Key k)								{ queue_.remove(k); }
    6.61 -		Key							null_key() 									{ return queue_.end(); }
    6.62  
    6.63 -		Time						current_time() const						{ return current_; }
    6.64 -		Time						audit_time() const;
    6.65 -		bool						reached_infinity() const					{ return queue_.empty() || (*queue_.top())->root_stack().empty(); }
    6.66 +        template<class Event_> 
    6.67 +        Key                         add(const Event_& e);
    6.68 +        template<class Event_> 
    6.69 +        Key                         add(const Function& f, const Event_& e);
    6.70 +        void                        process();
    6.71 +        void                        update(Key k, const Function& f);
    6.72 +        
    6.73 +        void                        remove(Key k)                               { queue_.remove(k); }
    6.74 +        Key                         null_key()                                  { return queue_.end(); }
    6.75 +
    6.76 +        Time                        current_time() const                        { return current_; }
    6.77 +        Time                        audit_time() const;
    6.78 +        bool                        reached_infinity() const                    { return queue_.empty() || (*queue_.top())->root_stack().empty(); }
    6.79 +        Time                        next_event_time() const                     { return queue_.empty() ? std::numeric_limits<Time>::infinity():(*queue_.top())->root_stack().top(); } 
    6.80          
    6.81          Event*                      top() const                                 { return *(queue_.top()); }
    6.82          unsigned                    size() const                                { return queue_.size(); }
    6.83 +        unsigned                    event_count() const                         { return count_; }
    6.84  
    6.85 -		std::ostream&				operator<<(std::ostream& out) const;
    6.86 +        std::ostream&               operator<<(std::ostream& out) const;
    6.87  
    6.88 -	private:
    6.89 -		void						update(Key i);
    6.90 +    private:
    6.91 +        void                        update(Key i);
    6.92  
    6.93 -	private:
    6.94 -		Time						current_;
    6.95 -		EventQueueS					queue_;
    6.96 +    private:
    6.97 +        Time                        current_;
    6.98 +        EventQueueS                 queue_;
    6.99 +        unsigned                    count_;
   6.100  };
   6.101  
   6.102  
   6.103 @@ -72,30 +77,30 @@
   6.104  template<class FuncKernel_, template<class Event> class EventComparison_>
   6.105  class Simulator<FuncKernel_, EventComparison_>::Event
   6.106  {
   6.107 -	public:
   6.108 -		typedef						FuncKernel_									FunctionKernel;
   6.109 -		typedef						typename FunctionKernel::RootStack		    RootStack;
   6.110 +    public:
   6.111 +        typedef                     FuncKernel_                                 FunctionKernel;
   6.112 +        typedef                     typename FunctionKernel::RootStack          RootStack;
   6.113  
   6.114          /// process() is called when the event is at the top of the queue 
   6.115          /// in the simulator.
   6.116 -		/// Returns true if the event needs to remain in the Simulator 
   6.117 -		/// (top of the root_stack() will be used for new time).
   6.118 -		virtual	bool				process(Simulator* s) const					=0;
   6.119 -		
   6.120 -		RootStack&					root_stack()								{ return root_stack_; }
   6.121 -		const RootStack&			root_stack() const							{ return root_stack_; }
   6.122 +        /// Returns true if the event needs to remain in the Simulator 
   6.123 +        /// (top of the root_stack() will be used for new time).
   6.124 +        virtual bool                process(Simulator* s) const                 =0;
   6.125 +        
   6.126 +        RootStack&                  root_stack()                                { return root_stack_; }
   6.127 +        const RootStack&            root_stack() const                          { return root_stack_; }
   6.128  
   6.129 -		bool						operator<(const Event& e) const				
   6.130 -		{ 
   6.131 -			if (root_stack().empty())
   6.132 -				return false;
   6.133 -			else if (e.root_stack().empty())
   6.134 -				return true;
   6.135 -			else
   6.136 -				return root_stack().top() < e.root_stack().top();
   6.137 -		}
   6.138 +        bool                        operator<(const Event& e) const             
   6.139 +        { 
   6.140 +            if (root_stack().empty())
   6.141 +                return false;
   6.142 +            else if (e.root_stack().empty())
   6.143 +                return true;
   6.144 +            else
   6.145 +                return root_stack().top() < e.root_stack().top();
   6.146 +        }
   6.147  
   6.148 -		virtual std::ostream&		operator<<(std::ostream& out) const			
   6.149 +        virtual std::ostream&       operator<<(std::ostream& out) const         
   6.150          { 
   6.151              out << "Event with " << root_stack_.size() << " roots"; 
   6.152              if (!root_stack_.empty()) out << "; top root: " << root_stack_.top();
   6.153 @@ -103,28 +108,30 @@
   6.154              return out;
   6.155          }
   6.156  
   6.157 -	private:
   6.158 -		RootStack					root_stack_;
   6.159 +    private:
   6.160 +        RootStack                   root_stack_;
   6.161  };
   6.162  
   6.163  /**
   6.164   * Compares elements pointed at by its arguments using the provided Comparison_ 
   6.165   * (which must not take any arguments during construction).
   6.166   */
   6.167 -template<class Comparison_>
   6.168 -class IndirectComparison: public std::binary_function<const typename Comparison_::first_argument_type*, 
   6.169 -													  const typename Comparison_::second_argument_type*, 
   6.170 -													  bool>
   6.171 +template<class FuncKernel_, template<class Event> class EventComparison_>
   6.172 +class Simulator<FuncKernel_, EventComparison_>::IndirectEventComparison: 
   6.173 +    public std::binary_function<const typename EventComparison::first_argument_type*, 
   6.174 +                                const typename EventComparison::second_argument_type*, 
   6.175 +                                bool>
   6.176  {
   6.177 -	public:
   6.178 -		typedef						Comparison_											Comparison;
   6.179 -		typedef						const typename Comparison::first_argument_type*		first_argument_type;
   6.180 -		typedef						const typename Comparison::second_argument_type*	second_argument_type;
   6.181 +    public:
   6.182 +        typedef                     EventComparison                                     Comparison;
   6.183 +        typedef                     const typename Comparison::first_argument_type*     first_argument_type;
   6.184 +        typedef                     const typename Comparison::second_argument_type*    second_argument_type;
   6.185  
   6.186 -		bool						operator()(first_argument_type e1, second_argument_type e2) const
   6.187 -		{ return Comparison()(*e1, *e2); }
   6.188 +        bool                        operator()(first_argument_type e1, second_argument_type e2) const
   6.189 +        { return Comparison()(*e1, *e2); }
   6.190  };
   6.191  
   6.192 +
   6.193  #include "simulator.hpp"
   6.194  
   6.195  #endif // __SIMULATOR_H__
     7.1 --- a/include/geometry/simulator.hpp	Fri Feb 05 11:49:52 2010 -0800
     7.2 +++ b/include/geometry/simulator.hpp	Sat Feb 06 23:22:43 2010 -0800
     7.3 @@ -88,10 +88,12 @@
     7.4  	current_ = e->root_stack().top(); e->root_stack().pop();
     7.5  	
     7.6      // Get the top element out of the queue, put it back depending on what process() says
     7.7 -    EventQueueS tmp; tmp.prepend(top, queue_);
     7.8 +    EventQueueS             tmp; tmp.prepend(top, queue_);
     7.9  
    7.10  	if (e->process(this))				{ queue_.prepend(top, tmp); update(top); }
    7.11  	else								{ delete e; }
    7.12 +
    7.13 +    ++count_;
    7.14  }
    7.15  
    7.16  template<class FuncKernel_, template<class Event> class EventComparison_>
     8.1 --- a/include/topology/lsvineyard.h	Fri Feb 05 11:49:52 2010 -0800
     8.2 +++ b/include/topology/lsvineyard.h	Sat Feb 06 23:22:43 2010 -0800
     8.3 @@ -1,6 +1,6 @@
     8.4  /**
     8.5   * Author: Dmitriy Morozov
     8.6 - * Department of Computer Science, Duke University, 2005 -- 2009
     8.7 + * Department of Computer Science, Duke University, 2005 -- 2010
     8.8   */
     8.9  
    8.10  #ifndef __LSVINEYARD_H__
    8.11 @@ -15,9 +15,9 @@
    8.12  
    8.13  #include <utilities/indirect.h>
    8.14  
    8.15 -#include <CGAL/Kinetic/Inexact_simulation_traits.h>
    8.16 -#include <CGAL/Kinetic/Sort.h>
    8.17 -#include <CGAL/Kinetic/Sort_visitor_base.h>
    8.18 +#include <geometry/simulator.h>
    8.19 +#include <geometry/kinetic-sort.h>
    8.20 +#include <geometry/linear-kernel.h>
    8.21  
    8.22  #include <boost/tuple/tuple.hpp>
    8.23  namespace b  = boost;
    8.24 @@ -37,37 +37,14 @@
    8.25          typedef                     Filtration_                                         LSFiltration;
    8.26          typedef                     typename LSFiltration::Index                        LSFIndex;
    8.27  
    8.28 -        class                       SortVisitor;
    8.29 -        typedef                     CGAL::Kinetic::Inexact_simulation_traits            Traits;
    8.30 -        typedef                     CGAL::Kinetic::Sort<Traits, SortVisitor>            Sort;
    8.31 -        typedef                     Traits::Simulator                                   Simulator;
    8.32 -        typedef                     Traits::Active_points_1_table                       ActivePointsTable;
    8.33 -        typedef                     ActivePointsTable::Key                              Key;
    8.34 -        class                       KineticVertexType
    8.35 -        {
    8.36 -            public:
    8.37 -                                        KineticVertexType(const Vertex& v):
    8.38 -                                            vertex_(v)                                              {}
    8.39 -
    8.40 -                Key                     kinetic_key() const                                         { return key_; }
    8.41 -                void                    set_kinetic_key(Key k)                                      { key_ = k; }
    8.42 -
    8.43 -                Vertex                  vertex() const                                              { return vertex_; }
    8.44 -                void                    set_vertex(Vertex v)                                        { vertex_ = v; }
    8.45 -
    8.46 -                LSFIndex                simplex_index() const                                       { return simplex_index_; }
    8.47 -                void                    set_simplex_index(LSFIndex i)                               { simplex_index_ = i; }
    8.48 -                
    8.49 -            private:
    8.50 -                Key                     key_;
    8.51 -                Vertex                  vertex_;
    8.52 -                LSFIndex                simplex_index_;
    8.53 -        };
    8.54 +        typedef                     LinearKernel<VertexValue>                           KineticKernel;
    8.55 +        typedef                     Simulator<KineticKernel>                            KineticSimulator;
    8.56 +        class                       KineticVertexType;
    8.57          class                       KineticVertexComparison;
    8.58 +        class                       TrajectoryExtractor;
    8.59          typedef                     typename OrderContainer<KineticVertexType>::Container        
    8.60                                                                                          VertexContainer;
    8.61          typedef                     typename VertexContainer::iterator                  VertexIndex;
    8.62 -        typedef                     std::map<Key, VertexIndex>                          KeyVertexMap;
    8.63  
    8.64          struct                      AttachmentData: public VineData                     
    8.65          {   
    8.66 @@ -103,7 +80,7 @@
    8.67                                                 const VertexEvaluator& veval = VertexEvaluator());
    8.68                                      ~LSVineyard();
    8.69  
    8.70 -        void                        compute_vineyard(const VertexEvaluator& veval, bool explicit_events = false);
    8.71 +        void                        compute_vineyard(const VertexEvaluator& veval);
    8.72          bool                        transpose_vertices(VertexIndex vi);
    8.73          
    8.74          const LSFiltration&         filtration() const                                  { return filtration_; }
    8.75 @@ -123,7 +100,7 @@
    8.76  
    8.77      public:
    8.78          // For Kinetic Sort
    8.79 -        void                        swap(Key a, Key b);
    8.80 +        void                        swap(VertexIndex a, KineticSimulator* simulator);
    8.81          
    8.82      private:
    8.83          void                        change_evaluator(Evaluator* eval);
    8.84 @@ -151,8 +128,6 @@
    8.85          Vnrd                        vineyard_;
    8.86          Evaluator*                  evaluator_;
    8.87          unsigned                    time_count_;
    8.88 -
    8.89 -        KeyVertexMap                kinetic_map_;
    8.90                  
    8.91  #if 0
    8.92      private:
    8.93 @@ -172,34 +147,48 @@
    8.94  };
    8.95  
    8.96  //BOOST_CLASS_EXPORT(LSVineyard)
    8.97 -        
    8.98 -// template<class V, class VE, class S, class C>
    8.99 -// class LSVineyard<V,VE,S,C>::KineticVertexType
   8.100 -// {
   8.101 -//     public:
   8.102 -//                                 KineticVertexType(const Vertex& v):
   8.103 -//                                     vertex_(v)                                              {}
   8.104 -
   8.105 -//         Key                     kinetic_key() const                                         { return key_; }
   8.106 -//         void                    set_kinetic_key(Key k)                                      { key_ = k; }
   8.107 -
   8.108 -//         Vertex                  vertex() const                                              { return vertex_; }
   8.109 -//         void                    set_vertex(Vertex v)                                        { vertex_ = v; }
   8.110 -
   8.111 -//         LSFIndex                simplex_index() const                                       { return simplex_index_; }
   8.112 -//         void                    set_simplex_index(iterator i)                               { simplex_index_ = i; }
   8.113 -        
   8.114 -//     private:
   8.115 -//         Key                     key_;
   8.116 -//         Vertex                  vertex_;
   8.117 -//         LSFIndex                simplex_index_;
   8.118 -// };
   8.119  
   8.120  template<class V, class VE, class S, class C>
   8.121  std::ostream& 
   8.122  operator<<(std::ostream& out, const typename LSVineyard<V,VE,S,C>::VertexIndex& vi)    
   8.123  { return out << vi->vertex(); }
   8.124 +
   8.125 +
   8.126 +template<class V, class VE, class S, class C>
   8.127 +class LSVineyard<V,VE,S,C>::KineticVertexType
   8.128 +{
   8.129 +    public:
   8.130 +                                KineticVertexType(const Vertex& v):
   8.131 +                                    vertex_(v)                                              {}
   8.132 +
   8.133 +        Vertex                  vertex() const                                              { return vertex_; }
   8.134 +        void                    set_vertex(Vertex v)                                        { vertex_ = v; }
   8.135 +
   8.136 +        LSFIndex                simplex_index() const                                       { return simplex_index_; }
   8.137 +        void                    set_simplex_index(LSFIndex i)                               { simplex_index_ = i; }
   8.138          
   8.139 +    private:
   8.140 +        Vertex                  vertex_;
   8.141 +        LSFIndex                simplex_index_;
   8.142 +};
   8.143 +
   8.144 +template<class V, class VE, class S, class C>
   8.145 +class LSVineyard<V,VE,S,C>::TrajectoryExtractor: public std::unary_function<VertexIndex, typename KineticSimulator::Function>
   8.146 +{
   8.147 +    public:
   8.148 +        typedef                 typename KineticSimulator::Function                         Function;
   8.149 +
   8.150 +                                TrajectoryExtractor(const VertexEvaluator& veval0, 
   8.151 +                                                    const VertexEvaluator& veval1):
   8.152 +                                    veval0_(veval0), veval1_(veval1)                        {}
   8.153 +                                
   8.154 +
   8.155 +        Function                operator()(VertexIndex i) const                             { VertexValue v0 = veval0_(i->vertex()), v1 = veval1_(i->vertex()); return Function(v0, v1 - v0); }
   8.156 +    
   8.157 +    private:
   8.158 +        const VertexEvaluator&  veval0_, veval1_;
   8.159 +};
   8.160 +
   8.161  template<class V, class VE, class S, class C>
   8.162  class LSVineyard<V,VE,S,C>::KineticVertexComparison: public std::binary_function<const KineticVertexType&, const KineticVertexType&, bool>
   8.163  {
   8.164 @@ -245,20 +234,6 @@
   8.165  };
   8.166  
   8.167  template<class V, class VE, class S, class C>
   8.168 -class LSVineyard<V,VE,S,C>::SortVisitor: public CGAL::Kinetic::Sort_visitor_base
   8.169 -{
   8.170 -    public:
   8.171 -                                SortVisitor(LSVineyard& v): 
   8.172 -                                    vineyard_(v)                                            {}
   8.173 -
   8.174 -        template<class Vertex_handle>
   8.175 -        void                    pre_swap(Vertex_handle a, Vertex_handle b) const            { vineyard_.swap(*a,*b); }
   8.176 -
   8.177 -    private:
   8.178 -        LSVineyard&             vineyard_;
   8.179 -};
   8.180 -
   8.181 -template<class V, class VE, class S, class C>
   8.182  class LSVineyard<V,VE,S,C>::DimensionFromIterator: std::unary_function<iterator, Dimension>
   8.183  {
   8.184      public:
     9.1 --- a/include/topology/lsvineyard.hpp	Fri Feb 05 11:49:52 2010 -0800
     9.2 +++ b/include/topology/lsvineyard.hpp	Sat Feb 06 23:22:43 2010 -0800
     9.3 @@ -111,50 +111,30 @@
     9.4  template<class V, class VE, class S, class F_>
     9.5  void                    
     9.6  LSVineyard<V,VE,S,F_>::
     9.7 -compute_vineyard(const VertexEvaluator& veval, bool explicit_events)
     9.8 +compute_vineyard(const VertexEvaluator& veval)
     9.9  {
    9.10 -    typedef Traits::Kinetic_kernel::Point_1                                 Point;
    9.11 -    typedef Traits::Kinetic_kernel::Function_kernel::Construct_function     CF; 
    9.12 -    typedef Traits::Kinetic_kernel::Motion_function                         F; 
    9.13 -    
    9.14 -    Traits tr(0,1);
    9.15 -    Simulator::Handle sp = tr.simulator_handle();
    9.16 -    ActivePointsTable::Handle apt = tr.active_points_1_table_handle();
    9.17 -    Sort sort(tr, SortVisitor(*this));
    9.18 +    typedef     KineticSort<VertexIndex, TrajectoryExtractor, KineticSimulator>       KineticSortDS;
    9.19      
    9.20      // Setup the (linear) trajectories
    9.21      rLog(rlLSVineyard, "Setting up trajectories");
    9.22 -    CF cf; 
    9.23 -    kinetic_map_.clear();
    9.24 +    KineticSimulator    simulator;
    9.25 +    TrajectoryExtractor traj(veval_, veval);
    9.26      
    9.27 -    // The reverse order takes advantage of how the kinetic sort arranges its elements 
    9.28 -    // (inserts before the upper bound) to avoid problems with duplicates; it's unfortunate 
    9.29 -    // that such dependence is necessary, it would be nice to eventually get rid of it.
    9.30 -    for (VertexIndex cur = boost::prior(vertices_.end()); cur != boost::prior(vertices_.begin()); --cur)
    9.31 -    {
    9.32 -        VertexValue val0 = veval_(cur->vertex());
    9.33 -        VertexValue val1 = veval(cur->vertex());
    9.34 -        rLog(rlLSVineyardDebug, "Vertex %d: %f -> %f", cur->vertex(), val0, val1);
    9.35 -        F x = cf(F::NT(val0), F::NT(val1 - val0));          // x = val0 + (val1 - val0)*t = (1-t)*val0 + t*val1
    9.36 -        Point p(x);
    9.37 -        vertices_.modify(cur,   b::bind(&KineticVertexType::set_kinetic_key, bl::_1, apt->insert(p)));    // cur->set_kinetic_key(apt->insert(p))
    9.38 -        kinetic_map_[cur->kinetic_key()] = cur;
    9.39 -        rLog(rlLSVineyardDebug, "Scheduling: %d %s", cur->vertex(), tostring(x).c_str());
    9.40 -    }
    9.41 +    KineticSortDS       sort(vertices_.begin(), vertices_.end(), 
    9.42 +                             boost::bind(&LSVineyard::swap, this, bl::_1, bl::_2),
    9.43 +                             &simulator, traj);
    9.44 +    AssertMsg(sort.audit(&simulator), "Sort audit should succeed");
    9.45      
    9.46      // Process all the events (compute the vineyard in the process)
    9.47 -    change_evaluator(new KineticEvaluator(*this, sp, apt, time_count_));
    9.48 -    if (explicit_events)
    9.49 +    change_evaluator(new KineticEvaluator(*this, simulator, time_count_, traj));
    9.50 +    while (!simulator.reached_infinity() && simulator.next_event_time() < 1)
    9.51      {
    9.52 -        while (sp->next_event_time() < 1)
    9.53 -        {
    9.54 -            rLog(rlLSVineyardDebug, "Next event time: %f", sp->next_event_time());
    9.55 -            sp->set_current_event_number(sp->current_event_number() + 1);
    9.56 -            rLog(rlLSVineyardDebug, "Processed event");
    9.57 -        }
    9.58 -    } else
    9.59 -        sp->set_current_time(1.0);
    9.60 -    rLog(rlLSVineyard, "Processed %d events", sp->current_event_number());
    9.61 +        rLog(rlLSVineyardDebug, "Next event time: %f", simulator.next_event_time());
    9.62 +        simulator.process();
    9.63 +        rLog(rlLSVineyardDebug, "Processed event");
    9.64 +    }
    9.65 +    rLog(rlLSVineyard, "Processed %d events", simulator.event_count());
    9.66 +    AssertMsg(sort.audit(&simulator), "Sort audit should succeed");
    9.67      
    9.68      veval_ = veval;
    9.69      change_evaluator(new StaticEvaluator(*this, ++time_count_));
    9.70 @@ -164,14 +144,15 @@
    9.71  template<class V, class VE, class S, class F>
    9.72  void                    
    9.73  LSVineyard<V,VE,S,F>::
    9.74 -swap(Key a, Key b)
    9.75 +swap(VertexIndex a, KineticSimulator* simulator)
    9.76  {
    9.77 +    VertexIndex b = boost::next(a);
    9.78      rLog(rlLSVineyardDebug, "Entered swap");
    9.79 -    VertexIndex ao = kinetic_map_[a], bo = kinetic_map_[b];
    9.80 -    rLog(rlLSVineyardDebug, "Vertices: %d %d compare %d", ao->vertex(), bo->vertex(), vcmp_(ao->vertex(), bo->vertex()));
    9.81 -    AssertMsg(!vcmp_(bo->vertex(), ao->vertex()), "In swap(a,b), a must precede b");
    9.82 -    transpose_vertices(ao);
    9.83 -    // AssertMsg(vcmp_(bo->vertex(), ao->vertex()), "In swap(a,b), b must precede a after the transposition");
    9.84 +    rLog(rlLSVineyardDebug, "Vertices: %d %d compare %d", a->vertex(), b->vertex(), vcmp_(a->vertex(), b->vertex()));
    9.85 +    AssertMsg(!vcmp_(b->vertex(), a->vertex()), "In swap(a,b), a must precede b");
    9.86 +    AssertMsg(a < b, "In swap(a,b), a must precede b");
    9.87 +    transpose_vertices(a);
    9.88 +    AssertMsg(b < a, "In swap(a,b), b must precede a after the transposition");
    9.89  }
    9.90  
    9.91  template<class V, class VE, class S, class F>
    9.92 @@ -318,27 +299,30 @@
    9.93  class LSVineyard<V,VE,S,C>::KineticEvaluator: public Evaluator
    9.94  {
    9.95      public:
    9.96 -                                KineticEvaluator(const LSVineyard& v, Simulator::Handle sp, ActivePointsTable::Handle apt, RealType time_offset): 
    9.97 -                                    vineyard_(v), sp_(sp), apt_(apt), time_offset_(time_offset)           {}
    9.98 +        typedef                 typename KineticSimulator::Time                             Time;
    9.99  
   9.100 -        virtual RealType        time() const                                                { return time_offset_ + CGAL::to_double(get_time()); }
   9.101 +                                KineticEvaluator(const LSVineyard& v, const KineticSimulator& sp, RealType time_offset, const TrajectoryExtractor& traj): 
   9.102 +                                    vineyard_(v), sp_(sp), 
   9.103 +                                    time_offset_(time_offset), traj_(traj)                  {}
   9.104 +
   9.105 +        virtual RealType        time() const                                                { return time_offset_ + get_time(); }
   9.106          virtual RealType        operator()(Index i) const                                   
   9.107          {
   9.108              rLog(rlLSVineyard, "%s (attached to %d): %s(%f) = %f", tostring(vineyard_.pfmap(i)).c_str(),
   9.109                                                                     i->attachment->vertex(),
   9.110 -                                                                   tostring(apt_->at(i->attachment->kinetic_key()).x()).c_str(),
   9.111 +                                                                   tostring(traj_(i->attachment)).c_str(),
   9.112                                                                     get_time(),
   9.113 -                                                                   CGAL::to_double(apt_->at(i->attachment->kinetic_key()).x()(get_time())));
   9.114 -            return CGAL::to_double(apt_->at(i->attachment->kinetic_key()).x()(get_time())); 
   9.115 +                                                                   traj_(i->attachment)(get_time()));
   9.116 +            return traj_(i->attachment)(get_time()); 
   9.117          }
   9.118          virtual Dimension       dimension(Index i) const                                    { return vineyard_.pfmap(i).dimension(); }
   9.119  
   9.120      private:
   9.121 -        Simulator::Time         get_time() const                                            { return sp_->current_time(); }
   9.122 +        Time                    get_time() const                                            { return sp_.current_time(); }
   9.123          
   9.124          const LSVineyard&           vineyard_;
   9.125 -        Simulator::Handle           sp_;
   9.126 -        ActivePointsTable::Handle   apt_;
   9.127 +        const KineticSimulator&     sp_;
   9.128 +        const TrajectoryExtractor&  traj_;
   9.129          RealType                    time_offset_;
   9.130  };
   9.131