--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/89ae955518665a61 Wed Feb 27 07:57:45 2008 -0500
@@ -0,0 +1,16 @@
+From artemis Tue Feb 26 23:24:44 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Tue, 26 Feb 2008 18:22:06 -0500
+State: new
+Subject: Get rid of intostring() and .c_str()
+Message-Id: <89ae955518665a61-0-artemis@metatron>
+
+Get rid of the need for intostring() (in addition to tostring()), and having to
+place .c_str() after tostring() in rLog calls.
+
+The former is necessary because of some problem with disambiguating which
+operator<<(ostream,T) to use when Event is being output, so intostring() calls
+T.operator<<(ostream) explicitly. This problem seems to exist only for Events.
+
+It should be possible to solve the latter by returning char* from tostring()
+rather than std::string.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/b5d8a8403ae3a0d2 Wed Feb 27 07:57:45 2008 -0500
@@ -0,0 +1,18 @@
+From artemis Thu Feb 21 09:54:54 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Thu, 21 Feb 2008 04:54:00 -0500
+State: fixed
+Subject: Segfault in tests/geometry/test-kinetic-sort.cpp
+Message-Id: <b5d8a8403ae3a0d2-0-artemis@metatron>
+
+test-kinetic-sort segfaults when we subscribe to geometry/simulator RLog channel.
+
+From artemis Fri Feb 22 18:05:55 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Fri, 22 Feb 2008 13:05:55 -0500
+Subject: properties changes (state)
+Message-Id: <b5d8a8403ae3a0d2-a0288cca74095157-artemis@metatron>
+References: <b5d8a8403ae3a0d2-0-artemis@metatron>
+In-Reply-To: <b5d8a8403ae3a0d2-0-artemis@metatron>
+
+state=fixed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/c420501cc5285bbc Wed Feb 27 07:57:45 2008 -0500
@@ -0,0 +1,31 @@
+From artemis Tue Feb 26 10:25:07 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Tue, 26 Feb 2008 05:22:56 -0500
+State: fixed
+Subject: Non-optimized CGAL runtime error
+Message-Id: <c420501cc5285bbc-0-artemis@metatron>
+
+If the code is compiled with optimizations off or with debug on, CGAL gives a
+runtime error.
+
+From artemis Tue Feb 26 13:00:23 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Tue, 26 Feb 2008 08:00:23 -0500
+Subject: properties changes (state)
+Message-Id: <c420501cc5285bbc-e983173b6cd399a6-artemis@metatron>
+References: <c420501cc5285bbc-0-artemis@metatron>
+In-Reply-To: <c420501cc5285bbc-0-artemis@metatron>
+
+state=fixed
+
+From artemis Tue Feb 26 13:01:57 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Tue, 26 Feb 2008 08:00:50 -0500
+Subject: Fixed by adding CGAL_NO_ASSERTIONS
+Message-Id: <c420501cc5285bbc-d75d85d67d421f7e-artemis@metatron>
+References: <c420501cc5285bbc-0-artemis@metatron>
+In-Reply-To: <c420501cc5285bbc-0-artemis@metatron>
+
+I suspect the bug in CGAL (since everything works fine with CGAL_NO_ASSERTIONS
+set), so I disabled the offending assertion (and all the rest of them) by
+setting a CXX flag in CMakeLists.txt.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/cd79223a108d3900 Wed Feb 27 07:57:45 2008 -0500
@@ -0,0 +1,9 @@
+From artemis Fri Feb 22 23:41:38 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Fri, 22 Feb 2008 18:39:31 -0500
+State: new
+Subject: Incorrect sign_at() in UPolynomial<double>
+Message-Id: <cd79223a108d3900-0-artemis@metatron>
+
+UPolynomial<double>::sign_at() reports -1 as the sign of -2*x + 4 at x = 2.
+Should be 0.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/ebda8db3f9908e33 Wed Feb 27 07:57:45 2008 -0500
@@ -0,0 +1,8 @@
+From artemis Mon Feb 25 16:30:10 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Mon, 25 Feb 2008 11:29:27 -0500
+State: new
+Subject: Efficient EventQueue
+Message-Id: <ebda8db3f9908e33-0-artemis@metatron>
+
+Change EventQueue to an efficient implementation, e.g., using a Fibonacci heap.
--- a/CMakeLists.txt Mon Feb 25 16:54:34 2008 -0500
+++ b/CMakeLists.txt Wed Feb 27 07:57:45 2008 -0500
@@ -39,6 +39,15 @@
${gmp_LIBRARY}
${gmpxx_LIBRARY}
${m_LIBRARY})
+set (cgal_cxxflags ${cgal_cxxflags}
+ -DCGAL_NO_ASSERTIONS)
+
+# 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)
--- a/examples/ar-vineyard/CMakeLists.txt Mon Feb 25 16:54:34 2008 -0500
+++ b/examples/ar-vineyard/CMakeLists.txt Wed Feb 27 07:57:45 2008 -0500
@@ -1,8 +1,12 @@
set (targets
ar-vineyard)
-
+
+find_library (boost_signals_LIBRARY NAMES boost_signals
+ PATHS ${Boost_LIBRARY_DIR})
+
add_definitions (${cgal_cxxflags})
foreach (t ${targets})
- add_executable (${t} ${t}.cpp)
- target_link_libraries (${t} ${libraries} ${cgal_libraries})
+ add_executable (${t} ${t}.cpp ${external_sources})
+ target_link_libraries (${t} ${libraries} ${cgal_libraries}
+ ${synaps_libraries} ${boost_signals_LIBRARY})
endforeach (t)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/ar-vineyard/ar-function-kernel.h Wed Feb 27 07:57:45 2008 -0500
@@ -0,0 +1,80 @@
+#ifndef __AR_FUNCTION_KERNEL_H__
+#define __AR_FUNCTION_KERNEL_H__
+
+
+#include <stack>
+#include <iostream>
+#include "ar-simplex3d.h"
+
+
+
+/**
+ * Represents function suitable for the FunctionKernel. Albeit very restrictive
+ * (it only supports subtraction), although that should suffice for KineticSort.
+ */
+class ARFunction
+{
+ public:
+ /// Represents the three forms of the involved functions. See Appendix B of LHI paper.
+ enum FunctionForm { none, rho, lambda, phi };
+
+ public:
+ ARFunction(FunctionForm f, ARSimplex3D* s):
+ form_(f), form2_(none),
+ simplex_(s), simplex2_(0) {}
+ ARFunction(const ARFunction& other):
+ form_(other.form_), form2_(other.form2_),
+ simplex_(other.simplex_),
+ simplex2_(other.simplex2_) {}
+
+ ARFunction& operator-=(const ARFunction& other) { form2_ = other.form_; simplex2_ = other.simplex_; return *this; }
+ ARFunction operator-(const ARFunction& other) { return (ARFunction(*this) -= other); }
+
+ FunctionForm form() const { return form_; }
+ FunctionForm form2() const { return form2_; }
+ ARSimplex3D* simplex() const { return simplex_; }
+ ARSimplex3D* simplex2() const { return simplex2_; }
+
+ std::ostream& operator<<(std::ostream& out) const { return (out << "(" << form_ << ", " << simplex_ << "), ("
+ << form2_ << ", " << simplex2_ << ")"); }
+
+ private:
+ FunctionForm form_, form2_; // the function is form_ - form2_
+ ARSimplex3D *simplex_, *simplex2_;
+};
+
+std::ostream&
+operator<<(std::ostream& out, const ARFunction& f)
+{ return f.operator<<(out); }
+
+/**
+ * Function kernel specialized at solving the kinds of functions involved in
+ * ARVineyard construction. We cannot use a polynomial kernel (which includes
+ * rational functions) that we already have because not all involved functions
+ * are rational.
+ */
+class ARFunctionKernel
+{
+ public:
+ typedef double FieldType;
+ typedef FieldType RootType;
+ typedef std::stack<RootType> RootStack;
+ typedef ARSimplex3D::RealValue SimplexFieldType;
+
+ typedef ARFunction Function;
+ typedef Function::FunctionForm FunctionForm;
+
+ public:
+ static void solve(const Function& f, RootStack& stack);
+ static RootType root(RootType r) { return r; }
+ static RootType root(SimplexFieldType r) { return CGAL::to_double(r); }
+ static int sign_at(const Function& f, RootType r);
+ static RootType between(RootType r1, RootType r2) { return (r1 + r2)/2; }
+ static int sign_at_negative_infinity(const Function& f);
+
+ static FieldType value_at(const Function& f, RootType v);
+};
+
+#include "ar-function-kernel.hpp"
+
+#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/ar-vineyard/ar-function-kernel.hpp Wed Feb 27 07:57:45 2008 -0500
@@ -0,0 +1,144 @@
+#include <utilities/log.h>
+#include <cmath>
+
+#if LOGGING
+static rlog::RLogChannel* rlARFunctionKernel = DEF_CHANNEL("ar/function-kernel", rlog::Log_Debug);
+static rlog::RLogChannel* rlARFunctionKernelValue = DEF_CHANNEL("ar/function-kernel/value", rlog::Log_Debug);
+#endif
+
+
+/* For now only rho and phi are implemented */
+void
+ARFunctionKernel::
+solve(const Function& f, RootStack& stack)
+{
+ AssertMsg(stack.empty(), "Stack must be empty before solving");
+ AssertMsg((f.form() != Function::none) && (f.form2() != Function::none), "Solving possible only for differences");
+
+ FunctionForm f1 = f.form(), f2 = f.form2();
+ const ARSimplex3D *s1 = f.simplex(),
+ *s2 = f.simplex2();
+ if (f1 < f2) { std::swap(f1,f2); std::swap(s1,s2); } // for simplicity enforce function order
+ // to handle fewer cases explicitly
+ AssertMsg(f1 != Function::lambda && f2 != Function::lambda, "Lambda not implemented yet");
+
+ rLog(rlARFunctionKernel, "Solve: function1 = (%i, %p), function2 = (%i, %p)",
+ f1, s1, f2, s2);
+
+ //if (f1 == Function::phi && f2 == Function::phi) return;
+ //if (f1 == Function::rho && f2 == Function::rho) return;
+
+ if (f1 == Function::phi && f2 == Function::rho)
+ {
+ SimplexFieldType r = s2->alpha() - s1->phi_const();
+ rLog(rlARFunctionKernel, " phi = rho => r^2 = %s (%lf)", tostring(r).c_str(), root(r));
+ stack.push(root(r));
+ }
+
+ if (f1 == Function::phi && f2 == Function::lambda)
+ {
+ rLog(rlARFunctionKernel, " phi = lambda");
+ SimplexFieldType r2 = (s2->rho() + s2->v() - s2->s() - s1->phi_const());
+ r2 *= r2;
+ r2 /= 4*s2->v();
+ r2 += s2->s();
+ if (r2 >= s2->s() + s2->v())
+ stack.push(root(r2));
+
+ SimplexFieldType r1 = s2->rho() - s1->phi_const();
+ if (r1 <= s2->s() + s2->v())
+ stack.push(root(r1));
+ }
+
+ // FIXME: this is far from complete!
+ if (f1 == Function::lambda && f2 == Function::lambda)
+ {
+ rLog(rlARFunctionKernel, " lambda = lambda");
+ if ((s1->s() + s1->v() < s2->s() + s2->v())) // let f1 be the function with larger break point
+ { std::swap(f1,f2); std::swap(s1,s2); }
+
+ if (s1->rho() > s2->rho())
+ {
+ RootType r = root(s2->s() + s2->v() + s1->rho() - s2->rho()) + 2*sqrt(root(s2->v()*(s1->rho() - s2->rho())));
+ if (r < root(s1->s() + s1->v()) && r > root(s2->s() + s2->v()))
+ stack.push(r);
+ }
+ }
+
+ if (f1 == Function::lambda && f2 == Function::rho)
+ {
+ rLog(rlARFunctionKernel, " lambda = rho");
+ // perhaps no solutions instead of an assertion is the right way to deal with this
+ AssertMsg(s2->alpha() > s1->rho(), "Rho_0^2 must be greater than Rho^2");
+
+ RootType r = sqrt(root(s2->v()*(s2->alpha() - s1->rho()))); // damn square roots
+ r *= 2;
+ r += root(s1->s() + s1->v() + s2->alpha() - s1->rho());
+ }
+ rLog(rlARFunctionKernel, " Stack size: %i", stack.size());
+ if (stack.size() > 0) rLog(rlARFunctionKernel, " Top: %lf", stack.top());
+}
+
+int
+ARFunctionKernel::
+sign_at(const Function& f, RootType r)
+{
+ FieldType v = value_at(f,r);
+ if (v > 0) return true;
+ else return false;
+}
+
+int
+ARFunctionKernel::
+sign_at_negative_infinity(const Function& f)
+{
+ FunctionForm f1 = f.form(), f2 = f.form2();
+ const ARSimplex3D *s1 = f.simplex(),
+ *s2 = f.simplex2();
+ int multiplier = 1;
+ if (f1 < f2) { std::swap(f1, f2); std::swap(s1, s2); multiplier = -1; }
+
+ AssertMsg(f1 != Function::lambda && f2 != Function::lambda, "Lambda not implemented yet");
+
+ if (f1 == Function::phi && f2 == Function::phi)
+ {
+ if (s1->phi_const() == s2->phi_const()) return 0;
+ if (s1->phi_const() > s2->phi_const()) return 1; // multiplier must be 1
+ else return -1;
+ }
+
+ if (f1 == Function::phi && f2 == Function::rho)
+ return -multiplier;
+
+ if (f1 == Function::rho && f2 == Function::rho)
+ {
+ if (s1->alpha() == s2->alpha()) return 0;
+ if (s1->alpha() > s2->alpha()) return 1; // multiplier must be 1
+ else return -1;
+ }
+
+ AssertMsg(false, "The case analysis should be exhaustive in sign at -infinity");
+ return false;
+}
+
+ARFunctionKernel::FieldType
+ARFunctionKernel::
+value_at(const Function& f, RootType v)
+{
+ FunctionForm f1 = f.form(), f2 = f.form2();
+ ARSimplex3D *s1 = f.simplex(),
+ *s2 = f.simplex2();
+
+ AssertMsg(f2 == Function::none && s2 == 0, "Value_at knows only about functions themselves, not their differences");
+ AssertMsg(f1 != Function::lambda, "Lambda not implemented yet");
+ rLog(rlARFunctionKernelValue, "Value_at: function = (%i, %p)", f1, s1);
+
+ if (f1 == Function::phi)
+ return v + root(s1->phi_const());
+
+ if (f1 == Function::rho)
+ return root(s1->alpha());
+
+ AssertMsg(false, "The case analysis should be exhaustive in value_at");
+ return 0;
+}
--- a/examples/ar-vineyard/ar-simplex3d.h Mon Feb 25 16:54:34 2008 -0500
+++ b/examples/ar-vineyard/ar-simplex3d.h Wed Feb 27 07:57:45 2008 -0500
@@ -28,7 +28,6 @@
typedef Delaunay::Facet Facet;
typedef Delaunay::Cell Cell;
typedef Delaunay::Cell_handle Cell_handle;
-typedef K::FT RealValue;
typedef Delaunay::Finite_vertices_iterator Vertex_iterator;
typedef Delaunay::Finite_edges_iterator Edge_iterator;
@@ -40,6 +39,7 @@
class ARSimplex3D: public SimplexWithVertices<Vertex_handle>
{
public:
+ typedef K::FT RealValue;
typedef std::map<ARSimplex3D, RealValue> SimplexPhiMap;
typedef SimplexWithVertices<Vertex_handle> Parent;
typedef Parent::VertexContainer VertexSet;
@@ -86,9 +86,11 @@
RealValue rho_; // rho_ is the squared radius of the smallest circumsphere
RealValue s_; // s_ is the squared distance from z to the affine hull of the simplex
RealValue v_; // v_ is the squared distance from z to the affine hull of the dual Voronoi cell
- RealValue phi_const_; // see LHI paper, Appendix 1
+ RealValue phi_const_; // see LHI paper, Appendices A and B
bool attached_;
+ // in paper's notation: s_ = v^2; v_ = d^2
+
};
typedef std::vector<ARSimplex3D> ARSimplex3DVector;
--- a/examples/ar-vineyard/ar-simplex3d.hpp Mon Feb 25 16:54:34 2008 -0500
+++ b/examples/ar-vineyard/ar-simplex3d.hpp Wed Feb 27 07:57:45 2008 -0500
@@ -1,3 +1,9 @@
+#include <utilities/log.h>
+
+#if LOGGING
+static rlog::RLogChannel* rlARSimplex3D = DEF_CHANNEL("ar/simplex3d", rlog::Log_Debug);
+#endif
+
ARSimplex3D::
ARSimplex3D(const ARSimplex3D& s): Parent(s)
{
@@ -209,7 +215,7 @@
operator<<(std::ostream& out) const
{
for (VertexSet::const_iterator cur = Parent::vertices().begin(); cur != Parent::vertices().end(); ++cur)
- out << **cur << ", ";
+ out << &(**cur) << ", ";
out << "value = " << value();
return out;
@@ -226,19 +232,19 @@
ARSimplex3D::SimplexPhiMap simplices;
for(Cell_iterator cur = Dt.finite_cells_begin(); cur != Dt.finite_cells_end(); ++cur)
update_simplex_phi_map(ARSimplex3D(*cur, z), simplices);
- std::cout << "Cells inserted" << std::endl;
+ rLog(rlARSimplex3D, "Cells inserted");
for(Vertex_iterator cur = Dt.finite_vertices_begin(); cur != Dt.finite_vertices_end(); ++cur)
simplices[ARSimplex3D(*cur, z)] = 0; // only one tetrahedron can have non-negative phi_const value
// (namely the one containing z); all other simplices will have a
// negative phi_const value, so 0 is safe
- std::cout << "Vertices inserted" << std::endl;
+ rLog(rlARSimplex3D, "Vertices inserted");
for(Facet_iterator cur = Dt.finite_facets_begin(); cur != Dt.finite_facets_end(); ++cur)
update_simplex_phi_map(ARSimplex3D(*cur, z, simplices), simplices);
- std::cout << "Facets inserted" << std::endl;
+ rLog(rlARSimplex3D, "Facets inserted");
for(Edge_iterator cur = Dt.finite_edges_begin(); cur != Dt.finite_edges_end(); ++cur)
update_simplex_phi_map(ARSimplex3D(*cur, z, simplices, Dt.incident_facets(*cur)), simplices);
- std::cout << "Edges inserted" << std::endl;
+ rLog(rlARSimplex3D, "Edges inserted");
// Sort simplices by their alpha values
alpha_order.resize(simplices.size()); ARSimplex3DVector::iterator out = alpha_order.begin();
--- a/examples/ar-vineyard/ar-vineyard.cpp Mon Feb 25 16:54:34 2008 -0500
+++ b/examples/ar-vineyard/ar-vineyard.cpp Wed Feb 27 07:57:45 2008 -0500
@@ -11,7 +11,11 @@
#ifdef LOGGING
rlog::RLogInit(argc, argv);
- stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
+ stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+ //stdoutLog.subscribeTo( RLOG_CHANNEL("ar/vineyard") );
+ //stdoutLog.subscribeTo( RLOG_CHANNEL("ar/function-kernel/value") );
+ //stdoutLog.subscribeTo( RLOG_CHANNEL("geometry/simulator") );
+
//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/filtration") );
//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/cycle") );
//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/vineyard") );
@@ -52,7 +56,7 @@
arv.compute_pairing();
// Compute vineyard
- arv.compute_vineyard(true);
+ arv.compute_vineyard();
std::cout << "Vineyard computed" << std::endl;
arv.vineyard()->save_edges(outfilename);
}
--- a/examples/ar-vineyard/ar-vineyard.h Mon Feb 25 16:54:34 2008 -0500
+++ b/examples/ar-vineyard/ar-vineyard.h Wed Feb 27 07:57:45 2008 -0500
@@ -6,70 +6,109 @@
#ifndef __AR_VINEYARD_H__
#define __AR_VINEYARD_H__
+#include <boost/signals.hpp>
+#include <boost/bind.hpp>
+#include <list>
+#include <vector>
+
#include "topology/conesimplex.h"
#include "topology/filtration.h"
+#include "geometry/kinetic-sort.h"
+#include "geometry/simulator.h"
-#include <CGAL/Kinetic/Inexact_simulation_traits.h>
-#include <CGAL/Kinetic/Event_base.h>
-#include <CGAL/Kinetic/Sort.h>
-#include <CGAL/Kinetic/Sort_visitor_base.h>
-
-#include <list>
#include "ar-simplex3d.h"
-
-#include <vector>
+#include "ar-function-kernel.h"
-class ARVineyardBase
-{
- public:
- /// \name CGAL Kinetic Sort types
- /// @{
- class SortVisitor;
- typedef CGAL::Kinetic::Inexact_simulation_traits Traits;
- typedef CGAL::Kinetic::Sort<Traits, SortVisitor> Sort;
- typedef Traits::Simulator Simulator;
- typedef Traits::Active_points_1_table ActivePointsTable;
- typedef ActivePointsTable::Key Key;
-
- typedef Traits::Kinetic_kernel::
- Function_kernel::Construct_function CF;
- typedef Traits::Kinetic_kernel::Motion_function F;
- /// @}
-
- class ARConeSimplex;
- class MembershipFunctionChangeEvent;
-};
-
-class ARVineyardBase::ARConeSimplex: public ConeSimplex<ARSimplex3D>
+template <class Simulator_>
+class ARConeSimplex3D: public ConeSimplex<ARSimplex3D>
{
public:
typedef ConeSimplex<ARSimplex3D> Parent;
typedef ARSimplex3D ARSimplex3D;
+
+ /// \name Simulator types
+ /// @{
+ typedef Simulator_ Simulator;
+ typedef typename Simulator::FunctionKernel FunctionKernel;
+ typedef typename FunctionKernel::Function Function;
+ /// @}
+
+ /// \name ThresholdSort types
+ /// @{
+ typedef std::list<Function> ThresholdList;
+ typedef typename ThresholdList::iterator ThresholdListIterator;
- ARConeSimplex(const ARSimplex3D& s, bool coned = false):
- Parent(s, coned) {}
+ struct ThresholdTrajectoryExtractor
+ { Function operator()(ThresholdListIterator i) const { return *i; } };
+
+ typedef KineticSort<ThresholdListIterator,
+ ThresholdTrajectoryExtractor, Simulator> ThresholdSort;
+ /// @}
- Key kinetic_key() const { return key_; }
- void set_kinetic_key(Key k) { key_ = k; }
+ typedef boost::signal<void (Simulator*)> NewMaxSignal;
+
+ public:
+ ARConeSimplex3D(const ARSimplex3D& s, bool coned = false);
+ ARConeSimplex3D(const ARConeSimplex3D& other): // need explicit copy-constructor because of the signal
+ Parent(other, other.coned()),
+ thresholds_(other.thresholds_) {}
+
+ const ThresholdList& thresholds() const { return thresholds_; }
+
+ NewMaxSignal& new_max_signal() { return new_max_signal_; }
+ const Function& max_threshold() const { return thresholds_.back(); }
+ void schedule_thresholds(Simulator* simulator);
+
+ // need explicit operator= because of the signal
+ ARConeSimplex3D& operator=(const ARConeSimplex3D& other) { Parent::operator=(other); thresholds_ = other.thresholds_; return *this; }
+ bool operator<(const ARConeSimplex3D& other) const { if (coned() ^ other.coned()) return !coned(); else return Parent::operator<(other); }
+
private:
- Key key_;
+ ThresholdList thresholds_;
+ ThresholdSort thresholds_sort_;
+ NewMaxSignal new_max_signal_;
+
+ void swap_thresholds(ThresholdListIterator i, Simulator* simulator);
};
-
-class ARVineyard: public ARVineyardBase
+/**
+ * Encapsulated filtration, and provides compute_vineyard() functionality.
+ */
+class ARVineyard
{
public:
typedef ARVineyard Self;
+
+ /// \name FunctionKernel and Simulator types
+ /// @{
+ typedef ARFunctionKernel FunctionKernel;
+ typedef FunctionKernel::Function Function;
+ typedef Simulator<FunctionKernel> Simulator;
+ /// @}
+
+ /// \name Filtration types
+ /// @{
+ typedef ARConeSimplex3D<Simulator> ARConeSimplex3D;
+ typedef Filtration<ARConeSimplex3D> Filtration;
+ typedef Filtration::Simplex Simplex;
+ typedef Filtration::Index Index;
+ typedef Filtration::Vineyard Vineyard;
+ typedef Vineyard::Evaluator Evaluator;
+ /// @}
- typedef Filtration<ARConeSimplex> ARFiltration;
- typedef ARFiltration::Simplex Simplex;
- typedef ARFiltration::Index Index;
- typedef ARFiltration::Vineyard Vineyard;
- typedef Vineyard::Evaluator Evaluator;
- typedef std::map<Key, Index> KeyIndexMap;
+ /// \name SimplexSort types
+ /// @{
+ struct SimplexTrajectoryExtractor
+ { Function operator()(Index i) const { return i->max_threshold(); } };
+
+ typedef KineticSort<Index, SimplexTrajectoryExtractor, Simulator> SimplexSort;
+ typedef SimplexSort::iterator SimplexSortIterator;
+ class ThresholdChangeSlot; // used to notify of change in max threshold
+ /// @}
+
typedef std::list<Point> PointList;
class StaticEvaluator;
@@ -80,26 +119,23 @@
~ARVineyard();
void compute_pairing();
- void compute_vineyard(bool explicit_events = false);
+ void compute_vineyard();
- const ARFiltration* filtration() const { return filtration_; }
+ const Filtration* filtration() const { return filtration_; }
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();
void change_evaluator(Evaluator* eval);
private:
- ARFiltration* filtration_;
+ Filtration* filtration_;
Vineyard* vineyard_;
Evaluator* evaluator_;
- KeyIndexMap kinetic_map_;
-
Point z_;
Delaunay dt_;
@@ -120,68 +156,38 @@
//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_;
+class ARVineyard::ThresholdChangeSlot
+{
+ public:
+ ThresholdChangeSlot(SimplexSortIterator iter, SimplexSort* sort):
+ iter_(iter), sort_(sort) { iter_->element->new_max_signal().connect(*this); }
+ void operator()(Simulator* simulator) { sort_->update_trajectory(iter_, simulator); }
+
+ private:
+ SimplexSortIterator iter_;
+ SimplexSort* sort_;
};
-std::ostream& operator<<(std::ostream& out, const ARVineyardBase::MembershipFunctionChangeEvent& e)
-{ return e.operator<<(out); }
-
class ARVineyard::StaticEvaluator: public Evaluator
{
public:
- StaticEvaluator(RealType t): time_(t) {}
+ StaticEvaluator() {}
- virtual RealType time() const { return time_; }
+ virtual RealType time() const { return 0; }
virtual RealType value(const Simplex& s) const { return s.value(); }
-
- private:
- RealType time_;
};
class ARVineyard::KineticEvaluator: public Evaluator
{
public:
- KineticEvaluator(Simulator::Handle sp,
- ActivePointsTable::Handle apt,
- RealType time_offset):
- sp_(sp), apt_(apt) {}
+ KineticEvaluator(Simulator* simulator):
+ simulator_(simulator) {}
- virtual RealType time() const { return CGAL::to_double(get_time()); }
- virtual RealType value(const Simplex& s) const { return CGAL::to_double(apt_->at(s.kinetic_key()).x()(get_time())); }
+ virtual RealType time() const { return simulator_->current_time(); }
+ virtual RealType value(const Simplex& s) const { return FunctionKernel::value_at(s.max_threshold(), time()); }
private:
- Simulator::Time get_time() const { return sp_->current_time(); }
-
- Simulator::Handle sp_;
- ActivePointsTable::Handle apt_;
-};
-
-
-class ARVineyardBase::SortVisitor: public CGAL::Kinetic::Sort_visitor_base
-{
- public:
- SortVisitor(ARVineyard* arv): arv_(arv) {}
-
- template<class Vertex_handle>
- void before_swap(Vertex_handle a, Vertex_handle b) const;
-
- private:
- ARVineyard* arv_;
+ Simulator* simulator_;
};
--- a/examples/ar-vineyard/ar-vineyard.hpp Mon Feb 25 16:54:34 2008 -0500
+++ b/examples/ar-vineyard/ar-vineyard.hpp Wed Feb 27 07:57:45 2008 -0500
@@ -1,26 +1,69 @@
#include <utilities/log.h>
/* Implementation */
+
+#ifdef LOGGING
+static rlog::RLogChannel* rlARVineyard = DEF_CHANNEL("ar/vineyard", rlog::Log_Debug);
+static rlog::RLogChannel* rlARVineyardComputing = DEF_CHANNEL("ar/vineyard/computing", rlog::Log_Debug);
+static rlog::RLogChannel* rlARVineyardSwap = DEF_CHANNEL("ar/vineyard/swap", rlog::Log_Debug);
+static rlog::RLogChannel* rlARVineyardConeSimplexSwap = DEF_CHANNEL("ar/vineyard/conesimplex/swap", rlog::Log_Debug);
+#endif
+
+
+template <class Simulator_>
+ARConeSimplex3D<Simulator_>::
+ARConeSimplex3D(const ARSimplex3D& s, bool coned): Parent(s, coned)
+{}
+template <class Simulator_>
+void
+ARConeSimplex3D<Simulator_>::
+swap_thresholds(ThresholdListIterator i, Simulator* simulator)
+{
+ rLog(rlARVineyardConeSimplexSwap, "Transposing %s and %s", tostring(*i).c_str(),
+ tostring(*boost::next(i)).c_str());
+ typename ThresholdList::iterator n = boost::next(i);
+ thresholds_.splice(i, thresholds_, n);
+ if (boost::next(i) == thresholds_.end())
+ new_max_signal_(simulator);
+}
+
+template <class Simulator_>
+void
+ARConeSimplex3D<Simulator_>::
+schedule_thresholds(Simulator* simulator)
+{
+ if (!coned()) thresholds_.push_back(Function(Function::rho, this));
+ else
+ {
+ thresholds_.push_back(Function(Function::phi, this));
+ thresholds_.push_back(Function(Function::rho, this));
+ }
+
+ thresholds_sort_.initialize(thresholds_.begin(), thresholds_.end(),
+ boost::bind(&ARConeSimplex3D::swap_thresholds, this, _1, _2), simulator);
+}
+
+
ARVineyard::
ARVineyard(const PointList& points, const Point& z): z_(z)
{
for (PointList::const_iterator cur = points.begin(); cur != points.end(); ++cur)
dt_.insert(*cur);
- std::cout << "Delaunay triangulation computed" << std::endl;
+ rLog(rlARVineyard, "Delaunay triangulation computed");
ARSimplex3DVector alpha_ordering;
fill_alpha_order(dt_, z_, alpha_ordering);
- std::cout << "Delaunay simplices: " << alpha_ordering.size() << std::endl;
+ rLog(rlARVineyard, "Delaunay simplices: %i", alpha_ordering.size());
- evaluator_ = new StaticEvaluator(0);
+ evaluator_ = new StaticEvaluator;
vineyard_ = new Vineyard(evaluator_);
- filtration_ = new ARFiltration(vineyard_);
+ filtration_ = new Filtration(vineyard_);
for (ARSimplex3DVector::const_iterator cur = alpha_ordering.begin(); cur != alpha_ordering.end(); ++cur)
{
- filtration_->append(*cur); // Delaunay simplex
- filtration_->append(ARConeSimplex(*cur, true)); // Coned off delaunay simplex
+ filtration_->append(ARConeSimplex3D(*cur)); // Delaunay simplex
+ filtration_->append(ARConeSimplex3D(*cur, true)); // Coned off delaunay simplex
}
}
@@ -39,98 +82,53 @@
filtration_->fill_simplex_index_map();
filtration_->pair_simplices(filtration_->begin(), filtration_->end());
vineyard_->start_vines(filtration_->begin(), filtration_->end());
- std::cout << "Simplices paired" << std::endl;
+ rLog(rlARVineyard, "Simplices paired");
}
void
ARVineyard::
-compute_vineyard(bool explicit_events)
+compute_vineyard()
{
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();
+ // Schedule thresholds
for (Index cur = filtration_->begin(); cur != filtration_->end(); ++cur)
- {
- F x = cf(F::NT(CGAL::to_double(cur->alpha())));
- Point_1 p(x);
- cur->set_kinetic_key(apt->insert(p));
- kinetic_map_[cur->kinetic_key()] = cur;
-
- if (!cur->coned()) continue; // non-coned simplices stay put, so we are done
+ cur->schedule_thresholds(&simulator);
- Time lambda_alpha = CGAL::to_double((cur->alpha() - cur->rho())); // when lambda becomes greater than alpha
- lambda_alpha += 2*CGAL::sqrt(CGAL::to_double(cur->s()*lambda_alpha));
- lambda_alpha += CGAL::to_double(cur->s() + cur->v());
-
- Time phi_alpha = CGAL::to_double(cur->alpha() - cur->phi_const());
-
- Time phi_lambda = CGAL::to_double(cur->rho() + cur->s() - cur->v() - cur->phi_const());
- phi_lambda *= phi_lambda;
- phi_lambda /= CGAL::to_double(4*cur->s());
- phi_lambda += CGAL::to_double(cur->v());
+ // Once thresholds are scheduled, we can initialize the simplex_sort
+ simplex_sort.initialize(filtration_->begin(), filtration_->end(),
+ boost::bind(&ARVineyard::swap, this, _1, _2), &simulator);
+ rLog(rlARVineyardComputing, "SimplexSort initialized");
- Time sv = CGAL::to_double(cur->s() + cur->v());
-
- if (true || phi_lambda < sv || phi_lambda < phi_alpha) // FIXME: remove true
- {
- sp->new_event(Time(phi_alpha),
- MembershipFunctionChangeEvent(cur->kinetic_key(),
- cf(F::NT(CGAL::to_double(cur->phi_const())), 1),
- apt)); // \phi^2 = r^2 + \phi_c^2
- std::cout << "Scheduled" << std::endl;
- } else
- std::cout << "Not scheduled" << std::endl;
-
-
- //sp->new_event(Time(...), MembershipFunctionChangeEvent(cur->kinetic_key()));
-
- std::cout << *cur << std::endl;
- std::cout << "lambda_alpha: " << lambda_alpha << std::endl;
- std::cout << "phi_alpha: " << phi_alpha << std::endl;
- std::cout << "phi_lambda: " << phi_lambda << std::endl;
- std::cout << "s^2 + v^2: " << sv << std::endl;
- std::cout << std::endl;
- }
+ // Connect signals and slots
+ std::vector<ThresholdChangeSlot> slots;
+ slots.reserve(filtration_->size());
+ for (SimplexSortIterator cur = simplex_sort.begin(); cur != simplex_sort.end(); ++cur)
+ slots.push_back(ThresholdChangeSlot(cur, &simplex_sort));
+ rLog(rlARVineyardComputing, "Signals and slots connected");
- // Process all the events (compute the vineyard in the process)
- // FIXME: the time should not be 1, but something like twice the radius of
- // the pointset as seen from z
- change_evaluator(new KineticEvaluator(sp, apt, 0));
- if (explicit_events)
- {
- while (sp->next_event_time() < 1)
- {
- std::cout << "Next event time: " << sp->next_event_time() << std::endl;
- sp->set_current_event_number(sp->current_event_number() + 1);
- std::cout << "Processed event" << std::endl;
- }
- } else
- sp->set_current_time(1.0);
- std::cout << "Processed " << sp->current_event_number() << " events" << std::endl;
+ // Simulate
+ change_evaluator(new KineticEvaluator(&simulator));
+ while(!simulator.reached_infinity())
+ {
+ rLog(rlARVineyardComputing, "Current time before: %lf", simulator.current_time());
+ simulator.process();
+ }
- //change_evaluator(new StaticEvaluator(1));
vineyard_->record_diagram(filtration_->begin(), filtration_->end());
}
void
ARVineyard::
-swap(Key a, Key b)
+swap(Index i, Simulator* simulator)
{
- Index ao = kinetic_map_[a], bo = kinetic_map_[b];
- AssertMsg(filtration_->get_trails_cmp()(ao, bo), "In swap(a,b), a must precede b");
- filtration_->transpose(ao);
- AssertMsg(filtration_->get_trails_cmp()(bo, ao), "In swap(a,b), b must precede a after the transposition");
+ rLog(rlARVineyardSwap, "Transposing %p and %p: %s and %s",
+ &(*i), &(*boost::next(i)),
+ tostring(*i).c_str(), tostring(*boost::next(i)).c_str());
+ filtration_->transpose(i);
}
void
@@ -143,30 +141,3 @@
evaluator_ = eval;
vineyard_->set_evaluator(evaluator_);
}
-
-void
-ARVineyardBase::MembershipFunctionChangeEvent::
-process() const
-{
- apt_->set(key_, function_);
- std::cout << "Updated for phi's dominance" << std::endl;
-}
-
-
-template<class Vertex_handle>
-void
-ARVineyardBase::SortVisitor::
-before_swap(Vertex_handle a, Vertex_handle b) const
-{
- std::cout << "Swapping elements" << *a << " and " << *b << std::endl;
- arv_->swap(*a,*b);
-}
-
-
-std::ostream&
-ARVineyardBase::MembershipFunctionChangeEvent::
-operator<<(std::ostream& out) const
-{
- return out << "Membership change" << std::endl;
-}
-
--- a/include/dionysus.h Mon Feb 25 16:54:34 2008 -0500
+++ b/include/dionysus.h Wed Feb 27 07:57:45 2008 -0500
@@ -1,6 +1,6 @@
/*
* Author: Dmitriy Morozov
- * Department of Computer Science, Duke University, 2006
+ * Department of Computer Science, Duke University, 2006--2008
*
* For now file exists only to store the main page of the documentation
*/
@@ -9,3 +9,16 @@
* Detailed description of Dionysus with references to all the classes,
* and code samples goes here.
*/
+
+/**
+ * \defgroup topology Topology Classes
+ */
+
+/**
+ * \defgroup geometry Geometry Classes
+ */
+
+/**
+ * \defgroup kinetic Kinetic Data Structures Classes
+ * \ingroup geometry
+ */
--- a/include/geometry/euclidean.h Mon Feb 25 16:54:34 2008 -0500
+++ b/include/geometry/euclidean.h Wed Feb 27 07:57:45 2008 -0500
@@ -14,6 +14,10 @@
#include "number-traits.h"
+/**
+ * Geometric Kernel. Defines operations on geometric primitives.
+ * \ingroup geometry
+ */
template<class NumberType_ = double>
class Kernel
{
@@ -56,7 +60,10 @@
};
-/* Point */
+/**
+ * Point class.
+ * \ingroup geometry
+ */
template<class NumberType_>
class Kernel<NumberType_>::Point: public VectorType
{
@@ -77,7 +84,10 @@
};
-/* Sphere */
+/**
+ * Sphere class.
+ * \ingroup geometry
+ */
template<class NumberType_>
class Kernel<NumberType_>::Sphere
{
--- a/include/geometry/kinetic-sort.h Mon Feb 25 16:54:34 2008 -0500
+++ b/include/geometry/kinetic-sort.h Wed Feb 27 07:57:45 2008 -0500
@@ -10,33 +10,37 @@
* Maintains elements of the given data structure in the sorted order assuming the elements follow
* trajectories given by TrajectoryExtractor_.
*
- * \arg SortDS_ should be forward and backward iterable, swaps are handles via SwapCallback
- * \arg TrajectoryExtractor_ applied to the iterator into SortDS_ should return a rational
- * function describing the
+ * \arg ElementIterator_ iterator over the underlying data structure that's kept in sorted order
+ * \arg TrajectoryExtractor_ applied to the iterator into SortDS_ should return a function
+ * (of type Simulator_::FunctionKernel::Function) describing the trajectory of the element
* \arg Simulator_ the Simulator type, e.g. Simulator. Note that KineticSort does not store
* a pointer to the Simulator (so a pointer is passed in each relevant operation)
+ * \arg Swap_ is called with an ElementIterator_ when a swap needs to be performed
+ *
+ * \ingroup kinetic
*/
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+template<class ElementIterator_, class TrajectoryExtractor_,
+ class Simulator_, class Swap_ = boost::function<void(ElementIterator_ pos, Simulator_* simulator)> >
class KineticSort
{
public:
typedef Simulator_ Simulator;
- typedef typename Simulator::PolynomialKernel PolynomialKernel;
- typedef SortDS_ SortDS;
+ typedef typename Simulator::FunctionKernel FunctionKernel;
+ typedef ElementIterator_ ElementIterator;
+ typedef Swap_ Swap;
typedef TrajectoryExtractor_ TrajectoryExtractor;
typedef typename Simulator::Key SimulatorKey;
- typedef typename SortDS::iterator SortDSIterator;
private:
/* Implementation */
struct Node
{
- SortDSIterator element;
+ ElementIterator element;
SimulatorKey swap_event_key;
- Node(SortDSIterator e, SimulatorKey k):
+ Node(ElementIterator e, SimulatorKey k):
element(e), swap_event_key(k) {}
};
@@ -44,16 +48,15 @@
public:
typedef typename NodeList::iterator iterator;
- typedef boost::function<void(SortDS*, SortDSIterator pos)>
- SwapCallback;
/// \name Core Functionality
/// @{
- KineticSort(SortDS* sort, Simulator* simulator, SwapCallback swap_callback);
+ KineticSort();
+ KineticSort(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator);
+ void initialize(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator);
- template<class InputIterator>
- void insert(iterator pos, InputIterator f, InputIterator l, Simulator* simulator);
+ void insert(iterator pos, ElementIterator f, ElementIterator l, Simulator* simulator);
void erase(iterator pos, Simulator* simulator);
void update_trajectory(iterator pos, Simulator* simulator);
@@ -62,6 +65,9 @@
bool audit(Simulator* simulator) const;
/// @}
+ iterator begin() { return list_.begin(); }
+ iterator end() { return list_.end(); }
+
private:
class SwapEvent;
void schedule_swaps(iterator b, iterator e, Simulator* s);
@@ -69,8 +75,7 @@
private:
NodeList list_;
- SortDS* sort_;
- SwapCallback swap_callback_;
+ Swap swap_;
};
#include "kinetic-sort.hpp"
--- a/include/geometry/kinetic-sort.hpp Mon Feb 25 16:54:34 2008 -0500
+++ b/include/geometry/kinetic-sort.hpp Wed Feb 27 07:57:45 2008 -0500
@@ -1,49 +1,78 @@
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
-KineticSort(SortDS* sort, Simulator* simulator, SwapCallback swap_callback):
- sort_(sort), swap_callback_(swap_callback)
+#include "utilities/log.h"
+#include "utilities/counter.h"
+
+#ifdef LOGGING
+static rlog::RLogChannel* rlKineticSort = DEF_CHANNEL("geometry/kinetic-sort", rlog::Log_Debug);
+static rlog::RLogChannel* rlKineticSortAudit = DEF_CHANNEL("geometry/kinetic-sort/audit", rlog::Log_Debug);
+static rlog::RLogChannel* rlKineticSortSchedule = DEF_CHANNEL("geometry/kinetic-sort/schedule", rlog::Log_Debug);
+static rlog::RLogChannel* rlKineticSortProcess = DEF_CHANNEL("geometry/kinetic-sort/process", rlog::Log_Debug);
+#endif // LOGGING
+
+#ifdef COUNTERS
+static Counter* cKineticSort = GetCounter("kinetic-sort");
+static Counter* cKineticSortSwap = GetCounter("kinetic-sort/swap");
+#endif // COUNTERS
+
+
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
+KineticSort()
+{}
+
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
+KineticSort(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator)
{
- for (SortDSIterator cur = sort->begin(); cur != sort->end(); ++cur)
+ initialize(b, e, swap, simulator);
+}
+
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
+void
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
+initialize(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator)
+{
+ swap_ = swap;
+ for (ElementIterator cur = b; cur != e; ++cur)
list_.push_back(Node(cur, simulator->null_key()));
schedule_swaps(list_.begin(), list_.end(), simulator);
}
-
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
-template<class InputIterator>
+/// Adds elements in the range [f,l) of the underlying data structure to the management by the KineticSort.
+/// pos must be a valid iterator, i.e., it cannot be end().
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
void
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
-insert(iterator pos, InputIterator f, InputIterator l, Simulator* simulator)
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
+insert(iterator pos, ElementIterator f, ElementIterator l, Simulator* simulator)
{
iterator previous = pos; --previous;
if (previous != list_.end()) simulator->remove(previous->swap_event_key);
- sort_->insert(pos->element, f, l);
+ ElementIterator cur = boost::next(previous)->element;
+ while(cur != pos->element)
+ list_.insert(pos->element, Node(cur++, simulator->null_key()));
- SortDSIterator cur = boost::next(previous)->element;
- while(cur != pos->element)
- list_.insert(pos->element, Node(cur++));
if (previous != list_.end())
schedule_swaps(previous, pos, simulator);
else
schedule_swaps(list_.begin(), pos, simulator);
}
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+/// Removes pos from the KineticSort structure. Assumes that the element of the underlying data structure
+/// that *pos refers to has already been removed.
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
void
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
erase(iterator pos, Simulator* simulator)
{
simulator->remove(pos->swap_event_key);
- sort_->erase(pos->element);
iterator prev = pos; --prev;
list_.erase(pos);
schedule_swaps(prev, simulator);
}
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
void
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
update_trajectory(iterator pos, Simulator* simulator)
{
iterator prev = boost::prior(pos);
@@ -61,15 +90,16 @@
}
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
void
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
swap(iterator pos, Simulator* simulator)
{
- swap_callback_(sort_, pos->element);
-
- // TODO: add assertion that boost::next(pos) != list_.end()
-
+ // TODO: AssertMsg(boost::next(pos) != list_.end(), "Cannot swap the last element");
+
+ Count(cKineticSortSwap);
+ swap_(pos->element, simulator);
+
// Remove events
iterator prev = boost::prior(pos);
if (prev != list_.end())
@@ -88,60 +118,63 @@
//audit(simulator);
}
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
bool
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
audit(Simulator* simulator) const
{
- typedef typename Simulator::RationalFunction RationalFunction;
+ typedef typename Simulator::Function Function;
typedef typename Simulator::Time Time;
Time t = simulator->audit_time();
- std::cout << "Auditing at " << t << std::endl;
+ rLog(rlKineticSortAudit, "Auditing at %s", tostring(t).c_str());
TrajectoryExtractor te;
typename NodeList::const_iterator next = list_.begin();
typename NodeList::const_iterator cur = next++;
- RationalFunction cur_trajectory = te(cur->element);
+ Function cur_trajectory = te(cur->element);
while (next != list_.end())
{
- (*(cur->swap_event_key))->print(std::cout << " ") << std::endl;
+ rLog(rlKineticSortAudit, " %s", intostring(**(cur->swap_event_key)).c_str());
- RationalFunction next_trajectory = te(next->element);
- std::cout << " Auditing: " << cur_trajectory << ", " << next_trajectory << std::endl;
- std::cout << " Difference: " << next_trajectory - cur_trajectory << std::endl;
- std::cout << " Sign at: " << t << ", " << PolynomialKernel::sign_at(next_trajectory - cur_trajectory, t) << std::endl;
- if (PolynomialKernel::sign_at(next_trajectory - cur_trajectory, t) == -1)
+ Function next_trajectory = te(next->element);
+ rLog(rlKineticSortAudit, " Auditing: %s, %s", tostring(cur_trajectory).c_str(),
+ tostring(next_trajectory).c_str());
+ rLog(rlKineticSortAudit, " Difference: %s", tostring(next_trajectory - cur_trajectory).c_str());
+ rLog(rlKineticSortAudit, " Sign at: %s, %s", tostring(t).c_str(),
+ tostring(FunctionKernel::sign_at(next_trajectory - cur_trajectory, t)).c_str());
+ if (FunctionKernel::sign_at(next_trajectory - cur_trajectory, t) == -1)
{
- std::cout << "Audit failed at " << *cur->element << ", " << *next->element << std::endl;
+ rError("Audit failed at %s, %s", tostring(*cur->element).c_str(),
+ tostring(*next->element).c_str());
return false;
}
cur_trajectory = next_trajectory;
cur = next++;
}
- if (cur != list_.end()) (*(cur->swap_event_key))->print(std::cout << " ") << std::endl;
+ if (cur != list_.end()) rLog(rlKineticSortAudit, " %s", intostring(**(cur->swap_event_key)).c_str());
return true;
}
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
void
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
schedule_swaps(iterator b, iterator e, Simulator* simulator)
{
- typedef typename Simulator::RationalFunction RationalFunction;
+ typedef typename Simulator::Function Function;
TrajectoryExtractor te;
iterator next = b;
iterator cur = next++;
- RationalFunction cur_trajectory = te(cur->element);
+ Function cur_trajectory = te(cur->element);
while (next != e)
{
- RationalFunction next_trajectory = te(next->element);
- std::cout << "Next trajectory: " << next_trajectory << std::endl;
+ Function next_trajectory = te(next->element);
+ rLog(rlKineticSortSchedule, "Next trajectory: %s", tostring(next_trajectory).c_str());
// TODO: add assertion that (next_trajectory - cur_trajectory)(s->curren_time()) > 0
cur->swap_event_key = simulator->add(next_trajectory - cur_trajectory, SwapEvent(this, cur));
cur = next++;
@@ -150,12 +183,12 @@
if (cur != e) schedule_swaps(cur, simulator);
}
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
void
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
schedule_swaps(iterator i, Simulator* simulator)
{
- typedef typename Simulator::RationalFunction RationalFunction;
+ typedef typename Simulator::Function Function;
if (i == list_.end()) return;
if (boost::next(i) == list_.end())
@@ -167,8 +200,8 @@
TrajectoryExtractor te;
iterator next = boost::next(i);
- RationalFunction i_trajectory = te(i->element);
- RationalFunction next_trajectory = te(next->element);
+ Function i_trajectory = te(i->element);
+ Function next_trajectory = te(next->element);
//std::cout << "Updating swaps for: " << i_trajectory << ", " << next_trajectory << std::endl;
//std::cout << "Difference: " << next_trajectory - i_trajectory << std::endl;
@@ -178,8 +211,8 @@
}
/* SwapEvent */
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
-class KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::SwapEvent: public Simulator::Event
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
+class KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::SwapEvent: public Simulator::Event
{
public:
typedef typename Simulator::Event Parent;
@@ -192,29 +225,28 @@
virtual bool process(Simulator* s) const;
void set_position(iterator i) { pos_ = i; }
iterator position() const { return pos_; }
- std::ostream& print(std::ostream& out) const;
+ std::ostream& operator<<(std::ostream& out) const;
private:
KineticSort* sort_;
iterator pos_;
};
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
bool
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::SwapEvent::
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::SwapEvent::
process(Simulator* s) const
{
- std::cout << "Swapping. Current time: " << s->current_time() << std::endl;
+ rLog(rlKineticSortProcess, "Swapping. Current time: %s", tostring(s->current_time()).c_str());
sort_->swap(pos_, s);
return true;
}
-template<class SortDS_, class TrajectoryExtractor_, class Simulator_>
+template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
std::ostream&
-KineticSort<SortDS_, TrajectoryExtractor_, Simulator_>::SwapEvent::
-print(std::ostream& out) const
+KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::SwapEvent::
+operator<<(std::ostream& out) const
{
- Parent::print(out) << ", SwapEvent at " << TrajectoryExtractor_()(position()->element);
+ Parent::operator<<(out) << "SwapEvent at " << TrajectoryExtractor_()(position()->element);
return out;
}
-
--- a/include/geometry/polynomial.h Mon Feb 25 16:54:34 2008 -0500
+++ b/include/geometry/polynomial.h Wed Feb 27 07:57:45 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>
@@ -21,6 +22,7 @@
public:
typedef typename SynapsTraits<T>::Polynomial Polynomial;
typedef RationalFunction<Polynomial> RationalFunction;
+ typedef RationalFunction Function;
typedef typename SynapsTraits<T>::Solver Solver;
typedef typename SynapsTraits<T>::RootType RootType;
@@ -30,20 +32,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 int 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 +56,8 @@
{
typedef ZZ CoefficientType;
typedef SYNAPS::UPolDse<CoefficientType> Polynomial;
- typedef SYNAPS::Algebraic<CoefficientType> Solver;
+ typedef Polynomial SolverPolynomial;
+ typedef SYNAPS::Algebraic<ZZ> Solver;
typedef Solver::root_t RootType;
static RootType root(const CoefficientType& r) { CoefficientType p[2] = {-r, 1}; return SYNAPS::solve(Polynomial(2, p), Solver(), 0);}
@@ -59,6 +65,24 @@
static int sign_at(const Polynomial& p,
const RootType& r) { return SYNAPS::ALGEBRAIC::sign_at(p, r); }
//static RootType between(const RootType& r1, const RootType& r2) { RootType r = r1; r += r2; r /= root(2); return r; }
+ static SolverPolynomial convert(const Polynomial& p) { return p; }
+};
+
+template<>
+struct SynapsTraits<QQ>
+{
+ typedef QQ CoefficientType;
+ typedef SYNAPS::UPolDse<CoefficientType> Polynomial;
+ typedef SYNAPS::Algebraic<ZZ> Solver;
+ typedef SYNAPS::UPolDse<ZZ> SolverPolynomial;
+ typedef Solver::root_t RootType;
+
+ static RootType root(const CoefficientType& r);
+ static unsigned int multiplicity(const RootType& r) { return r.multiplicity(); }
+ static int sign_at(const Polynomial& p,
+ const RootType& r) { return SYNAPS::ALGEBRAIC::sign_at(convert(p), r); }
+ //static RootType between(const RootType& r1, const RootType& r2) { RootType r = r1; r += r2; r /= root(2); return r; }
+ static SolverPolynomial convert(const Polynomial& p);
};
#include "polynomial.hpp"
--- a/include/geometry/polynomial.hpp Mon Feb 25 16:54:34 2008 -0500
+++ b/include/geometry/polynomial.hpp Wed Feb 27 07:57:45 2008 -0500
@@ -1,3 +1,5 @@
+#include <utilities/log.h>
+
template<class T>
void
UPolynomial<T>::
@@ -5,8 +7,10 @@
{
typedef SYNAPS::Seq<RootType> RootSeq;
- RootSeq seq_num = SYNAPS::solve(rf.numerator(), Solver());
- RootSeq seq_den = SYNAPS::solve(rf.denominator(), Solver());
+ AssertMsg(stack.empty(), "Stack must be empty before solve");
+
+ RootSeq seq_num = SYNAPS::solve(SynapsTraits<T>::convert(rf.numerator()), Solver());
+ RootSeq seq_den = SYNAPS::solve(SynapsTraits<T>::convert(rf.denominator()), Solver());
// TODO: assert that all roots in seq_den have positive multiplicity
// TODO: deal with multiplicities for the numerator
@@ -32,3 +36,41 @@
{
return SynapsTraits<T>::sign_at(rf.numerator(), r) * SynapsTraits<T>::sign_at(rf.denominator(), r);
}
+
+template<class T>
+int
+UPolynomial<T>::
+sign_at_negative_infinity(const RationalFunction& rf)
+{
+ const Polynomial& num = rf.numerator();
+ const Polynomial& den = rf.denominator();
+ int ndegree = num.get_degree();
+ int ddegree = den.get_degree();
+ if (ndegree == -1) return 0; // ndegree == -1 => num == 0, and 0 is 0 at -infinity
+ return !((((ndegree + 1) % 2 == 0) ^ (num[ndegree] > 0)) ^
+ (((ddegree + 1) % 2 == 0) ^ (den[ddegree] > 0))) ? 1:-1;
+}
+
+SynapsTraits<QQ>::RootType
+SynapsTraits<QQ>::
+root(const CoefficientType& r)
+{
+ ZZ p[2] = { -SYNAPS::numerator(r), SYNAPS::denominator(r) };
+ return SYNAPS::solve(SolverPolynomial(2, p), Solver(), 0);
+}
+
+SynapsTraits<QQ>::SolverPolynomial
+SynapsTraits<QQ>::
+convert(const Polynomial& p)
+{
+ SolverPolynomial result(1, p.size() - 1);
+ assert(result.size() == p.size());
+ ZZ denominator_product = 1;
+ for (unsigned int i = 0; i < p.size(); ++i)
+ denominator_product *= SYNAPS::denominator(p[i]);
+ for (unsigned int i = 0; i < p.size(); ++i)
+ result[i] = SYNAPS::numerator(p[i]) * denominator_product /
+ SYNAPS::denominator(p[i]);
+ return result;
+}
+
--- a/include/geometry/simulator.h Mon Feb 25 16:54:34 2008 -0500
+++ b/include/geometry/simulator.h Wed Feb 27 07:57:45 2008 -0500
@@ -2,20 +2,26 @@
#define __SIMULATOR_H__
#include "utilities/eventqueue.h"
-#include "polynomial.h"
template<class Comparison> class IndirectComparison;
-template<class PolyKernel_, class Simulator_> class Event;
-template<class PolyKernel_, template<class Event> class EventComparison_ = std::less>
+/**
+ * Simulator class. Keeps a queue of events. Infinity is reached if the Event
+ * at the front of the queue has an empty root stack. Keeps track of current time,
+ * Event addition, and processes events one by one. Degeneracies are handled by
+ * assuming that the FunctionKernel::Function responsible for the event must be
+ * positive before the Event occurs.
+ *
+ * \ingroup kinetic
+ */
+template<class FuncKernel_, template<class Event> class EventComparison_ = std::less>
class Simulator
{
public:
- typedef PolyKernel_ PolynomialKernel;
- typedef typename PolynomialKernel::Polynomial Polynomial;
- typedef typename PolynomialKernel::RationalFunction RationalFunction;
- typedef typename PolynomialKernel::RootStack RootStack;
- typedef typename PolynomialKernel::RootType RootType;
+ typedef FuncKernel_ FunctionKernel;
+ typedef typename FunctionKernel::Function Function;
+ typedef typename FunctionKernel::RootStack RootStack;
+ typedef typename FunctionKernel::RootType RootType;
typedef RootType Time;
class Event;
@@ -26,7 +32,7 @@
typedef typename EventQueue::const_iterator const_Key;
- Simulator(Time start = PolynomialKernel::root(0)):
+ Simulator(Time start = FunctionKernel::root(0)):
current_(start),
reached_infinity_(false) {}
@@ -34,9 +40,9 @@
template<class Event_>
Key add(const Event_& e);
template<class Event_>
- Key add(const RationalFunction& f, const Event_& e);
+ Key add(const Function& f, const Event_& e);
void process();
- void update(Key k, const RationalFunction& f);
+ void update(Key k, const Function& f);
void remove(Key k) { queue_.remove(k); }
Key null_key() { return queue_.end(); }
@@ -45,7 +51,7 @@
Time audit_time() const;
bool reached_infinity() const { return reached_infinity_; }
- std::ostream& print(std::ostream& out) const;
+ std::ostream& operator<<(std::ostream& out) const;
private:
void update(Key i);
@@ -56,15 +62,23 @@
bool reached_infinity_;
};
-template<class PolyKernel_, template<class Event> class EventComparison_>
-class Simulator<PolyKernel_, EventComparison_>::Event
+
+/**
+ * Base class for events. Stores a root stack, subclasses need to define process().
+ * Event with an empty root stack compares greater than any other Event,
+ * pushing those events to the end of the queue.
+ */
+template<class FuncKernel_, template<class Event> class EventComparison_>
+class Simulator<FuncKernel_, EventComparison_>::Event
{
public:
- typedef PolyKernel_ PolynomialKernel;
- typedef typename PolynomialKernel::RootStack RootStack;
+ typedef FuncKernel_ FunctionKernel;
+ typedef typename FunctionKernel::RootStack RootStack;
+ /// process() is called when the event is at the top of the queue
+ /// in the simulator.
/// Returns true if the event needs to remain in the Simulator
- /// (top of the root_stack() will be used for new time)
+ /// (top of the root_stack() will be used for new time).
virtual bool process(Simulator* s) const =0;
RootStack& root_stack() { return root_stack_; }
@@ -80,15 +94,22 @@
return root_stack().top() < e.root_stack().top();
}
- int sign_before() const { return root_stack().top().sign_low(); }
- int sign_after() const { return root_stack().top().sign_high(); }
-
- virtual std::ostream& print(std::ostream& out) const { return out << "Event with " << root_stack_.size() << " roots"; }
+ virtual std::ostream& operator<<(std::ostream& out) const
+ {
+ out << "Event with " << root_stack_.size() << " roots";
+ if (!root_stack_.empty()) out << "; top root: " << root_stack_.top();
+ out << ", ";
+ return out;
+ }
private:
RootStack root_stack_;
};
+/**
+ * Compares elements pointed at by its arguments using the provided Comparison_
+ * (which must not take any arguments during construction).
+ */
template<class Comparison_>
class IndirectComparison: public std::binary_function<const typename Comparison_::first_argument_type*,
const typename Comparison_::second_argument_type*,
--- a/include/geometry/simulator.hpp Mon Feb 25 16:54:34 2008 -0500
+++ b/include/geometry/simulator.hpp Wed Feb 27 07:57:45 2008 -0500
@@ -1,82 +1,137 @@
-template<class PolyKernel_, template<class Event> class EventComparison_>
+#include "utilities/log.h"
+#include "utilities/counter.h"
+
+#ifdef LOGGING
+static rlog::RLogChannel* rlSimulator = DEF_CHANNEL("geometry/simulator", rlog::Log_Debug);
+
+#endif // LOGGING
+
+#ifdef COUNTERS
+static Counter* cSimulatorProcess = GetCounter("simulator/process");
+#endif // COUNTERS
+
+
+template<class FuncKernel_, template<class Event> class EventComparison_>
template<class Event_>
-typename Simulator<PolyKernel_, EventComparison_>::Key
-Simulator<PolyKernel_, EventComparison_>::
+typename Simulator<FuncKernel_, EventComparison_>::Key
+Simulator<FuncKernel_, EventComparison_>::
add(const Event_& e)
{
Event* ee = new Event_(e);
return queue_.push(ee);
}
-template<class PolyKernel_, template<class Event> class EventComparison_>
+template<class FuncKernel_, template<class Event> class EventComparison_>
template<class Event_>
-typename Simulator<PolyKernel_, EventComparison_>::Key
-Simulator<PolyKernel_, EventComparison_>::
-add(const RationalFunction& f, const Event_& e)
+typename Simulator<FuncKernel_, EventComparison_>::Key
+Simulator<FuncKernel_, EventComparison_>::
+add(const Function& f, const Event_& e)
{
Event* ee = new Event_(e);
- //std::cout << "Solving: " << f << std::endl;
- PolynomialKernel::solve(f, ee->root_stack());
+ rLog(rlSimulator, "Solving: %s", tostring(f).c_str());
+ int sign = FunctionKernel::sign_at_negative_infinity(f);
+ rLog(rlSimulator, "Sign at -infinity: %i", sign);
+ if (sign != 0)
+ {
+ FunctionKernel::solve(f, ee->root_stack());
+ rLog(rlSimulator, "Got solution with root stack size: %i", ee->root_stack().size());
+ }
+
while (!ee->root_stack().empty() && ee->root_stack().top() < current_time())
+ {
ee->root_stack().pop();
- //std::cout << "Pushing: " << ee->root_stack().top() << std::endl;
+ sign *= -1;
+ }
+ if (sign == -1)
+ {
+ AssertMsg(ee->root_stack().top() == current_time(),
+ "If sign is negative, we must be in the degenerate case");
+ ee->root_stack().pop();
+ }
+
+ if (ee->root_stack().empty())
+ rLog(rlSimulator, "Pushing event with empty root stack");
+ else
+ {
+ rLog(rlSimulator, "Root stack size: %i", ee->root_stack().size());
+ rLog(rlSimulator, "Pushing: %s", tostring(ee->root_stack().top()).c_str());
+ }
return queue_.push(ee);
}
-template<class PolyKernel_, template<class Event> class EventComparison_>
+template<class FuncKernel_, template<class Event> class EventComparison_>
void
-Simulator<PolyKernel_, EventComparison_>::
-update(Key k, const RationalFunction& f)
+Simulator<FuncKernel_, EventComparison_>::
+update(Key k, const Function& f)
{
Event* ee = *k;
ee->root_stack() = RootStack(); // no clear() in std::stack
- PolynomialKernel::solve(f, ee->root_stack());
+ FunctionKernel::solve(f, ee->root_stack());
while (!ee->root_stack().empty() && ee->root_stack().top() < current_time())
ee->root_stack().pop();
update(k);
}
-template<class PolyKernel_, template<class Event> class EventComparison_>
+template<class FuncKernel_, template<class Event> class EventComparison_>
void
-Simulator<PolyKernel_, EventComparison_>::
+Simulator<FuncKernel_, EventComparison_>::
process()
{
- std::cout << "Queue size: " << queue_.size() << std::endl;
+ Count(cSimulatorProcess);
+ rLog(rlSimulator, "Queue size: %i", queue_.size());
Key top = queue_.top();
Event* e = *top;
+ rLog(rlSimulator, "Processing event: %s", intostring(*e).c_str());
if (e->root_stack().empty()) { reached_infinity_ = true; return; }
else { current_ = e->root_stack().top(); e->root_stack().pop(); }
- if (e->process(this)) update(top);
- else { queue_.pop(); delete e; }
+ // Get the top element out of the queue, put it back depending on what process() says
+ EventQueue tmp; tmp.prepend(top, queue_);
+
+ if (e->process(this)) { queue_.prepend(top, tmp); update(top); }
+ else { delete e; }
}
-template<class PolyKernel_, template<class Event> class EventComparison_>
+template<class FuncKernel_, template<class Event> class EventComparison_>
void
-Simulator<PolyKernel_, EventComparison_>::
+Simulator<FuncKernel_, EventComparison_>::
update(Key i)
{
queue_.update(i);
}
-template<class PolyKernel_, template<class Event> class EventComparison_>
-typename Simulator<PolyKernel_, EventComparison_>::Time
-Simulator<PolyKernel_, EventComparison_>::
+template<class FuncKernel_, template<class Event> class EventComparison_>
+typename Simulator<FuncKernel_, EventComparison_>::Time
+Simulator<FuncKernel_, EventComparison_>::
audit_time() const
{
const_Key top = queue_.top();
Event* e = *top;
if (e->root_stack().empty()) return current_ + 1;
- else return PolynomialKernel::between(e->root_stack().top(), current_);
+ else return FunctionKernel::between(e->root_stack().top(), current_);
}
-template<class PolyKernel_, template<class Event> class EventComparison_>
+template<class FuncKernel_, template<class Event> class EventComparison_>
std::ostream&
-Simulator<PolyKernel_, EventComparison_>::
-print(std::ostream& out) const
+Simulator<FuncKernel_, EventComparison_>::
+operator<<(std::ostream& out) const
{
out << "Simulator: " << std::endl;
return queue_.print(out, " ");
}
+
+template<class FuncKernel_, template<class Event> class EventComparison_>
+std::ostream&
+operator<<(std::ostream& out, const Simulator<FuncKernel_, EventComparison_>& s)
+{
+ return s.operator<<(out);
+}
+
+template<class FuncKernel_, template<class Event> class EventComparison_>
+std::ostream&
+operator<<(std::ostream& out, const typename Simulator<FuncKernel_, EventComparison_>::Event& e)
+{
+ return e.operator<<(out);
+}
--- a/include/topology/cycle.h Mon Feb 25 16:54:34 2008 -0500
+++ b/include/topology/cycle.h Wed Feb 27 07:57:45 2008 -0500
@@ -16,6 +16,8 @@
* The actual order of the elements is defined by OrderCmp. Instances of those
* classes are not stored in Cycle for efficiency, and are passed as arguments to those methods
* that require them.
+ *
+ * \ingroup topology
*/
template <class Itm, class OrderCmp, class ConsistencyCmp = OrderCmp>
class Cycle: public List<Itm>
--- a/include/topology/filtration.h Mon Feb 25 16:54:34 2008 -0500
+++ b/include/topology/filtration.h Wed Feb 27 07:57:45 2008 -0500
@@ -22,6 +22,8 @@
* and provides pair_simplices() method that computes the RU-decomposition
* for the simplex order stored in the filtration. Iterators remain valid
* through all the operations.
+ *
+ * \ingroup topology
*/
template<class Smplx, class FltrSmplx = FiltrationSimplex<Smplx>, class Vnrd = Vineyard<FltrSmplx> >
class Filtration: public FltrSmplx::Container
@@ -37,6 +39,8 @@
typedef typename FiltrationSimplex::Container FiltrationContainer;
typedef typename FiltrationContainer::Index Index;
typedef typename FiltrationContainer::const_Index const_Index;
+ typedef Index iterator;
+ typedef const_Index const_iterator;
/// @}
/// \name Cycles and Trails
--- a/include/topology/filtrationcontainer.h Mon Feb 25 16:54:34 2008 -0500
+++ b/include/topology/filtrationcontainer.h Wed Feb 27 07:57:45 2008 -0500
@@ -13,6 +13,8 @@
* FiltrationContainer class. Serves as a parent of Filtration that
* describes the container functionality. Used by FiltrationSimplex
* to get Cycle representation.
+ *
+ * \ingroup topology
*/
template<class FltrSmplx>
class FiltrationContainer: public ConsistencyList<FltrSmplx>
--- a/include/topology/filtrationsimplex.h Mon Feb 25 16:54:34 2008 -0500
+++ b/include/topology/filtrationsimplex.h Wed Feb 27 07:57:45 2008 -0500
@@ -21,6 +21,8 @@
/**
* Evaluator is a base class for the structures that are able to return a value
* given a simplex.
+ *
+ * \ingroup topology
*/
template<class Smplx>
class Evaluator
@@ -37,6 +39,8 @@
/**
* FiltrationSimplex stores information needed for the RU-decomposition:
* cycle (column of R), trail (row of U), and pair.
+ *
+ * \ingroup topology
*/
template<class Smplx>
class FiltrationSimplex: public Smplx
--- a/include/topology/simplex.h Mon Feb 25 16:54:34 2008 -0500
+++ b/include/topology/simplex.h Wed Feb 27 07:57:45 2008 -0500
@@ -20,6 +20,8 @@
* SimplexWithVertices is a basic simplex class. It stores vertices of a given type,
* and knows how to compute its own boundary. It should probably be used as a base
* class for any explicit simplex representation.
+ *
+ * \ingroup topology
*/
template<class V>
class SimplexWithVertices
@@ -83,6 +85,8 @@
/**
* SimplexWithValue explicitly adds a RealType value to the SimplexWithVertices.
+ *
+ * \ingroup topology
*/
template<class Vert>
class SimplexWithValue: public SimplexWithVertices<Vert>
@@ -130,6 +134,8 @@
/**
* SimplexWithAttachment stores the vertex to which the simplex is attached (meant for lower-star filtrations)
+ *
+ * \ingroup topology
*/
template<typename V>
class SimplexWithAttachment: public SimplexWithVertices<V>
--- a/include/topology/vineyard.h Mon Feb 25 16:54:34 2008 -0500
+++ b/include/topology/vineyard.h Wed Feb 27 07:57:45 2008 -0500
@@ -23,6 +23,8 @@
/**
* Vineyard class. Keeps track of vines and knees. switched() is the key function called
* by filtration when pairing switches after a Filtration::transpose().
+ *
+ * \ingroup topology
*/
template<class FltrSmplx>
class Vineyard
@@ -65,6 +67,8 @@
/**
* Knee class stores the knee in R^3 as well as the cycle that is associated with the Simplex starting from the Knee.
+ *
+ * \ingroup topology
*/
template<class S>
class Knee
--- a/include/utilities/eventqueue.h Mon Feb 25 16:54:34 2008 -0500
+++ b/include/utilities/eventqueue.h Wed Feb 27 07:57:45 2008 -0500
@@ -34,6 +34,9 @@
void pop() { assert(!empty()); queue_.erase(queue_.begin()); }
void remove(iterator i) { queue_.erase(i); }
void update(iterator i);
+ void prepend(iterator i,
+ EventQueue& other) { queue_.splice(queue_.begin(), other.queue_, i); }
+ ///< intended for temporary storage of elements from other queues
iterator end() { return queue_.end(); }
const_iterator end() const { return queue_.end(); }
@@ -64,7 +67,7 @@
print(std::ostream& out, const std::string& prefix) const
{
for (typename QueueRepresentation::const_iterator cur = queue_.begin(); cur != queue_.end(); ++cur)
- (*cur)->print(out << prefix) << std::endl;
+ (*cur)->operator<<(out << prefix) << std::endl;
return out;
}
--- a/include/utilities/log.h Mon Feb 25 16:54:34 2008 -0500
+++ b/include/utilities/log.h Wed Feb 27 07:57:45 2008 -0500
@@ -12,6 +12,8 @@
template<class T>
std::string tostring(const T& t) { std::ostringstream out; out << t; return out.str(); }
+template<class T>
+std::string intostring(const T& t) { std::ostringstream out; t.operator<<(out); return out.str(); }
#define AssertMsg(cond, message, ...) do { if (!(cond)) { rError(message, ##__VA_ARGS__); rAssertSilent(cond); } } while (0)
--- a/tests/CMakeLists.txt Mon Feb 25 16:54:34 2008 -0500
+++ b/tests/CMakeLists.txt Wed Feb 27 07:57:45 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 Wed Feb 27 07:57:45 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} ${libraries})
+endforeach (t ${targets})
+
--- a/tests/geometry/Makefile Mon Feb 25 16:54:34 2008 -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 Mon Feb 25 16:54:34 2008 -0500
+++ b/tests/geometry/euclidean.cpp Wed Feb 27 07:57:45 2008 -0500
@@ -1,4 +1,4 @@
-#include "euclidean.h"
+#include "geometry/euclidean.h"
#include <vector>
#include <iostream>
#include <cmath>
--- a/tests/geometry/polynomial.cpp Mon Feb 25 16:54:34 2008 -0500
+++ b/tests/geometry/polynomial.cpp Wed Feb 27 07:57:45 2008 -0500
@@ -1,5 +1,5 @@
-#include <euclidean.h>
-#include <polynomial.h>
+#include <geometry/euclidean.h>
+#include <geometry/polynomial.h>
#include <vector>
#include <iostream>
@@ -75,8 +75,6 @@
while (!roots.empty()) { std::cout << roots.top() << std::endl; roots.pop(); }
}
- return 0;
-
// Edges
{
PointContainer vertices(2);
--- a/tests/geometry/test-eventqueue.cpp Mon Feb 25 16:54:34 2008 -0500
+++ b/tests/geometry/test-eventqueue.cpp Wed Feb 27 07:57:45 2008 -0500
@@ -1,4 +1,4 @@
-#include <eventqueue.h>
+#include <utilities/eventqueue.h>
#include <functional>
#include <iostream>
--- a/tests/geometry/test-kinetic-sort.cpp Mon Feb 25 16:54:34 2008 -0500
+++ b/tests/geometry/test-kinetic-sort.cpp Wed Feb 27 07:57:45 2008 -0500
@@ -1,25 +1,29 @@
-#include <polynomial.h>
-#include <simulator.h>
-#include <kinetic-sort.h>
+#include <geometry/polynomial.h>
+#include <geometry/simulator.h>
+#include <geometry/kinetic-sort.h>
#include <iostream>
#include <boost/utility.hpp>
+#include <boost/bind.hpp>
-//typedef double FieldType;
+#include <utilities/log.h>
+
+typedef double FieldType;
//typedef ZZ FieldType;
-typedef QQ FieldType;
+//typedef QQ FieldType;
typedef UPolynomial<FieldType> PolyKernel;
typedef PolyKernel::Polynomial Polynomial;
typedef std::list<Polynomial> SortDS;
+typedef SortDS::iterator SortDSIterator;
typedef Simulator<PolyKernel> SimulatorFT;
class TrajectoryExtractor
{
public:
- Polynomial operator()(SortDS::iterator i) const { return *i; }
+ Polynomial operator()(SortDSIterator i) const { return *i; }
};
-typedef KineticSort<SortDS, TrajectoryExtractor, SimulatorFT> KineticSortDS;
+typedef KineticSort<SortDSIterator, TrajectoryExtractor, SimulatorFT> KineticSortDS;
struct EvaluatedComparison: public std::binary_function<const Polynomial&, const Polynomial&, bool>
{
@@ -29,21 +33,33 @@
FieldType vv;
};
-void swap(SortDS* s, SortDS::iterator i)
+void swap(SortDS* s, SortDSIterator i, SimulatorFT* simulator)
{
std::cout << "Swapping " << *i << " " << *boost::next(i) << std::endl;
s->splice(i, *s, boost::next(i));
}
-int main()
+int main(int argc, char** argv)
{
+#ifdef LOGGING
+ rlog::RLogInit(argc, argv);
+
+ stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+ stdoutLog.subscribeTo( RLOG_CHANNEL("geometry/simulator") );
+ stdoutLog.subscribeTo( RLOG_CHANNEL("geometry/kinetic-sort") );
+#endif
+
SimulatorFT simulator;
SortDS list;
// Insert polynomials and sort the list for current time
list.push_back(Polynomial("x^3 - 3"));
list.push_back(Polynomial("x^2 - 2*x - 2"));
+ list.push_back(Polynomial("x^2 - 2*x - 2"));
list.push_back(Polynomial("2*x - 4"));
+ list.push_back(Polynomial("x"));
+ list.push_back(Polynomial("-x + 4"));
+ list.push_back(Polynomial("2"));
//list.sort(EvaluatedComparison(simulator.current_time()));
list.sort(EvaluatedComparison(0));
@@ -52,15 +68,23 @@
std::cout << *cur << std::endl;
// Setup kinetic sort
- KineticSortDS ks(&list, &simulator, swap);
+ KineticSortDS ks(list.begin(), list.end(), boost::bind(swap, &list, _1, _2), &simulator);
- while(!simulator.reached_infinity() && simulator.current_time() < 1)
+ std::cout << "Examining " << simulator;
+
+ while(!simulator.reached_infinity())
{
- //std::cout << "Current time before: " << simulator.current_time() << std::endl;
- if (!ks.audit(&simulator)) return 1;
- //simulator.print(std::cout << "Auditing ");
+ std::cout << "Current time before: " << simulator.current_time() << std::endl;
+ //if (!ks.audit(&simulator)) return 1;
+ ks.audit(&simulator);
+ std::cout << "Examining " << simulator;
simulator.process();
- //std::cout << "Current time after: " << simulator.current_time() << std::endl;
+ std::cout << "Current time after: " << simulator.current_time() << std::endl;
}
+ ks.audit(&simulator);
+ std::cout << "Examining " << simulator;
+ std::cout << "Done at " << simulator.current_time() << std::endl;
+ for (SortDS::const_iterator cur = list.begin(); cur != list.end(); ++cur)
+ std::cout << " " << *cur << std::endl;
}
--- a/tests/geometry/test-linalg.cpp Mon Feb 25 16:54:34 2008 -0500
+++ b/tests/geometry/test-linalg.cpp Wed Feb 27 07:57:45 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;