--- a/CMakeLists.txt Sun Nov 25 09:12:36 2007 -0500
+++ b/CMakeLists.txt Tue Feb 19 06:56:49 2008 -0500
@@ -40,6 +40,13 @@
${gmpxx_LIBRARY}
${m_LIBRARY})
+# SYNAPS
+add_definitions (-DBOOST_UBLAS_TYPE_CHECK=0)
+find_library (synaps_LIBRARY NAMES synaps)
+set (synaps_libraries ${synaps_LIBRARY}
+ ${gmp_LIBRARY}
+ ${gmpxx_LIBRARY})
+
# Debugging
if (debug)
if (optimize)
--- a/examples/ar-vineyard/CMakeLists.txt Sun Nov 25 09:12:36 2007 -0500
+++ b/examples/ar-vineyard/CMakeLists.txt Tue Feb 19 06:56:49 2008 -0500
@@ -3,6 +3,6 @@
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} ${synaps_libraries})
endforeach (t)
--- a/examples/ar-vineyard/ar-simplex3d.h Sun Nov 25 09:12:36 2007 -0500
+++ b/examples/ar-vineyard/ar-simplex3d.h Tue Feb 19 06:56:49 2008 -0500
@@ -86,7 +86,7 @@
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_;
};
--- a/examples/ar-vineyard/ar-vineyard.h Sun Nov 25 09:12:36 2007 -0500
+++ b/examples/ar-vineyard/ar-vineyard.h Tue Feb 19 06:56:49 2008 -0500
@@ -8,68 +8,80 @@
#include "topology/conesimplex.h"
#include "topology/filtration.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 "geometry/kinetic-sort.h"
#include <list>
#include "ar-simplex3d.h"
-#include <vector>
-
-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>
+class ARConeSimplex: public ConeSimplex<ARSimplex3D>
{
public:
typedef ConeSimplex<ARSimplex3D> Parent;
typedef ARSimplex3D ARSimplex3D;
+ typedef Filtration<ARConeSimplex> Filtration;
+ typedef Filtration::Index Index;
+
+ /// \name Polynomial Kernel types
+ /// @{
+ typedef double FieldType;
+ typedef UPolynomial<FieldType> PolyKernel;
+ typedef PolyKernel::Polynomial Polynomial;
+ typedef Simulator<PolyKernel> Simulator;
+
+ typedef KineticSort<Index, SimplexTrajectoryExtractor, Simulator>
+ SimplexSort;
+ typedef SimplexSort::iterator SimplexSortIterator;
+ typedef SimplexSortIterator Key;
+ /// @}
+
+ /// \name Kinetic Sort types
+ /// @{
+ typedef std::list<Polynomial> ThresholdList;
+
+ struct ThresholdTrajectoryExtractor
+ { Polynomial operator()(ThresholdList::iterator i) const { return *i; } }
+ struct SimplexTrajectoryExtractor
+ { Polynomial operator()(Index i) const { i->thresholds().front(); }
+
+ typedef KineticSort<ThresholdList::iterator, ThresholdTrajectoryExtractor, Simulator>
+ ThresholdSort;
+ /// @}
ARConeSimplex(const ARSimplex3D& s, bool coned = false):
Parent(s, coned) {}
Key kinetic_key() const { return key_; }
void set_kinetic_key(Key k) { key_ = k; }
+ ThresholdList& thresholds() { return thresholds_; }
+
+ void schedule_thresholds(SimplexSort* sort, Simulator* simulator);
+
private:
Key key_;
+ ThresholdList thresholds_;
+ ThresholdSort thresholds_sort_;
+
+ void swap_thresholds(SimplexSort* sort, ThresholdList::iterator i, Simulator* simulator);
};
-class ARVineyard: public ARVineyardBase
+class ARVineyard
{
public:
typedef ARVineyard Self;
- typedef Filtration<ARConeSimplex> ARFiltration;
+ typedef ARConeSimplex::Filtration ARFiltration;
typedef ARFiltration::Simplex Simplex;
typedef ARFiltration::Index Index;
typedef ARFiltration::Vineyard Vineyard;
typedef Vineyard::Evaluator Evaluator;
- typedef std::map<Key, Index> KeyIndexMap;
+ typedef ARConeSimplex::Simulator Simulator;
+ typedef ARConeSimplex::SimplexSort SimplexSort;
+
+
typedef std::list<Point> PointList;
class StaticEvaluator;
@@ -86,8 +98,7 @@
const Vineyard* vineyard() const { return vineyard_; }
public:
- // For Kinetic Sort
- void swap(Key a, Key b);
+ void swap(Index i, Simulator* simulator); ///< For kinetic sort
private:
void add_simplices();
@@ -98,8 +109,6 @@
Vineyard* vineyard_;
Evaluator* evaluator_;
- KeyIndexMap kinetic_map_;
-
Point z_;
Delaunay dt_;
@@ -121,26 +130,6 @@
//BOOST_CLASS_EXPORT(ARVineyard)
-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_;
-};
-
-std::ostream& operator<<(std::ostream& out, const ARVineyardBase::MembershipFunctionChangeEvent& e)
-{ return e.operator<<(out); }
-
class ARVineyard::StaticEvaluator: public Evaluator
{
public:
@@ -172,19 +161,6 @@
};
-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_;
-};
-
-
#include "ar-vineyard.hpp"
#endif // __AR_VINEYARD_H__
--- a/examples/ar-vineyard/ar-vineyard.hpp Sun Nov 25 09:12:36 2007 -0500
+++ b/examples/ar-vineyard/ar-vineyard.hpp Tue Feb 19 06:56:49 2008 -0500
@@ -2,8 +2,27 @@
/* Implementation */
+void
+ARConeSimplex::
+swap_thresholds(SimplexSort* sort, ThresholdList::iterator i, Simulator* simulator)
+{
+ typename ThresholdList::iterator n = boost::next(i);
+ tl->splice(i, *tl, n);
+ if (n == tl->begin())
+ sort->update_trajectory(kinetic_key(), simulator);
+}
+
+void
+ARConeSimplex::
+schedule_thresholds(SimplexSort* sort, Simulator* simulator)
+{
+ thresholds_sort_.initialize(thresholds_.begin(), thresholds_.end(),
+ boost::bind(&ARConeSimplex::swap_thresholds, this, sort, _1, _2), simulator);
+}
+
+
ARVineyard::
-ARVineyard(const PointList& points, const Point& z): z_(z)
+ARVineyard(const PointList& points, const Point& z): simplex_sort_(0), z_(z)
{
for (PointList::const_iterator cur = points.begin(); cur != points.end(); ++cur)
dt_.insert(*cur);
@@ -19,7 +38,7 @@
filtration_ = new ARFiltration(vineyard_);
for (ARSimplex3DVector::const_iterator cur = alpha_ordering.begin(); cur != alpha_ordering.end(); ++cur)
{
- filtration_->append(*cur); // Delaunay simplex
+ filtration_->append(ARConeSimplex(*cur)); // Delaunay simplex
filtration_->append(ARConeSimplex(*cur, true)); // Coned off delaunay simplex
}
}
@@ -48,24 +67,13 @@
{
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;
+ Simulator 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();
+ // Set 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;
+ cur->thresholds.push_back(ARConeSimplex::Polynomial(CGAL::to_double(cur->alpha())));
if (!cur->coned()) continue; // non-coned simplices stay put, so we are done
@@ -101,7 +109,16 @@
std::cout << "phi_lambda: " << phi_lambda << std::endl;
std::cout << "s^2 + v^2: " << sv << std::endl;
std::cout << std::endl;
+
+ cur->set_kinetic_key();
}
+
+
+ // Once thresholds are set (and sorted), we can initialize the simplex_sort
+ simplex_sort.initialize(filtration_.begin(), filtration_.end(),
+ boost::bind(&ARVineyard::swap, this, _1, _2), &simulator);
+
+
// Process all the events (compute the vineyard in the process)
// FIXME: the time should not be 1, but something like twice the radius of
--- a/include/geometry/kinetic-sort.h Sun Nov 25 09:12:36 2007 -0500
+++ b/include/geometry/kinetic-sort.h Tue Feb 19 06:56:49 2008 -0500
@@ -10,33 +10,35 @@
* 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 ElementIterator_ iterator over the underlying data structure that's kept in sorted order
* \arg TrajectoryExtractor_ applied to the iterator into SortDS_ should return a rational
* function describing the
* \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
*/
-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 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 +46,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 +63,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 +73,7 @@
private:
NodeList list_;
- SortDS* sort_;
- SwapCallback swap_callback_;
+ Swap swap_;
};
#include "kinetic-sort.hpp"
--- a/include/geometry/kinetic-sort.hpp Sun Nov 25 09:12:36 2007 -0500
+++ b/include/geometry/kinetic-sort.hpp Tue Feb 19 06:56:49 2008 -0500
@@ -1,49 +1,62 @@
-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)
+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 +74,15 @@
}
-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");
+
+ swap_(pos->element, simulator);
+
// Remove events
iterator prev = boost::prior(pos);
if (prev != list_.end())
@@ -88,9 +101,9 @@
//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;
@@ -126,9 +139,9 @@
}
-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;
@@ -150,9 +163,9 @@
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;
@@ -178,8 +191,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;
@@ -199,9 +212,9 @@
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;
@@ -209,9 +222,9 @@
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::
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::SwapEvent::
print(std::ostream& out) const
{
Parent::print(out) << ", SwapEvent at " << TrajectoryExtractor_()(position()->element);
--- a/include/geometry/polynomial.h Sun Nov 25 09:12:36 2007 -0500
+++ b/include/geometry/polynomial.h Tue Feb 19 06:56:49 2008 -0500
@@ -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>
@@ -30,20 +31,23 @@
static RootType root(const T& r) { return SynapsTraits<T>::root(r); }
static int sign_at(const RationalFunction& rf, const RootType& r);
static RootType between(const RootType& r1, const RootType& r2) { return SynapsTraits<T>::between(r1,r2); }
+ static bool sign_at_negative_infinity(const RationalFunction& 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 +55,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 +64,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 Sun Nov 25 09:12:36 2007 -0500
+++ b/include/geometry/polynomial.hpp Tue Feb 19 06:56:49 2008 -0500
@@ -5,8 +5,8 @@
{
typedef SYNAPS::Seq<RootType> RootSeq;
- RootSeq seq_num = SYNAPS::solve(rf.numerator(), Solver());
- RootSeq seq_den = SYNAPS::solve(rf.denominator(), Solver());
+ 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
@@ -32,3 +32,40 @@
{
return SynapsTraits<T>::sign_at(rf.numerator(), r) * SynapsTraits<T>::sign_at(rf.denominator(), r);
}
+
+template<class T>
+bool
+UPolynomial<T>::
+sign_at_negative_infinity(const RationalFunction& rf)
+{
+ const Polynomial& num = rf.numerator();
+ const Polynomial& den = rf.denominator();
+ unsigned int ndegree = num.get_degree();
+ unsigned int ddegree = den.get_degree();
+ return (((ndegree + 1) % 2 == 0) ^ (num[ndegree] > 0)) ^
+ (((ddegree + 1) % 2 == 0) ^ (den[ddegree] > 0));
+}
+
+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 Sun Nov 25 09:12:36 2007 -0500
+++ b/include/geometry/simulator.h Tue Feb 19 06:56:49 2008 -0500
@@ -5,8 +5,14 @@
#include "polynomial.h"
template<class Comparison> class IndirectComparison;
-template<class PolyKernel_, class Simulator_> class Event;
+/**
+ * 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 RationalFunction responsible for the event must be positive
+ * before the Event occurs.
+ */
template<class PolyKernel_, template<class Event> class EventComparison_ = std::less>
class Simulator
{
@@ -56,6 +62,12 @@
bool reached_infinity_;
};
+
+/**
+ * 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 PolyKernel_, template<class Event> class EventComparison_>
class Simulator<PolyKernel_, EventComparison_>::Event
{
@@ -80,15 +92,16 @@
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"; }
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 Sun Nov 25 09:12:36 2007 -0500
+++ b/include/geometry/simulator.hpp Tue Feb 19 06:56:49 2008 -0500
@@ -17,8 +17,13 @@
Event* ee = new Event_(e);
//std::cout << "Solving: " << f << std::endl;
PolynomialKernel::solve(f, ee->root_stack());
+ bool sign = PolynomialKernel::sign_at_negative_infinity(f);
while (!ee->root_stack().empty() && ee->root_stack().top() < current_time())
+ {
ee->root_stack().pop();
+ sign = !sign;
+ }
+ if (sign) ee->root_stack().pop(); // TODO: double-check the logic
//std::cout << "Pushing: " << ee->root_stack().top() << std::endl;
return queue_.push(ee);
}
--- a/include/topology/filtration.h Sun Nov 25 09:12:36 2007 -0500
+++ b/include/topology/filtration.h Tue Feb 19 06:56:49 2008 -0500
@@ -37,6 +37,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
--- a/tests/CMakeLists.txt Sun Nov 25 09:12:36 2007 -0500
+++ b/tests/CMakeLists.txt Tue Feb 19 06:56:49 2008 -0500
@@ -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 Tue Feb 19 06:56:49 2008 -0500
@@ -0,0 +1,12 @@
+set (targets
+ euclidean
+ polynomial
+ test-eventqueue
+ test-kinetic-sort
+ test-linalg)
+
+foreach (t ${targets})
+ add_executable (${t} ${t}.cpp)
+ target_link_libraries (${t} ${synaps_libraries})
+endforeach (t ${targets})
+
--- a/tests/geometry/Makefile Sun Nov 25 09:12:36 2007 -0500
+++ /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 Sun Nov 25 09:12:36 2007 -0500
+++ b/tests/geometry/euclidean.cpp Tue Feb 19 06:56:49 2008 -0500
@@ -1,4 +1,4 @@
-#include "euclidean.h"
+#include "geometry/euclidean.h"
#include <vector>
#include <iostream>
#include <cmath>
--- a/tests/geometry/polynomial.cpp Sun Nov 25 09:12:36 2007 -0500
+++ b/tests/geometry/polynomial.cpp Tue Feb 19 06:56:49 2008 -0500
@@ -1,5 +1,5 @@
-#include <euclidean.h>
-#include <polynomial.h>
+#include <geometry/euclidean.h>
+#include <geometry/polynomial.h>
#include <vector>
#include <iostream>
--- a/tests/geometry/test-eventqueue.cpp Sun Nov 25 09:12:36 2007 -0500
+++ b/tests/geometry/test-eventqueue.cpp Tue Feb 19 06:56:49 2008 -0500
@@ -1,4 +1,4 @@
-#include <eventqueue.h>
+#include <utilities/eventqueue.h>
#include <functional>
#include <iostream>
--- a/tests/geometry/test-kinetic-sort.cpp Sun Nov 25 09:12:36 2007 -0500
+++ b/tests/geometry/test-kinetic-sort.cpp Tue Feb 19 06:56:49 2008 -0500
@@ -1,25 +1,27 @@
-#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;
+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,7 +31,7 @@
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));
@@ -44,6 +46,9 @@
list.push_back(Polynomial("x^3 - 3"));
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 +57,14 @@
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)
+ while(!simulator.reached_infinity() && simulator.current_time() < 4)
{
- //std::cout << "Current time before: " << simulator.current_time() << std::endl;
- if (!ks.audit(&simulator)) return 1;
+ std::cout << "Current time before: " << simulator.current_time() << std::endl;
+ //if (!ks.audit(&simulator)) return 1;
//simulator.print(std::cout << "Auditing ");
simulator.process();
- //std::cout << "Current time after: " << simulator.current_time() << std::endl;
+ std::cout << "Current time after: " << simulator.current_time() << std::endl;
}
-
}
--- a/tests/geometry/test-linalg.cpp Sun Nov 25 09:12:36 2007 -0500
+++ b/tests/geometry/test-linalg.cpp Tue Feb 19 06:56:49 2008 -0500
@@ -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;