--- a/README Wed May 23 10:48:04 2007 -0400
+++ b/README Sun Jul 01 00:00:00 2007 -0400
@@ -5,6 +5,7 @@
boost - great set of C++ libraries
Doxygen - for building documentation
libcwd - for debugging only (is not needed by default)
+ synaps - for solving polynomial (for kinetic kernel), which in turn requires GMP
Configuration
The path to CGAL's Makefile is expected to be set in $CGAL_MAKEFILE, the rest
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/NOTES Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,69 @@
+== Predicates/Constructions
+ * circumsphere, side of circumsphere
+ * barycentric coordinates
+ * degenerate simplex
+ * projection on a k-flat, squared distance to the k-flat
+ * signed square distance to an oriented (d-1)-flat (?)
+
+
+== Computing circumcenter and circumradius of a k-simplex in d dimensions.
+
+ * Ken Clarkson's Gram-Schmidt-based algorithm works perfectly over any field.
+ The idea is as follows.
+
+ Given points p_i, the circumcenter c satisfies (p_i-c)^2 = r^2. By
+ translation, we can assume that some p_i is at the origin, so that c^c=r^2,
+ and that equation becomes p_i \cdot c - p_i^2/2 = 0. That is, [c -1/2] is
+ normal to [p_i p_i^2]. Also, c must be an affine combination of the p_i,
+ so if we start with a point [c' 0] with c' in the affine hull of the p_i,
+ and remove all multiples of all [p_i p_i^2] from it, via Gram-Schmidt for
+ example, and then scale the result so that the last entry is -1/2, we get
+ c (alternatively we can start with [0 1]). Note that if the simplex is
+ degenerate (its affine hull is less than k-dimensional), the last entry of
+ c is 0.
+
+ P = {p_0, ..., p_k}
+ p'_i = p_i - p_0 # make p_0 the origin
+ p''_i = [p'_i p'_i^2] # lift the points
+ c' = 0 # c is in the affine hull of p'_i
+ c'' = [c' 1] # c'' is not in the affine hull of p'_i^2
+
+ for i = 1 to k: # create an orthogonal basis for aff(p''_i)
+ a''_i = p''_i
+ for j = k-1 downto 1:
+ a''_i -= a''_j * (a''_i \cdot c''_j)/c''_j^2
+
+ for i = 1 to k:
+ c'' -= a''_i * (c'' \cdot a''_i) / (a''_i^2)
+ c' = c''[1..d]/(-2 * c''[d+1]) # scale c'' down
+
+ * Ulfar Erlingsson's, Erich Kaltofen's and David Musser's "Generic
+ Gram-Schmidt Orthogonalization by Exact Division" (ISSAC'96) gives an
+ algorithm to compute Gram-Schmidt for vector space with integral domain
+ coefficients. I need to understand the paper better to see the advantages
+ over normalizing rationals at every step of the computation above
+ (the paper's introduction says that it avoids "costly GCD computations").
+
+ * Using QR decomposition (or any other orthogonal decomposition), and
+ then using the formula of the circumsphere on the matrix R (which is
+ just the simplex rotated into a k-dimensional space):
+
+ * Only (d-2) determinants of a k by k matrix are required which is nice.
+
+ * All (implemented) computation of the QR decomposition that I have found
+ so far requires square roots (for norms), divisions, and often
+ comparisons. While divisions are acceptable, square roots and
+ comparisons aren't.
+
+ * There is a string of papers in late 80s, early 90s on square root free
+ and division free algorithms for QR decomposition out of VLSI community.
+ There is a paper that claims to put all such techniques into one
+ framework. There is also a paper that says that all such techniques don't
+ work, and square root implementation needs to be reconsidered. (authors
+ to look up if searching for these papers are K.J.R. Liu and E.
+ Frantzeskakis)
+
+ * Using the formulas from Herbert's Weighted Alpha Shapes paper which are
+ derived from formulas of the planes that are lifted Voronoi cells. This
+ approach requires (d-k)*(d+2) determinants of k by k matrices, which is not
+ great.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/euclidean.h Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,105 @@
+/*
+ * Author: Dmitriy Morozov
+ * Department of Computer Science, Duke University, 2005 -- 2007
+ */
+
+#ifndef __EUCLEDIAN_H__
+#define __EUCLEDIAN_H__
+
+#include <utility>
+#include <vector>
+#include <algorithm>
+
+#include "linalg.h"
+#include "number-traits.h"
+
+
+template<class NumberType_ = double>
+class Kernel
+{
+ public:
+ typedef unsigned int DimensionType;
+ typedef NumberType_ NumberType;
+ typedef LinearAlgebra<NumberType> LinearAlgebra;
+ typedef typename LinearAlgebra::MatrixType MatrixType;
+ typedef typename LinearAlgebra::VectorType VectorType;
+
+ class Point;
+ class Sphere;
+ typedef std::vector<const Point*> PointContainer;
+
+
+ Kernel(DimensionType dimension);
+
+ DimensionType dimension() const { return dimension_; }
+ Point origin() const { return origin_; }
+
+ /** Returns matrix describing the equation of a circumsphere of points */
+ Sphere circumsphere(const PointContainer& points) const;
+
+ /** Returns squared radius of the circumsphere */
+ NumberType circumradius(const PointContainer& points) const;
+
+ /** Returns center of the circumsphere */
+ Point circumcenter(const PointContainer& points) const;
+
+ /** The result is positive if points[0] lies outside the circumsphere of points,
+ 0 if points[0] is on the circumsphere, and negative if it's inside */
+ NumberType side_of_circumsphere(const PointContainer& points, const Point& p) const;
+
+ private:
+ NumberType& normalize(NumberType& n) const;
+ Point& normalize(Point& p) const;
+
+ DimensionType dimension_;
+ Point origin_;
+};
+
+
+/* Point */
+template<class NumberType_>
+class Kernel<NumberType_>::Point: public VectorType
+{
+ public:
+ typedef VectorType Parent;
+ typedef NumberType_ NumberType;
+
+ Point(DimensionType d): Parent(d) {}
+ template<class Vec> Point(const Vec& v): Parent(v) {}
+ Point(const Point& p, const NumberType& pp);
+
+
+ //operator VectorType() const { return *this; }
+
+ NumberType squared_distance(const Point& p) const;
+
+ using Parent::size;
+};
+
+
+/* Sphere */
+template<class NumberType_>
+class Kernel<NumberType_>::Sphere
+{
+ public:
+ Sphere(const Point& center,
+ const NumberType& squared_radius):
+ center_(center), squared_radius_(squared_radius)
+ {}
+
+ /** The result is positive if p lies outside the sphere,
+ 0 if p is on the sphere, and negative if it's inside */
+ NumberType side_of(const Point& p) const;
+
+ const Point& center() const { return center_; }
+ const NumberType& squared_radius() const { return squared_radius_; }
+
+ private:
+ Point center_;
+ NumberType squared_radius_;
+};
+
+
+#include "euclidean.hpp"
+
+#endif // __EUCLEDIAN_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/euclidean.hpp Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,141 @@
+/* --- Point --- */
+template<class NumberType_>
+Kernel<NumberType_>::
+Kernel(DimensionType dimension): dimension_(dimension), origin_(dimension)
+{
+ for (unsigned int i = 0; i < dimension; ++i)
+ origin_(i) = NumberType(0);
+}
+
+
+template<class NumberType_>
+Kernel<NumberType_>::Point::
+Point(const Point& p, const NumberType& pp): Parent(p.size() + 1)
+{
+ using boost::numeric::ublas::vector_range;
+ using boost::numeric::ublas::range;
+
+ vector_range<VectorType> vr(*this, range(0, size() - 1));
+ vr = p;
+ (*this)(size() - 1) = pp;
+}
+
+template<class NumberType_>
+typename Kernel<NumberType_>::NumberType
+Kernel<NumberType_>::Point::
+squared_distance(const Point& p) const
+{
+ return boost::numeric::ublas::inner_prod(*this - p, *this - p);
+}
+
+template<class NumberType_>
+std::istream&
+operator>>(std::istream& in, typename Kernel<NumberType_>::Point& p)
+{
+ for (unsigned int i = 0; i < p.size(); ++i)
+ in >> p[i];
+ return in;
+}
+
+template<class NumberType_>
+std::ostream&
+operator<<(std::ostream& out, typename Kernel<NumberType_>::Point& p)
+{
+ for (unsigned int i = 0; p.size(); ++i)
+ out << p[i];
+ return out;
+}
+
+/* --- Sphere --- */
+template<class NumberType_>
+typename Kernel<NumberType_>::NumberType
+Kernel<NumberType_>::Sphere::
+side_of(const Point& p) const
+{ return p.squared_distance(center_) - squared_radius(); }
+
+
+/* --- Kernel --- */
+template<class NumberType_>
+typename Kernel<NumberType_>::Sphere
+Kernel<NumberType_>::
+circumsphere(const PointContainer& points) const
+{
+ if (points.size() == 0) return Sphere(origin(), NumberType(0)); // FIXME: should this be an assertion instead?
+
+ using boost::numeric::ublas::inner_prod;
+ using boost::numeric::ublas::vector_range;
+ using boost::numeric::ublas::range;
+
+ std::vector<Point> basis;
+ basis.reserve(points.size() - 1);
+ for (unsigned int i = 1; i < points.size(); ++i)
+ {
+ Point pt = *points[i] - *points[0];
+ basis.push_back(Point(pt, inner_prod(pt, pt)));
+ for (unsigned int j = 0; j < i - 1; ++j)
+ {
+ basis[i-1] -= basis[j] * inner_prod(basis[i-1], basis[j]) / inner_prod(basis[j], basis[j]);
+ normalize(basis[i-1]);
+ }
+ }
+ Point clifted(origin(), NumberType(1));
+ for (unsigned int j = 0; j < basis.size(); ++j)
+ {
+ clifted -= basis[j] * inner_prod(clifted, basis[j]) / inner_prod(basis[j], basis[j]);
+ normalize(clifted);
+ }
+
+ Point center(vector_range<VectorType>(clifted, range(0, dimension())));
+ center /= NumberType_(-2) * clifted(dimension());
+ NumberType squared_radius = inner_prod(center,center);
+ center += *points[0];
+
+ normalize(center);
+ normalize(squared_radius);
+
+ return Sphere(center, squared_radius);
+}
+
+template<class NumberType_>
+typename Kernel<NumberType_>::NumberType
+Kernel<NumberType_>::
+circumradius(const PointContainer& points) const
+{
+ return circumsphere(points).squared_radius();
+}
+
+template<class NumberType_>
+typename Kernel<NumberType_>::Point
+Kernel<NumberType_>::
+circumcenter(const PointContainer& points) const
+{
+ return circumsphere(points).center();
+}
+
+template<class NumberType_>
+typename Kernel<NumberType_>::NumberType
+Kernel<NumberType_>::
+side_of_circumsphere(const PointContainer& points, const Point& p) const
+{
+ Sphere s = circumsphere(points);
+ return s.side_of(p);
+}
+
+/* --- Kernel Private --- */
+template<class NumberType_>
+typename Kernel<NumberType_>::NumberType&
+Kernel<NumberType_>::
+normalize(NumberType& n) const
+{
+ return number_traits<typename Point::NumberType>::normalize(n);
+}
+
+template<class NumberType_>
+typename Kernel<NumberType_>::Point&
+Kernel<NumberType_>::
+normalize(Point& p) const
+{
+ for (unsigned int i = 0; i < p.size(); ++i)
+ number_traits<typename Point::NumberType>::normalize(p(i));
+ return p;
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/eventqueue.h Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,71 @@
+#ifndef __EVENTQUEUE_H__
+#define __EVENTQUEUE_H__
+
+#include <list>
+#include <functional>
+#include <boost/utility.hpp>
+
+#include <iostream>
+#include <cassert> // TODO: switch to internal debugging system
+#include <string>
+
+// TODO: change inefficient list-based implementation to something heap-based
+// Need a queue that supports deleting arbitrary items (given by iterator),
+// and maintaining correctness of iterators when different operations (notably
+// remove and update are performed)
+
+template<class Event_, class EventComparison_>
+class EventQueue
+{
+
+ public:
+ typedef Event_ Event;
+ typedef EventComparison_ EventComparison;
+
+ typedef std::list<Event> QueueRepresentation;
+ typedef typename QueueRepresentation::iterator iterator;
+ typedef typename QueueRepresentation::const_iterator const_iterator;
+
+ EventQueue() {}
+
+ const_iterator top() const { assert(!empty()); return queue_.begin(); }
+ iterator top() { assert(!empty()); return queue_.begin(); }
+ iterator push(Event e) { queue_.push_front(e); iterator i = top(); update(i); return i; }
+ void pop() { assert(!empty()); queue_.erase(queue_.begin()); }
+ void remove(iterator i) { queue_.erase(i); }
+ void update(iterator i);
+
+ iterator end() { return queue_.end(); }
+ const_iterator end() const { return queue_.end(); }
+ bool empty() const { return queue_.empty(); }
+ size_t size() const { return queue_.size(); }
+
+ std::ostream& print(std::ostream& out, const std::string& prefix) const;
+
+ private:
+ QueueRepresentation queue_;
+};
+
+
+template<class Event_, class EventComparison_>
+void
+EventQueue<Event_, EventComparison_>::
+update(iterator i)
+{
+ QueueRepresentation tmp;
+ tmp.splice(tmp.end(), queue_, i);
+ iterator pos = std::find_if(queue_.begin(), queue_.end(), std::not1(std::bind2nd(EventComparison(), *i)));
+ queue_.splice(pos, tmp);
+}
+
+template<class Event_, class EventComparison_>
+std::ostream&
+EventQueue<Event_, EventComparison_>::
+print(std::ostream& out, const std::string& prefix) const
+{
+ for (typename QueueRepresentation::const_iterator cur = queue_.begin(); cur != queue_.end(); ++cur)
+ (*cur)->print(out << prefix) << std::endl;
+ return out;
+}
+
+#endif // __EVENTQUEUE_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/kinetic-sort.h Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,78 @@
+#ifndef __KINETIC_SORT_H__
+#define __KINETIC_SORT_H__
+
+#include <list>
+#include <boost/function.hpp>
+#include <boost/utility.hpp>
+#include <iostream>
+
+/**
+ * Maintains elements of the given data structure in the sorted order assuming the elements follow
+ * trajectories given by TrajectoryExtractor_.
+ *
+ * \arg SortDS_ should be forward and backward iterable, swaps are handles via SwapCallback
+ * \arg TrajectoryExtractor_ applied to the iterator into SortDS_ should return a rational
+ * function describing the
+ * \arg 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)
+ */
+template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+class KineticSort
+{
+ public:
+ typedef Simulator_ Simulator;
+ typedef typename Simulator::PolynomialKernel PolynomialKernel;
+ typedef SortDS_ SortDS;
+ typedef TrajectoryExtractor_ TrajectoryExtractor;
+
+ typedef typename Simulator::Key SimulatorKey;
+ typedef typename SortDS::iterator SortDSIterator;
+
+
+ private:
+ /* Implementation */
+ struct Node
+ {
+ SortDSIterator element;
+ SimulatorKey swap_event_key;
+
+ Node(SortDSIterator e, SimulatorKey k):
+ element(e), swap_event_key(k) {}
+ };
+
+ typedef std::list<Node> NodeList;
+
+ 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);
+
+ template<class InputIterator>
+ void insert(iterator pos, InputIterator f, InputIterator l, Simulator* simulator);
+ void erase(iterator pos, Simulator* simulator);
+ void update_trajectory(iterator pos, Simulator* simulator);
+
+ void swap(iterator pos, Simulator* simulator);
+
+ bool audit(Simulator* simulator) const;
+ /// @}
+
+ private:
+ class SwapEvent;
+ void schedule_swaps(iterator b, iterator e, Simulator* s);
+ void schedule_swaps(iterator i, Simulator* s);
+
+ private:
+ NodeList list_;
+ SortDS* sort_;
+ SwapCallback swap_callback_;
+};
+
+#include "kinetic-sort.hpp"
+
+#endif // __KINETIC_SORT_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/kinetic-sort.hpp Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,220 @@
+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)
+{
+ for (SortDSIterator cur = sort->begin(); cur != sort->end(); ++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>
+void
+KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+insert(iterator pos, InputIterator f, InputIterator l, Simulator* simulator)
+{
+ iterator previous = pos; --previous;
+ if (previous != list_.end()) simulator->remove(previous->swap_event_key);
+
+ sort_->insert(pos->element, f, l);
+
+ 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_>
+void
+KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+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_>
+void
+KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+update_trajectory(iterator pos, Simulator* simulator)
+{
+ iterator prev = boost::prior(pos);
+ if (prev != list_.end())
+ {
+ simulator->remove(prev->swap_event_key);
+ schedule_swaps(prev, simulator);
+ }
+
+ if (boost::next(pos) != list_.end())
+ {
+ simulator->remove(pos->swap_event_key);
+ schedule_swaps(pos, simulator);
+ }
+}
+
+
+template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+void
+KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+swap(iterator pos, Simulator* simulator)
+{
+ swap_callback_(sort_, pos->element);
+
+ // TODO: add assertion that boost::next(pos) != list_.end()
+
+ // Remove events
+ iterator prev = boost::prior(pos);
+ if (prev != list_.end())
+ simulator->remove(prev->swap_event_key);
+ iterator next = boost::next(pos);
+ simulator->remove(next->swap_event_key);
+
+ // Swap
+ list_.splice(pos, list_, next);
+
+ // update events
+ next->swap_event_key = pos->swap_event_key;
+ static_cast<SwapEvent*>(*(next->swap_event_key))->set_position(next);
+ schedule_swaps(prev, simulator);
+ schedule_swaps(pos, simulator);
+ //audit(simulator);
+}
+
+template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+bool
+KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+audit(Simulator* simulator) const
+{
+ typedef typename Simulator::RationalFunction RationalFunction;
+ typedef typename Simulator::Time Time;
+
+ Time t = simulator->audit_time();
+ std::cout << "Auditing at " << t << std::endl;
+
+ TrajectoryExtractor te;
+
+ typename NodeList::const_iterator next = list_.begin();
+ typename NodeList::const_iterator cur = next++;
+ RationalFunction cur_trajectory = te(cur->element);
+ while (next != list_.end())
+ {
+ (*(cur->swap_event_key))->print(std::cout << " ") << std::endl;
+
+ RationalFunction next_trajectory = te(next->element);
+ std::cout << " Auditing: " << cur_trajectory << ", " << next_trajectory << std::endl;
+ std::cout << " Difference: " << next_trajectory - cur_trajectory << std::endl;
+ std::cout << " Sign at: " << t << ", " << PolynomialKernel::sign_at(next_trajectory - cur_trajectory, t) << std::endl;
+ if (PolynomialKernel::sign_at(next_trajectory - cur_trajectory, t) == -1)
+ {
+ std::cout << "Audit failed at " << *cur->element << ", " << *next->element << std::endl;
+ return false;
+ }
+
+ cur_trajectory = next_trajectory;
+ cur = next++;
+ }
+ if (cur != list_.end()) (*(cur->swap_event_key))->print(std::cout << " ") << std::endl;
+ return true;
+}
+
+
+template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+void
+KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+schedule_swaps(iterator b, iterator e, Simulator* simulator)
+{
+ typedef typename Simulator::RationalFunction RationalFunction;
+
+ TrajectoryExtractor te;
+
+ iterator next = b;
+ iterator cur = next++;
+ RationalFunction cur_trajectory = te(cur->element);
+ while (next != e)
+ {
+ RationalFunction next_trajectory = te(next->element);
+ std::cout << "Next trajectory: " << next_trajectory << std::endl;
+ // TODO: add assertion that (next_trajectory - cur_trajectory)(s->curren_time()) > 0
+ cur->swap_event_key = simulator->add(next_trajectory - cur_trajectory, SwapEvent(this, cur));
+ cur = next++;
+ cur_trajectory = next_trajectory;
+ }
+ if (cur != e) schedule_swaps(cur, simulator);
+}
+
+template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+void
+KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+schedule_swaps(iterator i, Simulator* simulator)
+{
+ typedef typename Simulator::RationalFunction RationalFunction;
+
+ if (i == list_.end()) return;
+ if (boost::next(i) == list_.end())
+ {
+ i->swap_event_key = simulator->add(SwapEvent(this, i));
+ return;
+ }
+
+ TrajectoryExtractor te;
+
+ iterator next = boost::next(i);
+ RationalFunction i_trajectory = te(i->element);
+ RationalFunction 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;
+
+ i->swap_event_key = simulator->add(next_trajectory - i_trajectory, SwapEvent(this, i));
+ //i->swap_event_key = simulator->add(next_trajectory, SwapEvent(this, i));
+}
+
+/* SwapEvent */
+template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+class KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::SwapEvent: public Simulator::Event
+{
+ public:
+ typedef typename Simulator::Event Parent;
+
+ SwapEvent(const SwapEvent& e):
+ sort_(e.sort_), pos_(e.pos_) {}
+ SwapEvent(KineticSort* sort, iterator pos):
+ sort_(sort), pos_(pos) {}
+
+ virtual bool process(Simulator* s) const;
+ void set_position(iterator i) { pos_ = i; }
+ iterator position() const { return pos_; }
+ std::ostream& print(std::ostream& out) const;
+
+ private:
+ KineticSort* sort_;
+ iterator pos_;
+};
+
+template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+bool
+KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::SwapEvent::
+process(Simulator* s) const
+{
+ std::cout << "Swapping. Current time: " << s->current_time() << std::endl;
+ sort_->swap(pos_, s);
+ return true;
+}
+
+template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+std::ostream&
+KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::SwapEvent::
+print(std::ostream& out) const
+{
+ Parent::print(out) << ", SwapEvent at " << TrajectoryExtractor_()(position()->element);
+ return out;
+}
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/linalg.h Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,31 @@
+#ifndef __LINALG_H__
+#define __LINALG_H__
+
+#include <vector>
+
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/lu.hpp>
+#include <boost/numeric/ublas/io.hpp>
+
+template<class ValueType_>
+class LinearAlgebra
+{
+ public:
+ typedef ValueType_ ValueType;
+ typedef boost::numeric::ublas::matrix<ValueType> MatrixType;
+ typedef boost::numeric::ublas::vector<ValueType> VectorType;
+
+
+ /* Currently don't need any of this */
+ static ValueType determinant(const MatrixType& a);
+ static void solve(const MatrixType& a, const VectorType& b, VectorType& x);
+
+ private:
+ template<class TriangularType_>
+ static ValueType determinant(const boost::numeric::ublas::triangular_adaptor<MatrixType, TriangularType_>& t);
+};
+
+#include "linalg.hpp"
+
+#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/linalg.hpp Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,38 @@
+template<class ValueType_>
+typename LinearAlgebra<ValueType_>::ValueType
+LinearAlgebra<ValueType_>::
+determinant(const MatrixType& a)
+{
+ using namespace boost::numeric::ublas;
+ MatrixType m = a;
+ lu_factorize(m);
+ return determinant(triangular_adaptor<MatrixType, upper>(m));
+}
+
+template<class ValueType_>
+void
+LinearAlgebra<ValueType_>::
+solve(const MatrixType& a, const VectorType& b, VectorType& x)
+{
+ using namespace boost::numeric::ublas;
+ MatrixType m = a;
+ x = b;
+ lu_factorize(m);
+ //lu_substitute(m, x);
+ inplace_solve(m,x,unit_lower_tag());
+ inplace_solve(m,x,upper_tag());
+}
+
+
+/* Private */
+template<class ValueType_>
+template<class TriangularType_>
+typename LinearAlgebra<ValueType_>::ValueType
+LinearAlgebra<ValueType_>::
+determinant(const boost::numeric::ublas::triangular_adaptor<MatrixType, TriangularType_>& t)
+{
+ ValueType res = ValueType(1);
+ for (typename MatrixType::size_type i = 0; i < t.size1(); ++i)
+ res *= t(i,i);
+ return res;
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/number-traits.h Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,13 @@
+#ifndef __NUMBER_TRAITS_H__
+#define __NUMBER_TRAITS_H__
+
+template<class NumberType_>
+class number_traits
+{
+ public:
+ typedef NumberType_ NumberType;
+
+ static NumberType& normalize(NumberType& n) { return n; }
+};
+
+#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/polynomial.h Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,66 @@
+#ifndef __POLYNOMIAL_H__
+#define __POLYNOMIAL_H__
+
+#include <synaps/upol.h>
+#include <synaps/upol/gcd.h>
+#include <synaps/usolve/Algebraic.h>
+#include <synaps/usolve/algebraic/sign_fct.h>
+#include <synaps/usolve/bezier/SlvBzStd.h>
+//#include <synaps/usolve/Sturm.h>
+//#include <synaps/arithm/Infinity.h>
+
+#include <stack>
+
+#include "rational-function.h"
+
+template<class T> class SynapsTraits;
+
+template<class T>
+class UPolynomial
+{
+ public:
+ typedef typename SynapsTraits<T>::Polynomial Polynomial;
+ typedef RationalFunction<Polynomial> RationalFunction;
+
+ typedef typename SynapsTraits<T>::Solver Solver;
+ typedef typename SynapsTraits<T>::RootType RootType;
+ typedef std::stack<RootType> RootStack;
+
+ static void solve(const RationalFunction& rf, RootStack& stack);
+ static 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); }
+};
+
+template<class T>
+struct SynapsTraits ///< Suitable for double
+{
+ typedef T CoefficientType;
+ typedef SYNAPS::UPolDse<CoefficientType> Polynomial;
+ typedef SYNAPS::SlvBzStd<CoefficientType> Solver;
+ typedef T RootType;
+
+ 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; }
+};
+
+template<>
+struct SynapsTraits<ZZ>
+{
+ typedef ZZ CoefficientType;
+ typedef SYNAPS::UPolDse<CoefficientType> Polynomial;
+ typedef SYNAPS::Algebraic<CoefficientType> 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);}
+ 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(p, r); }
+ //static RootType between(const RootType& r1, const RootType& r2) { RootType r = r1; r += r2; r /= root(2); return r; }
+};
+
+#include "polynomial.hpp"
+
+#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/polynomial.hpp Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,34 @@
+template<class T>
+void
+UPolynomial<T>::
+solve(const RationalFunction& rf, RootStack& stack)
+{
+ typedef SYNAPS::Seq<RootType> RootSeq;
+
+ RootSeq seq_num = SYNAPS::solve(rf.numerator(), Solver());
+ RootSeq seq_den = SYNAPS::solve(rf.denominator(), Solver());
+
+ // TODO: assert that all roots in seq_den have positive multiplicity
+ // TODO: deal with multiplicities for the numerator
+ for (typename RootSeq::const_reverse_iterator cur = seq_num.rbegin();
+ cur != seq_num.rend();
+ ++cur)
+ {
+ if (SynapsTraits<T>::multiplicity(*cur) % 2 != 0)
+ {
+ if (!stack.empty() && stack.top() == *cur) // get rid of even multiplicities
+ // TODO: add logging information for this
+ stack.pop();
+ else
+ stack.push(*cur);
+ }
+ }
+}
+
+template<class T>
+int
+UPolynomial<T>::
+sign_at(const RationalFunction& rf, const RootType& r)
+{
+ return SynapsTraits<T>::sign_at(rf.numerator(), r) * SynapsTraits<T>::sign_at(rf.denominator(), r);
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/rational-function.h Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,107 @@
+#ifndef __RATIONAL_FUNCTION_H__
+#define __RATIONAL_FUNCTION_H__
+
+#include <iostream>
+#include "number-traits.h"
+
+template <class Polynomial_>
+class RationalFunction
+{
+ public:
+ typedef Polynomial_ Polynomial;
+ typedef typename Polynomial::coeff_t CoefficientType;
+ typedef typename Polynomial::value_type ValueType;
+
+ /// \name Constructors
+ /// @{
+ RationalFunction():
+ numerator_(CoefficientType(0)),
+ denominator_(CoefficientType(1)) {}
+ RationalFunction(const Polynomial& p):
+ numerator_(p),
+ denominator_(CoefficientType(1)) {}
+ RationalFunction(const Polynomial& num, const Polynomial& denom):
+ numerator_(num), denominator_(denom) { normalize(); }
+ RationalFunction(const RationalFunction& other):
+ numerator_(other.numerator_),
+ denominator_(other.denominator_) { normalize(); }
+ /// @}
+
+ /// \name Operators
+ /// @{
+ RationalFunction operator-() const;
+ RationalFunction operator+(const RationalFunction& o) const;
+ RationalFunction operator-(const RationalFunction& o) const;
+ RationalFunction operator*(const RationalFunction& o) const;
+ RationalFunction operator/(const RationalFunction& o) const;
+ RationalFunction operator+(const CoefficientType& a) const;
+ RationalFunction operator-(const CoefficientType& a) const;
+ RationalFunction operator*(const CoefficientType& a) const;
+ RationalFunction operator/(const CoefficientType& a) const;
+ /// @}
+
+ /// \name Modifiers
+ /// @{
+ RationalFunction& operator+=(const RationalFunction& o);
+ RationalFunction& operator-=(const RationalFunction& o);
+ RationalFunction& operator*=(const RationalFunction& o);
+ RationalFunction& operator/=(const RationalFunction& o);
+ /// @}
+
+ /// \name Assignment
+ /// @{
+ RationalFunction& operator=(const RationalFunction& o);
+ //RationalFunction& operator=(const Polynomial& o);
+ /// @}
+
+ /// \name Evaluation
+ /// @{
+ ValueType operator()(const ValueType& t) const;
+ bool operator==(const RationalFunction& o) const;
+ bool operator!=(const RationalFunction& o) const { return !operator==(o); }
+ /// @}
+
+ /// \name Accessors
+ /// @{
+ const Polynomial& numerator() const { return numerator_; }
+ const Polynomial& denominator() const { return denominator_; }
+ /// @}
+
+ RationalFunction& normalize();
+
+ private:
+ Polynomial numerator_, denominator_;
+};
+
+template<class P>
+struct number_traits<RationalFunction<P> >
+{
+ typedef RationalFunction<P> NumberType;
+ static NumberType& normalize(NumberType& n) { return n.normalize(); }
+};
+
+
+template<class Polynomial_>
+std::ostream&
+operator<<(std::ostream& out, const RationalFunction<Polynomial_>& r)
+{ return out << r.numerator() << " / " << r.denominator(); }// << ", gcd: " << gcd(r.numerator(), r.denominator()); }
+
+template<class Polynomial_>
+inline RationalFunction<Polynomial_>
+operator*(const typename RationalFunction<Polynomial_>::CoefficientType& a, const RationalFunction<Polynomial_>& r)
+{ return (r * a); }
+
+template<class Polynomial_>
+inline RationalFunction<Polynomial_>
+operator+(const typename RationalFunction<Polynomial_>::CoefficientType& a, const RationalFunction<Polynomial_>& r)
+{ return (r + a); }
+
+template<class Polynomial_>
+inline RationalFunction<Polynomial_>
+operator-(const typename RationalFunction<Polynomial_>::NT& a, const RationalFunction<Polynomial_>& r)
+{ return -(r - a); }
+
+
+#include "rational-function.hpp"
+
+#endif // __RATIONAL_FUNCTION_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/rational-function.hpp Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,150 @@
+/* Operators */
+template<class P>
+RationalFunction<P>
+RationalFunction<P>::
+operator-() const
+{ return RationalFunction(-numerator_,denominator_); }
+
+template<class P>
+RationalFunction<P>
+RationalFunction<P>::
+operator+(const RationalFunction& o) const
+{ return RationalFunction(numerator_*o.denominator_ + o.numerator_*denominator_, denominator_*o.denominator_); }
+
+template<class P>
+RationalFunction<P>
+RationalFunction<P>::
+operator-(const RationalFunction& o) const
+{ RationalFunction tmp(*this); tmp -= o; return tmp; }
+//{ return RationalFunction(numerator_*o.denominator_ - o.numerator_*denominator_, denominator_*o.denominator_); }
+
+template<class P>
+RationalFunction<P>
+RationalFunction<P>::
+operator*(const RationalFunction& o) const
+{ return RationalFunction(numerator_*o.numerator_, denominator_*o.denominator_); }
+
+template<class P>
+RationalFunction<P>
+RationalFunction<P>::
+operator/(const RationalFunction& o) const
+{ return RationalFunction(numerator_*o.denominator_, denominator_*o.numerator_); }
+
+template<class P>
+RationalFunction<P>
+RationalFunction<P>::
+operator+(const typename RationalFunction<P>::CoefficientType& a) const
+{ return RationalFunction(numerator_ + a*denominator_, denominator_); }
+
+template<class P>
+RationalFunction<P>
+RationalFunction<P>::
+operator-(const typename RationalFunction<P>::CoefficientType& a) const
+{ return operator+(-a); }
+
+template<class P>
+RationalFunction<P>
+RationalFunction<P>::
+operator*(const typename RationalFunction<P>::CoefficientType& a) const
+{ return RationalFunction(a*numerator_, denominator_); }
+
+template<class P>
+RationalFunction<P>
+RationalFunction<P>::
+operator/(const typename RationalFunction<P>::CoefficientType& a) const
+{ return RationalFunction(numerator_, a*denominator_); }
+
+template<class P>
+RationalFunction<P>&
+RationalFunction<P>::
+operator+=(const RationalFunction& o)
+{
+ numerator_ *= o.denominator_;
+ numerator_ += o.numerator_*denominator_;
+ denominator_ *= o.denominator_;
+ return *this;
+}
+
+template<class P>
+RationalFunction<P>&
+RationalFunction<P>::
+operator-=(const RationalFunction& o)
+{
+ numerator_ *= o.denominator_;
+ numerator_ -= o.numerator_*denominator_;
+ denominator_ *= o.denominator_;
+ return *this;
+}
+
+template<class P>
+RationalFunction<P>&
+RationalFunction<P>::
+operator*=(const RationalFunction& o)
+{
+ numerator_ *= o.numerator_;
+ denominator_ *= o.denominator_;
+ return *this;
+}
+
+template<class P>
+RationalFunction<P>&
+RationalFunction<P>::
+operator/=(const RationalFunction& o)
+{
+ numerator_ *= o.denominator_;
+ denominator_ *= o.numerator_;
+ return *this;
+}
+
+template<class P>
+RationalFunction<P>&
+RationalFunction<P>::
+operator=(const RationalFunction& o)
+{
+ numerator_ = o.numerator_;
+ denominator_ = o.denominator_;
+ return *this;
+}
+
+#if 0
+template<class P>
+RationalFunction<P>&
+RationalFunction<P>::
+operator=(const Polynomial& o)
+{
+ numerator_ = o;
+ denominator_ = 1;
+ return *this;
+}
+#endif
+
+/* Evaluation */
+template<class P>
+typename RationalFunction<P>::ValueType
+RationalFunction<P>::
+operator()(const typename RationalFunction<P>::ValueType& t) const
+{ return numerator_(t)/denominator_(t); }
+
+template<class P>
+bool
+RationalFunction<P>::
+operator==(const RationalFunction& o) const
+{ return (numerator_ == o.numerator_) && (denominator_ == o.denominator_); }
+
+template<class P>
+RationalFunction<P>&
+RationalFunction<P>::
+normalize()
+{
+#if 0
+ std::cout << "This: " << std::flush << this << std::endl;
+ std::cout << "Numerator address: " << std::flush << &numerator_ << std::endl;
+ std::cout << "Denominator address: " << std::flush << &denominator_ << std::endl;
+ std::cout << "Normalizing (numerator): " << std::flush << numerator_ << std::endl;
+ std::cout << "Normalizing (denominator): " << std::flush << denominator_ << std::endl;
+#endif
+ Polynomial divisor = gcd(numerator_, denominator_);
+ numerator_ /= divisor;
+ denominator_ /= divisor;
+ return *this;
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/simulator.h Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,108 @@
+#ifndef __SIMULATOR_H__
+#define __SIMULATOR_H__
+
+#include "eventqueue.h"
+#include "polynomial.h"
+
+template<class Comparison> class IndirectComparison;
+template<class PolyKernel_, class Simulator_> class Event;
+
+template<class PolyKernel_, template<class Event> class EventComparison_ = std::less>
+class Simulator
+{
+ public:
+ typedef PolyKernel_ PolynomialKernel;
+ typedef typename PolynomialKernel::Polynomial Polynomial;
+ typedef typename PolynomialKernel::RationalFunction RationalFunction;
+ typedef typename PolynomialKernel::RootStack RootStack;
+ typedef typename PolynomialKernel::RootType RootType;
+ typedef RootType Time;
+
+ class Event;
+ typedef EventComparison_<Event> EventComparison;
+ typedef EventQueue<Event*, IndirectComparison<EventComparison> >
+ EventQueue;
+ typedef typename EventQueue::iterator Key;
+ typedef typename EventQueue::const_iterator const_Key;
+
+
+ Simulator(Time start = PolynomialKernel::root(0)):
+ current_(start),
+ reached_infinity_(false) {}
+
+
+ template<class Event_>
+ Key add(const Event_& e);
+ template<class Event_>
+ Key add(const RationalFunction& f, const Event_& e);
+ void process();
+ void update(Key k, const RationalFunction& f);
+
+ void remove(Key k) { queue_.remove(k); }
+ Key null_key() { return queue_.end(); }
+
+ Time current_time() const { return current_; }
+ Time audit_time() const;
+ bool reached_infinity() const { return reached_infinity_; }
+
+ std::ostream& print(std::ostream& out) const;
+
+ private:
+ void update(Key i);
+
+ private:
+ Time current_;
+ EventQueue queue_;
+ bool reached_infinity_;
+};
+
+template<class PolyKernel_, template<class Event> class EventComparison_>
+class Simulator<PolyKernel_, EventComparison_>::Event
+{
+ public:
+ typedef PolyKernel_ PolynomialKernel;
+ typedef typename PolynomialKernel::RootStack RootStack;
+
+ /// 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();
+ }
+
+ 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_;
+};
+
+template<class Comparison_>
+class IndirectComparison: public std::binary_function<const typename Comparison_::first_argument_type*,
+ const typename Comparison_::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;
+
+ bool operator()(first_argument_type e1, second_argument_type e2) const
+ { return Comparison()(*e1, *e2); }
+};
+
+#include "simulator.hpp"
+
+#endif // __SIMULATOR_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/simulator.hpp Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,82 @@
+template<class PolyKernel_, template<class Event> class EventComparison_>
+template<class Event_>
+typename Simulator<PolyKernel_, EventComparison_>::Key
+Simulator<PolyKernel_, EventComparison_>::
+add(const Event_& e)
+{
+ Event* ee = new Event_(e);
+ return queue_.push(ee);
+}
+
+template<class PolyKernel_, template<class Event> class EventComparison_>
+template<class Event_>
+typename Simulator<PolyKernel_, EventComparison_>::Key
+Simulator<PolyKernel_, EventComparison_>::
+add(const RationalFunction& f, const Event_& e)
+{
+ Event* ee = new Event_(e);
+ //std::cout << "Solving: " << f << std::endl;
+ PolynomialKernel::solve(f, ee->root_stack());
+ while (!ee->root_stack().empty() && ee->root_stack().top() < current_time())
+ ee->root_stack().pop();
+ //std::cout << "Pushing: " << ee->root_stack().top() << std::endl;
+ return queue_.push(ee);
+}
+
+template<class PolyKernel_, template<class Event> class EventComparison_>
+void
+Simulator<PolyKernel_, EventComparison_>::
+update(Key k, const RationalFunction& f)
+{
+ Event* ee = *k;
+ ee->root_stack() = RootStack(); // no clear() in std::stack
+ PolynomialKernel::solve(f, ee->root_stack());
+ while (!ee->root_stack().empty() && ee->root_stack().top() < current_time())
+ ee->root_stack().pop();
+ update(k);
+}
+
+template<class PolyKernel_, template<class Event> class EventComparison_>
+void
+Simulator<PolyKernel_, EventComparison_>::
+process()
+{
+ std::cout << "Queue size: " << queue_.size() << std::endl;
+ Key top = queue_.top();
+ Event* e = *top;
+
+ if (e->root_stack().empty()) { reached_infinity_ = true; return; }
+ else { current_ = e->root_stack().top(); e->root_stack().pop(); }
+
+ if (e->process(this)) update(top);
+ else { queue_.pop(); delete e; }
+}
+
+template<class PolyKernel_, template<class Event> class EventComparison_>
+void
+Simulator<PolyKernel_, EventComparison_>::
+update(Key i)
+{
+ queue_.update(i);
+}
+
+template<class PolyKernel_, template<class Event> class EventComparison_>
+typename Simulator<PolyKernel_, EventComparison_>::Time
+Simulator<PolyKernel_, EventComparison_>::
+audit_time() const
+{
+ const_Key top = queue_.top();
+ Event* e = *top;
+
+ if (e->root_stack().empty()) return current_ + 1;
+ else return PolynomialKernel::between(e->root_stack().top(), current_);
+}
+
+template<class PolyKernel_, template<class Event> class EventComparison_>
+std::ostream&
+Simulator<PolyKernel_, EventComparison_>::
+print(std::ostream& out) const
+{
+ out << "Simulator: " << std::endl;
+ return queue_.print(out, " ");
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/geometry/Makefile Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,26 @@
+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}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/geometry/NOTES Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,5 @@
+Write an audit for the data structure.
+
+-DBOOST_UBLAS_TYPE_CHECK=0 is necessary since otherwise Boost's uBLAS gives a compile-time error
+when it doesn't find sqrt for polynomials. It doesn't actually need it (e.g., for LU-decomposition,
+but it looks for it during type check). (File a bug report with Boost?)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/geometry/euclidean.cpp Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,86 @@
+#include "euclidean.h"
+#include <vector>
+#include <iostream>
+#include <cmath>
+
+typedef Kernel<double> K;
+typedef K::Point Point;
+typedef K::Sphere Sphere;
+typedef K::PointContainer PointContainer;
+typedef K::MatrixType MatrixType;
+
+int main()
+{
+ K k(3);
+ std::vector<Point> points(6, k.origin());
+ points[0][0] = 0; points[0][1] = 0; points[0][2] = 0;
+ points[1][0] = 0; points[1][1] = 2; points[1][2] = 0;
+ points[2][0] = 0; points[2][1] = 0; points[2][2] = 5;
+ points[3][0] = 1; points[3][1] = 1; points[3][2] = 1;
+ points[4][0] = 0; points[4][1] = 1; points[4][2] = 0;
+ points[5][0] = 1; points[5][1] = 0; points[5][2] = 0;
+
+
+ // Edges
+ {
+ PointContainer vertices(2);
+ vertices[0] = &points[0]; vertices[1] = &points[2];
+ std::cout << "{0, 2}:" << std::endl;
+ Sphere s = k.circumsphere(vertices);
+ std::cout << "Circumsphere: " << s.center() << " " << s.squared_radius() << std::endl;
+ std::cout << "Side of: " << k.side_of_circumsphere(vertices, *vertices[1]) << std::endl;
+ std::cout << std::endl;
+
+ vertices[0] = &points[0]; vertices[1] = &points[3];
+ std::cout << "{0, 3}:" << std::endl;
+ s = k.circumsphere(vertices);
+ std::cout << "Circumsphere: " << s.center() << " " << s.squared_radius() << std::endl;
+ std::cout << "Side of: " << k.side_of_circumsphere(vertices, *vertices[1]) << std::endl;
+ std::cout << std::endl;
+ }
+
+ // Triangles
+ {
+ PointContainer vertices(3);
+ vertices[0] = &points[0]; vertices[1] = &points[3]; vertices[2] = &points[1];
+ std::cout << "{0, 3, 1}:" << std::endl;;
+ Sphere s = k.circumsphere(vertices);
+ std::cout << "Circumsphere: " << s.center() << " " << s.squared_radius() << std::endl;
+ std::cout << "Side of: " << k.side_of_circumsphere(vertices, *vertices[1]) << std::endl;
+ std::cout << std::endl;
+
+ vertices[0] = &points[0]; vertices[1] = &points[4]; vertices[2] = &points[5];
+ std::cout << "{0, 4, 5}:" << std::endl;
+ s = k.circumsphere(vertices);
+ std::cout << "Circumsphere: " << s.center() << " " << s.squared_radius() << std::endl;
+ std::cout << "Side of: " << k.side_of_circumsphere(vertices, *vertices[1]) << std::endl;
+ std::cout << std::endl;
+ }
+
+ // Tetrahedron
+ {
+ PointContainer vertices(4);
+ vertices[0] = &points[3]; vertices[1] = &points[1]; vertices[2] = &points[2]; vertices[3] = &points[0];
+ std::cout << "{3, 1, 2, 0}:" << std::endl;
+ Sphere s = k.circumsphere(vertices);
+ std::cout << "Circumsphere: " << s.center() << " " << s.squared_radius() << std::endl;
+ std::cout << "Side of: " << k.side_of_circumsphere(vertices, *vertices[1]) << std::endl;
+
+ std::cout << s.center().squared_distance(points[0]) << std::endl;
+ std::cout << std::endl;
+ }
+
+ {
+ PointContainer vertices(3);
+ vertices[0] = &points[3]; vertices[1] = &points[1]; vertices[2] = &points[2];
+ std::cout << "{3, 1, 2}:" << std::endl;
+ Sphere s = k.circumsphere(vertices);
+ std::cout << "Circumsphere: " << s.center() << " " << s.squared_radius() << std::endl;
+ std::cout << "Side of: " << k.side_of_circumsphere(vertices, points[0]) << std::endl;
+
+ std::cout << "Distance: " << points[0].squared_distance(s.center()) << std::endl;
+
+ std::cout << s.center().squared_distance(points[0]) << std::endl;
+ std::cout << std::endl;
+ }
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/geometry/polynomial.cpp Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,140 @@
+#include <euclidean.h>
+#include <polynomial.h>
+
+#include <vector>
+#include <iostream>
+#include <cmath>
+
+typedef UPolynomial<ZZ> PolynomialKernel;
+typedef PolynomialKernel::Polynomial Polynomial;
+typedef PolynomialKernel::RationalFunction RationalF;
+
+typedef Kernel<RationalF> K;
+typedef K::Point Point;
+typedef K::Sphere Sphere;
+typedef K::PointContainer PointContainer;
+typedef K::MatrixType MatrixType;
+
+int main()
+{
+ K k(3);
+ std::vector<Point> points(7, k.origin());
+ points[0][0] = Polynomial(0); points[0][1] = Polynomial(0); points[0][2] = Polynomial(0);
+ points[1][0] = Polynomial(0); points[1][1] = Polynomial("x+2"); points[1][2] = Polynomial(0);
+ points[2][0] = Polynomial(0); points[2][1] = Polynomial(0); points[2][2] = Polynomial("x^2 + 5");
+ points[3][0] = Polynomial("x^3"); points[3][1] = Polynomial(1); points[3][2] = Polynomial("x");
+ points[4][0] = Polynomial(0); points[4][1] = Polynomial("x^2 + 2*x + 5"); points[4][2] = Polynomial(0);
+ points[5][0] = Polynomial("x^3 + 3*x + 7"); points[5][1] = Polynomial(0); points[5][2] = Polynomial(0);
+ points[6][0] = Polynomial(0); points[6][1] = Polynomial("x + 6"); points[6][2] = Polynomial("x");
+
+ // Solving polynomials
+ {
+ PolynomialKernel::RootStack roots;
+ std::cout << "Solving " << points[5][0] << ": " << std::endl;
+ PolynomialKernel::solve(points[5][0], roots);
+ while (!roots.empty()) { std::cout << roots.top() << std::endl; roots.pop(); }
+ }
+
+ {
+ Polynomial p("x^3 - 2*x + 1");
+ PolynomialKernel::RootStack roots;
+ std::cout << "Solving " << p << ": " << std::endl;
+ PolynomialKernel::solve(p, roots);
+ while (!roots.empty()) { std::cout << roots.top() << std::endl; roots.pop(); }
+ }
+
+#if 0
+ // FIXME: explore
+ {
+ UPolynomial<QQ>::Polynomial p("1.2*x + 3.67");
+ UPolynomial<QQ>::RootStack roots;
+ UPolynomial<QQ>::solve(p, roots);
+ while (!roots.empty()) { std::cout << roots.top() << std::endl; roots.pop(); }
+ }
+#endif
+
+ {
+ RationalF r1 = Polynomial("2*x - 4");
+ RationalF r2 = Polynomial("x^3 - 3");
+ RationalF r3 = Polynomial("x^2 - 3*x^3");
+ std::cout << r2 - r1 << std::endl;
+ std::cout << RationalF(Polynomial("2*x"), Polynomial(1)*Polynomial(1)) << std::endl;
+
+ PolynomialKernel::RootStack roots;
+ std::cout << "Solving " << (r2 - r1) << ": " << std::endl;
+ PolynomialKernel::solve(r2 - r1, roots);
+ while (!roots.empty()) { std::cout << roots.top() << std::endl; roots.pop(); }
+
+ std::cout << "Solving " << (r3 - r1) << ": " << std::endl;
+ PolynomialKernel::solve(r3 - r1, roots);
+ while (!roots.empty()) { std::cout << roots.top() << std::endl; roots.pop(); }
+
+ std::cout << "Solving " << (r3 - r2) << ": " << std::endl;
+ PolynomialKernel::solve(r3 - r2, roots);
+ //std::cout << "Sign of r3 at " << roots.top() << " is " << PolynomialKernel::sign_at(r3, roots.top()) << std::endl;
+ while (!roots.empty()) { std::cout << roots.top() << std::endl; roots.pop(); }
+ }
+
+ return 0;
+
+ // Edges
+ {
+ PointContainer vertices(2);
+ vertices[0] = &points[0]; vertices[1] = &points[2];
+ std::cout << "{0, 2}:" << std::endl;
+ Sphere s = k.circumsphere(vertices);
+ std::cout << "Circumsphere: " << s.center() << " " << s.squared_radius() << std::endl;
+ std::cout << "Side of: " << k.side_of_circumsphere(vertices, *vertices[1]) << std::endl;
+
+ vertices[0] = &points[0]; vertices[1] = &points[3];
+ std::cout << "{0, 3}:" << std::endl;
+ s = k.circumsphere(vertices);
+ std::cout << "Circumsphere: " << s.center() << " " << s.squared_radius() << std::endl;
+ std::cout << "Side of: " << k.side_of_circumsphere(vertices, *vertices[1]) << std::endl;
+ }
+
+#if 1
+ // Triangles
+ {
+ PointContainer vertices(3);
+ vertices[0] = &points[0]; vertices[1] = &points[3]; vertices[2] = &points[1];
+ std::cout << "{0, 3, 1}:" << std::endl;;
+ Sphere s = k.circumsphere(vertices);
+ std::cout << "Circumsphere: " << s.center() << " " << s.squared_radius() << std::endl;
+ std::cout << "Side of: " << k.side_of_circumsphere(vertices, *vertices[1]) << std::endl;
+
+ vertices[0] = &points[0]; vertices[1] = &points[4]; vertices[2] = &points[5];
+ std::cout << "{0, 4, 5}:" << std::endl;
+ s = k.circumsphere(vertices);
+ std::cout << "Circumsphere: " << s.center() << " " << s.squared_radius() << std::endl;
+ std::cout << "Side of: " << k.side_of_circumsphere(vertices, *vertices[1]) << std::endl;
+
+ // Degenerate
+ vertices[0] = &points[0]; vertices[1] = &points[1]; vertices[2] = &points[6];
+ std::cout << "{0, 1, 6}:" << std::endl;
+ s = k.circumsphere(vertices);
+ std::cout << "Circumsphere: " << s.center() << " " << s.squared_radius() << std::endl;
+ std::cout << "Side of: " << k.side_of_circumsphere(vertices, *vertices[1]) << std::endl;
+ }
+
+ // Tetrahedron
+ {
+ PointContainer vertices(4);
+ vertices[0] = &points[3]; vertices[1] = &points[1]; vertices[2] = &points[2]; vertices[3] = &points[0];
+ std::cout << "{3, 1, 2, 0}:" << std::endl;
+ Sphere s = k.circumsphere(vertices);
+ std::cout << "Circumsphere: " << s.center() << " " << s.squared_radius() << std::endl;
+ std::cout << "Side of: " << k.side_of_circumsphere(vertices, *vertices[1]) << std::endl;
+ }
+
+ // Tetrahedron
+ {
+ PointContainer vertices(3);
+ vertices[0] = &points[3]; vertices[1] = &points[1]; vertices[2] = &points[2];
+ std::cout << "{3, 1, 2}:" << std::endl;
+ Sphere s = k.circumsphere(vertices);
+ std::cout << "Circumsphere: " << s.center() << ", radius: " << s.squared_radius() << std::endl;
+ std::cout << "Side of: " << k.side_of_circumsphere(vertices, points[0]) << std::endl;
+ }
+#endif
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/geometry/test-eventqueue.cpp Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,27 @@
+#include <eventqueue.h>
+#include <functional>
+#include <iostream>
+
+int main()
+{
+ typedef EventQueue<int, std::less<int> > EQ;
+ typedef EQ::iterator iterator;
+
+ EQ queue;
+
+ iterator i = queue.push(4);
+ queue.push(2);
+ queue.push(7);
+ iterator j = queue.push(6);
+ queue.push(5);
+
+ *i = 8;
+ queue.update(i);
+ queue.remove(j);
+
+ while (!queue.empty())
+ {
+ std::cout << *queue.top() << std::endl;
+ queue.pop();
+ }
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/geometry/test-kinetic-sort.cpp Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,66 @@
+#include <polynomial.h>
+#include <simulator.h>
+#include <kinetic-sort.h>
+#include <iostream>
+
+#include <boost/utility.hpp>
+
+//typedef double FieldType;
+//typedef ZZ FieldType;
+typedef QQ FieldType;
+typedef UPolynomial<FieldType> PolyKernel;
+typedef PolyKernel::Polynomial Polynomial;
+typedef std::list<Polynomial> SortDS;
+typedef Simulator<PolyKernel> SimulatorFT;
+
+class TrajectoryExtractor
+{
+ public:
+ Polynomial operator()(SortDS::iterator i) const { return *i; }
+};
+
+typedef KineticSort<SortDS, TrajectoryExtractor, SimulatorFT> KineticSortDS;
+
+struct EvaluatedComparison: public std::binary_function<const Polynomial&, const Polynomial&, bool>
+{
+ EvaluatedComparison(FieldType v): vv(v) {}
+ bool operator()(const Polynomial& p1, const Polynomial& p2)
+ { return p1(vv) < p2(vv); }
+ FieldType vv;
+};
+
+void swap(SortDS* s, SortDS::iterator i)
+{
+ std::cout << "Swapping " << *i << " " << *boost::next(i) << std::endl;
+ s->splice(i, *s, boost::next(i));
+}
+
+int main()
+{
+ SimulatorFT simulator;
+ SortDS list;
+
+ // Insert polynomials and sort the list for current time
+ list.push_back(Polynomial("x^3 - 3"));
+ list.push_back(Polynomial("x^2 - 2*x - 2"));
+ list.push_back(Polynomial("2*x - 4"));
+ //list.sort(EvaluatedComparison(simulator.current_time()));
+ list.sort(EvaluatedComparison(0));
+
+ // Print out the list
+ for (SortDS::const_iterator cur = list.begin(); cur != list.end(); ++cur)
+ std::cout << *cur << std::endl;
+
+ // Setup kinetic sort
+ KineticSortDS ks(&list, &simulator, swap);
+
+ while(!simulator.reached_infinity() && simulator.current_time() < 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;
+ }
+
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/geometry/test-linalg.cpp Sun Jul 01 00:00:00 2007 -0400
@@ -0,0 +1,32 @@
+#include "linalg.h"
+#include <iostream>
+
+#include <synaps/upol.h>
+#include <synaps/upol/gcd.h>
+#include "rational-function.h"
+
+
+typedef UPolDse<double> Polynomial;
+typedef RationalFunction<Polynomial> RationalF;
+
+//typedef LinearAlgebra<double> LinAlg;
+//typedef LinearAlgebra<Polynomial> LinAlg;
+typedef LinearAlgebra<RationalF> LinAlg;
+
+int main()
+{
+ LinAlg::MatrixType a(2,2);
+ a(0,0) = Polynomial("3*x"); a(1,0) = Polynomial(4); a(0,1) = Polynomial(0); a(1,1) = Polynomial(7);
+
+ std::cout << a << std::endl;
+ std::cout << LinAlg::determinant(a) << std::endl;
+ std::cout << Polynomial("3*x^2 + 4*x") / Polynomial("3*x^2") << std::endl;
+
+ LinAlg::VectorType b(2), x;
+ b(0) = Polynomial(4); b(1) = Polynomial(3);
+ LinAlg::solve(a, b, x);
+ std::cout << x << std::endl;
+
+ x(0).normalize(); x(1).normalize();
+ std::cout << x << std::endl;
+}