--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/89ae955518665a61 Fri Apr 04 21:14:44 2008 -0400
@@ -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 Fri Apr 04 21:14:44 2008 -0400
@@ -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 Fri Apr 04 21:14:44 2008 -0400
@@ -0,0 +1,71 @@
+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.
+
+From artemis Wed Mar 19 16:47:55 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Wed, 19 Mar 2008 12:47:24 -0400
+Subject: Fixed in 0a18d6902a55
+Message-Id: <c420501cc5285bbc-db67b15857938f7b-artemis@metatron>
+References: <c420501cc5285bbc-0-artemis@metatron>
+In-Reply-To: <c420501cc5285bbc-0-artemis@metatron>
+
+Fixed by dealing with infinite simplices in 0a18d6902a55
+
+From artemis Wed Mar 19 16:49:36 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Wed, 19 Mar 2008 12:49:01 -0400
+Subject: Reopened for ar-vineyard
+Message-Id: <c420501cc5285bbc-8f74ccde6b3f0bea-artemis@metatron>
+References: <c420501cc5285bbc-0-artemis@metatron>
+In-Reply-To: <c420501cc5285bbc-0-artemis@metatron>
+
+Fix for ar-vineyard like for alphashapes
+
+From artemis Wed Mar 19 16:49:47 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Wed, 19 Mar 2008 12:49:47 -0400
+Subject: properties changes (state)
+Message-Id: <c420501cc5285bbc-ef5f506b45aba335-artemis@metatron>
+References: <c420501cc5285bbc-0-artemis@metatron>
+In-Reply-To: <c420501cc5285bbc-0-artemis@metatron>
+
+state=open
+
+From artemis Wed Mar 19 17:17:47 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Wed, 19 Mar 2008 13:17:30 -0400
+Subject: Fixed for ar-vineyard
+Message-Id: <c420501cc5285bbc-feb676745fc16b8c-artemis@metatron>
+References: <c420501cc5285bbc-8f74ccde6b3f0bea-artemis@metatron>
+In-Reply-To: <c420501cc5285bbc-8f74ccde6b3f0bea-artemis@metatron>
+
+Fixed for ar-vineyard.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/c664b2f69b5f6ea3 Fri Apr 04 21:14:44 2008 -0400
@@ -0,0 +1,29 @@
+From artemis Thu Feb 28 10:34:23 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Thu, 28 Feb 2008 05:33:34 -0500
+State: fixed
+Subject: Two simulators
+Message-Id: <c664b2f69b5f6ea3-0-artemis@metatron>
+
+Consider using two simulators instead of one to make sure that trajectory
+changes are processed before simplex swaps.
+
+From artemis Sat Mar 1 09:45:05 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Sat, 01 Mar 2008 04:45:05 -0500
+Subject: properties changes (state)
+Message-Id: <c664b2f69b5f6ea3-f8bf16bdff01098d-artemis@metatron>
+References: <c664b2f69b5f6ea3-0-artemis@metatron>
+In-Reply-To: <c664b2f69b5f6ea3-0-artemis@metatron>
+
+state=fixed
+
+From artemis Sat Mar 1 09:46:20 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Sat, 01 Mar 2008 04:45:29 -0500
+Subject: Fixed in f236c7d659d0
+Message-Id: <c664b2f69b5f6ea3-a1742f0eb7b5e1c0-artemis@metatron>
+References: <c664b2f69b5f6ea3-0-artemis@metatron>
+In-Reply-To: <c664b2f69b5f6ea3-0-artemis@metatron>
+
+Fixed in f236c7d659d0.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/cd79223a108d3900 Fri Apr 04 21:14:44 2008 -0400
@@ -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/d2ab07329c3588ca Fri Apr 04 21:14:44 2008 -0400
@@ -0,0 +1,29 @@
+From artemis Wed Feb 27 21:31:41 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Wed, 27 Feb 2008 16:30:57 -0500
+State: fixed
+Subject: Segfault with topology/vineyard log enabled
+Message-Id: <d2ab07329c3588ca-0-artemis@metatron>
+
+The code segfaults in Vineyard::record_knee() if we subscribe to
+topology/vineyard log.
+
+From artemis Sat Mar 1 09:46:46 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Sat, 01 Mar 2008 04:46:46 -0500
+Subject: properties changes (state)
+Message-Id: <d2ab07329c3588ca-8c184b421f538482-artemis@metatron>
+References: <d2ab07329c3588ca-0-artemis@metatron>
+In-Reply-To: <d2ab07329c3588ca-0-artemis@metatron>
+
+state=fixed
+
+From artemis Sat Mar 1 09:47:07 2008
+From: Dmitriy Morozov <morozov@cs.duke.edu>
+Date: Sat, 01 Mar 2008 04:46:51 -0500
+Subject: Fixed in abba2950aced
+Message-Id: <d2ab07329c3588ca-fa65c35dfff6366f-artemis@metatron>
+References: <d2ab07329c3588ca-0-artemis@metatron>
+In-Reply-To: <d2ab07329c3588ca-0-artemis@metatron>
+
+Fixed in abba2950aced
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.issues/ebda8db3f9908e33 Fri Apr 04 21:14:44 2008 -0400
@@ -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 Fri Apr 04 16:17:28 2008 -0400
+++ b/CMakeLists.txt Fri Apr 04 21:14:44 2008 -0400
@@ -5,6 +5,7 @@
option (debug "Build Dionysus with debugging on" OFF)
option (optimize "Build Dionysus with optimization" ON)
option (use_dsrpdb "Build examples that use DSR-PDB" OFF)
+option (use_synaps "Build examples that use SYNAPS" OFF)
# Find everything that's always required
find_package (Boost REQUIRED)
@@ -43,6 +44,15 @@
${gmpxx_LIBRARY}
${m_LIBRARY})
+# SYNAPS
+if (use_synaps)
+ add_definitions (-DBOOST_UBLAS_TYPE_CHECK=0)
+ find_library (synaps_LIBRARY NAMES synaps)
+ set (synaps_libraries ${synaps_LIBRARY}
+ ${gmp_LIBRARY}
+ ${gmpxx_LIBRARY})
+endif (use_synaps)
+
# Debugging
if (debug)
if (optimize)
--- a/examples/ar-vineyard/CMakeLists.txt Fri Apr 04 16:17:28 2008 -0400
+++ b/examples/ar-vineyard/CMakeLists.txt Fri Apr 04 21:14:44 2008 -0400
@@ -1,8 +1,15 @@
set (targets
ar-vineyard)
-
+
+find_library (boost_signals_LIBRARY NAMES boost_signals
+ PATHS ${Boost_LIBRARY_DIR})
+find_library (boost_program_options_LIBRARY NAME boost_program_options
+ 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}
+ ${boost_signals_LIBRARY}
+ ${boost_program_options_LIBRARY})
endforeach (t)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/ar-vineyard/ar-function-kernel.h Fri Apr 04 21:14:44 2008 -0400
@@ -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 Fri Apr 04 21:14:44 2008 -0400
@@ -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 Fri Apr 04 16:17:28 2008 -0400
+++ b/examples/ar-vineyard/ar-simplex3d.h Fri Apr 04 21:14:44 2008 -0400
@@ -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;
@@ -56,10 +56,11 @@
ARSimplex3D(const Edge& e);
ARSimplex3D(const Edge& e, const Point& z, SimplexPhiMap& simplices,
- Facet_circulator facet_bg);
+ const Delaunay& Dt, Facet_circulator facet_bg);
ARSimplex3D(const Facet& f);
- ARSimplex3D(const Facet& f, const Point& z, const SimplexPhiMap& simplices);
+ ARSimplex3D(const Facet& f, const Point& z, const SimplexPhiMap& simplices,
+ const Delaunay& Dt);
ARSimplex3D(const Cell& c, const Point& z);
@@ -86,9 +87,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 Fri Apr 04 16:17:28 2008 -0400
+++ b/examples/ar-vineyard/ar-simplex3d.hpp Fri Apr 04 21:14:44 2008 -0400
@@ -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)
{
@@ -35,19 +41,15 @@
}
ARSimplex3D::
-ARSimplex3D(const Edge& e, const Point& z, SimplexPhiMap& simplices, Facet_circulator facet_bg)
+ARSimplex3D(const Edge& e, const Point& z, SimplexPhiMap& simplices, const Delaunay& Dt, Facet_circulator facet_bg)
{
Cell_handle c = e.first;
Parent::add(c->vertex(e.second));
Parent::add(c->vertex(e.third));
Facet_circulator cur = facet_bg;
+ while (Dt.is_infinite(*cur)) ++cur;
SimplexPhiMap::const_iterator cur_iter = simplices.find(ARSimplex3D(*cur));
- while (cur_iter == simplices.end())
- {
- ++cur;
- cur_iter = simplices.find(ARSimplex3D(*cur));
- }
RealValue min = cur_iter->first.alpha();
RealValue phi_const_min = cur_iter->first.phi_const();
@@ -63,10 +65,8 @@
int i0 = (*cur).first->index(*v1);
int i1 = (*cur).first->index(*v2);
int i = 6 - i0 - i1 - (*cur).second;
- Point p3 = (*cur).first->vertex(i)->point();
- cur_iter = simplices.find(ARSimplex3D(*cur));
- if (cur_iter == simplices.end()) // cur is infinite
+ if (Dt.is_infinite(cur->first->vertex(i)))
{
++cur; continue;
// FIXME: what do we do with infinite cofaces (i.e., what
@@ -74,9 +74,11 @@
// infinite?
}
+ Point p3 = (*cur).first->vertex(i)->point();
if (CGAL::side_of_bounded_sphere(p1, p2, p3) == CGAL::ON_BOUNDED_SIDE)
attached_ = true;
+ SimplexPhiMap::const_iterator cur_iter = simplices.find(ARSimplex3D(*cur));
RealValue val = cur_iter->first.alpha();
if (val < min) min = val;
RealValue phi_const_val = cur_iter->first.phi_const();
@@ -99,7 +101,7 @@
s_ = CGAL::squared_distance(z, K::Segment_3(p1,p2).supporting_line());
- Point origin;
+ Point origin(0,0,0);
Point cc = origin + ((p1 - origin) + (p2 - origin))/2; // CGAL is funny
v_ = CGAL::squared_distance(z, cc) - s_;
}
@@ -114,7 +116,7 @@
}
ARSimplex3D::
-ARSimplex3D(const Facet& f, const Point& z, const SimplexPhiMap& simplices)
+ARSimplex3D(const Facet& f, const Point& z, const SimplexPhiMap& simplices, const Delaunay& Dt)
{
Cell_handle c = f.first;
for (int i = 0; i < 4; ++i)
@@ -133,29 +135,33 @@
rho_ = squared_radius(p1, p2, p3);
attached_ = false;
- if (CGAL::side_of_bounded_sphere(p1, p2, p3,
+ if (!Dt.is_infinite(c->vertex(f.second)) &&
+ CGAL::side_of_bounded_sphere(p1, p2, p3,
c->vertex(f.second)->point()) == CGAL::ON_BOUNDED_SIDE)
attached_ = true;
- else if (CGAL::side_of_bounded_sphere(p1, p2, p3,
+ else if (!Dt.is_infinite(o->vertex(oi)) &&
+ CGAL::side_of_bounded_sphere(p1, p2, p3,
o->vertex(oi)->point()) == CGAL::ON_BOUNDED_SIDE)
attached_ = true;
else
alpha_ = rho_;
- SimplexPhiMap::const_iterator c_iter = simplices.find(ARSimplex3D(*c,z));
- SimplexPhiMap::const_iterator o_iter = simplices.find(ARSimplex3D(*o,z));
- if (c_iter == simplices.end()) // c is infinite
+ if (Dt.is_infinite(c))
{
+ SimplexPhiMap::const_iterator o_iter = simplices.find(ARSimplex3D(*o,z));
if (attached_) alpha_ = o_iter->first.alpha();
phi_const_ = o_iter->first.phi_const(); // FIXME: it's probably the other way around
}
- else if (o_iter == simplices.end()) // o is infinite
+ else if (Dt.is_infinite(o))
{
+ SimplexPhiMap::const_iterator c_iter = simplices.find(ARSimplex3D(*c,z));
if (attached_) alpha_ = c_iter->first.alpha();
phi_const_ = c_iter->first.phi_const(); // FIXME: it's probably the other way around
}
else
{
+ SimplexPhiMap::const_iterator o_iter = simplices.find(ARSimplex3D(*o,z));
+ SimplexPhiMap::const_iterator c_iter = simplices.find(ARSimplex3D(*c,z));
if (attached_) alpha_ = std::min(c_iter->first.alpha(), o_iter->first.alpha());
phi_const_ = std::min(c_iter->first.phi_const(), o_iter->first.phi_const());
}
@@ -208,8 +214,9 @@
ARSimplex3D::
operator<<(std::ostream& out) const
{
+ out << this << ": ";
for (VertexSet::const_iterator cur = Parent::vertices().begin(); cur != Parent::vertices().end(); ++cur)
- out << **cur << ", ";
+ out << &(**cur) << ", ";
out << "value = " << value();
return out;
@@ -226,19 +233,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;
+ update_simplex_phi_map(ARSimplex3D(*cur, z, simplices, Dt), simplices);
+ 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;
+ update_simplex_phi_map(ARSimplex3D(*cur, z, simplices, Dt, Dt.incident_facets(*cur)), simplices);
+ 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 Fri Apr 04 16:17:28 2008 -0400
+++ b/examples/ar-vineyard/ar-vineyard.cpp Fri Apr 04 21:14:44 2008 -0400
@@ -4,37 +4,77 @@
#include <iostream>
#include <fstream>
+#include <string>
+#include <vector>
+
+#include <boost/program_options.hpp>
+namespace po = boost::program_options;
int main(int argc, char** argv)
{
#ifdef LOGGING
rlog::RLogInit(argc, argv);
+ stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+#endif
- stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
- //stdoutLog.subscribeTo( RLOG_CHANNEL("topology/filtration") );
- //stdoutLog.subscribeTo( RLOG_CHANNEL("topology/cycle") );
- //stdoutLog.subscribeTo( RLOG_CHANNEL("topology/vineyard") );
- //stdoutLog.subscribeTo( RLOG_CHANNEL("topology/filtration/transpositions") );
- //stdoutLog.subscribeTo( RLOG_CHANNEL("topology/lowerstar") );
+ std::string infilename;
+ double zx,zy,zz,r;
+ std::string outfilename;
+
+ po::options_description hidden("Hidden options");
+ hidden.add_options()
+ ("input-file", po::value<std::string>(&infilename), "Points file")
+ ("x", po::value<double>(&zx), "x")
+ ("y", po::value<double>(&zy), "y")
+ ("z", po::value<double>(&zz), "z")
+ ("r", po::value<double>(&r), "r")
+ ("output-file", po::value<std::string>(&outfilename), "Vineyard edges output");
+
+ std::vector<std::string> log_channels;
+ po::options_description visible("Allowed options");
+ visible.add_options()
+ ("help,h", "produce help message");
+#if LOGGING
+ visible.add_options()
+ ("log,l", po::value< std::vector<std::string> >(&log_channels),
+ "log channels to turn on");
+#endif
+
+ po::positional_options_description p;
+ p.add("input-file", 1).add("x", 1).add("y", 1).add("z", 1).add("r", 1).add("output-file", 1);
+
+ po::options_description all; all.add(visible).add(hidden);
+
+ po::variables_map vm;
+ po::store(po::command_line_parser(argc, argv).
+ options(all).positional(p).run(), vm);
+ po::notify(vm);
+
+#if LOGGING
+ for (std::vector<std::string>::const_iterator cur = log_channels.begin(); cur != log_channels.end(); ++cur)
+ stdoutLog.subscribeTo( RLOG_CHANNEL(cur->c_str()) );
+ /* Interesting channels
+ * "ar/vineyard", "ar/function-kernel/value", "geometry/simulator",
+ * "topology/filtration", "topology/cycle", "topology/vineyard",
+ * "topology/filtration/transpositions", "topology/lowerstar"
+ */
#endif
// Read command-line arguments
- if (argc < 6)
+ if (vm.count("help") || !vm.count("input-file") || !vm.count("output-file")
+ || !vm.count("x") || !vm.count("y") || !vm.count("z") || !vm.count("r"))
{
- std::cout << "Usage: ar-vineyard POINTS X Y Z OUTFILENAME" << std::endl;
+ std::cout << "Usage: ar-vineyard [OPTIONS] POINTS X Y Z R OUTFILENAME" << std::endl;
std::cout << " POINTS - filename containing points" << std::endl;
std::cout << " X,Y,Z - center-point z at which to compute the vineyard" << std::endl;
+ std::cout << " R - maximum radius" << std::endl;
std::cout << " OUTFILENAME - filename for the resulting vineyard" << std::endl;
+ std::cout << visible << std::endl;
std::cout << std::endl;
std::cout << "Computes an (alpha,r)-vineyard of the given pointset around the given point." << std::endl;
- exit(0);
+ return 1;
}
- std::string infilename = argv[1];
- double zx,zy,zz; std::istringstream(argv[2]) >> zx;
- std::istringstream(argv[3]) >> zy; std::istringstream(argv[4]) >> zz;
- std::string outfilename = argv[5];
-
// Read in the point set and compute its Delaunay triangulation
std::ifstream in(infilename.c_str());
@@ -52,7 +92,7 @@
arv.compute_pairing();
// Compute vineyard
- arv.compute_vineyard(true);
+ arv.compute_vineyard(r);
std::cout << "Vineyard computed" << std::endl;
arv.vineyard()->save_edges(outfilename);
}
--- a/examples/ar-vineyard/ar-vineyard.h Fri Apr 04 16:17:28 2008 -0400
+++ b/examples/ar-vineyard/ar-vineyard.h Fri Apr 04 21:14:44 2008 -0400
@@ -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 Parent& p): Parent(p) {} // crucial for boundary() to work correctly
+ 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; }
+
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> SimulatorAR;
+ /// @}
+
+ /// \name Filtration types
+ /// @{
+ typedef ARConeSimplex3D<SimulatorAR> ConeSimplex3D;
+ typedef Filtration<ConeSimplex3D> FiltrationAR;
+ typedef FiltrationAR::Simplex Simplex;
+ typedef FiltrationAR::Index Index;
+ typedef FiltrationAR::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, SimulatorAR> 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(double max_radius);
- const ARFiltration* filtration() const { return filtration_; }
+ const FiltrationAR* filtration() const { return filtration_; }
const Vineyard* vineyard() const { return vineyard_; }
public:
- // For Kinetic Sort
- void swap(Key a, Key b);
+ void swap(Index i, SimulatorAR* simulator); ///< For kinetic sort
private:
void add_simplices();
void change_evaluator(Evaluator* eval);
private:
- ARFiltration* filtration_;
+ FiltrationAR* filtration_;
Vineyard* vineyard_;
Evaluator* evaluator_;
- KeyIndexMap kinetic_map_;
-
Point z_;
Delaunay dt_;
@@ -120,68 +156,58 @@
//BOOST_CLASS_EXPORT(ARVineyard)
+#ifdef COUNTERS
+static Counter* cARVineyardTrajectoryKnee = GetCounter("ar/vineyard/trajectoryknee");
+#endif
-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,
+ Vineyard* vineyard, SimulatorAR* sort_simulator):
+ iter_(iter), sort_(sort), vineyard_(vineyard), sort_simulator_(sort_simulator) { iter_->element->new_max_signal().connect(*this); }
+ void operator()(SimulatorAR* simulator)
+ {
+ Count(cARVineyardTrajectoryKnee);
+ sort_->update_trajectory(iter_, sort_simulator_);
+ if (iter_->element->sign())
+ vineyard_->record_knee(iter_->element);
+ else
+ vineyard_->record_knee(iter_->element->pair());
+ }
+
+ private:
+ SimplexSortIterator iter_;
+ SimplexSort* sort_; // could make both of these static
+ Vineyard* vineyard_; // currently inefficient since there is
+ // only one SimplexSort and one Vineyard,
+ // but each is stored in every slot
+ SimulatorAR* sort_simulator_;
};
-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(SimulatorAR* simplex_sort_simulator,
+ SimulatorAR* trajectory_sort_simulator):
+ simplex_sort_simulator_(simplex_sort_simulator),
+ trajectory_sort_simulator_(trajectory_sort_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 std::max(simplex_sort_simulator_->current_time(), trajectory_sort_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_;
+ SimulatorAR* simplex_sort_simulator_;
+ SimulatorAR* trajectory_sort_simulator_;
};
--- a/examples/ar-vineyard/ar-vineyard.hpp Fri Apr 04 16:17:28 2008 -0400
+++ b/examples/ar-vineyard/ar-vineyard.hpp Fri Apr 04 21:14:44 2008 -0400
@@ -1,26 +1,68 @@
#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* rlARVineyardThresholdSwap = DEF_CHANNEL("ar/vineyard/threshold/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(rlARVineyardThresholdSwap, "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 FiltrationAR(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(ConeSimplex3D(*cur)); // Delaunay simplex
+ filtration_->append(ConeSimplex3D(*cur, true)); // Coned off delaunay simplex
}
}
@@ -39,98 +81,64 @@
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(double max_radius)
{
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;
+ SimulatorAR simplex_sort_simulator, trajectory_sort_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(&trajectory_sort_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), &simplex_sort_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, vineyard_, &simplex_sort_simulator));
+ rLog(rlARVineyardComputing, "Signals and slots connected");
+ rLog(rlARVineyardComputing, "SimplexSort size: %i", simplex_sort_simulator.size());
+ rLog(rlARVineyardComputing, "TrajectorySort size: %i", trajectory_sort_simulator.size());
- // 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(&simplex_sort_simulator, &trajectory_sort_simulator));
+ while ((simplex_sort_simulator.current_time() < max_radius || trajectory_sort_simulator.current_time() < max_radius) && !(simplex_sort_simulator.reached_infinity() && trajectory_sort_simulator.reached_infinity()))
+ {
+ if (*(simplex_sort_simulator.top()) < *(trajectory_sort_simulator.top()))
+ {
+ rLog(rlARVineyardComputing, "Current time before: %lf (simplex sort)", simplex_sort_simulator.current_time());
+ rLog(rlARVineyardComputing, "Top event: %s (simplex sort)", intostring(*(simplex_sort_simulator.top())).c_str());
+ simplex_sort_simulator.process();
+ } else
+ {
+ rLog(rlARVineyardComputing, "Current time before: %lf (trajectory sort)", trajectory_sort_simulator.current_time());
+ rLog(rlARVineyardComputing, "Top event: %s (trajectory sort)", intostring(*(trajectory_sort_simulator.top())).c_str());
+ trajectory_sort_simulator.process();
+ }
+ }
- //change_evaluator(new StaticEvaluator(1));
vineyard_->record_diagram(filtration_->begin(), filtration_->end());
}
void
ARVineyard::
-swap(Key a, Key b)
+swap(Index i, SimulatorAR* 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:", &(*i), &(*boost::next(i)));
+ rLog(rlARVineyardSwap, " %s and", tostring(*i).c_str());
+ rLog(rlARVineyardSwap, " %s", tostring(*boost::next(i)).c_str());
+ filtration_->transpose(i);
}
void
@@ -143,30 +151,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 Fri Apr 04 16:17:28 2008 -0400
+++ b/include/dionysus.h Fri Apr 04 21:14:44 2008 -0400
@@ -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 Fri Apr 04 16:17:28 2008 -0400
+++ b/include/geometry/euclidean.h Fri Apr 04 21:14:44 2008 -0400
@@ -14,15 +14,19 @@
#include "number-traits.h"
+/**
+ * Geometric Kernel. Defines operations on geometric primitives.
+ * \ingroup geometry
+ */
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;
+ typedef LinearAlgebra<NumberType> LinearAlgebraK;
+ typedef typename LinearAlgebraK::MatrixType MatrixType;
+ typedef typename LinearAlgebraK::VectorType VectorType;
class Point;
class Sphere;
@@ -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 Fri Apr 04 16:17:28 2008 -0400
+++ b/include/geometry/kinetic-sort.h Fri Apr 04 21:14:44 2008 -0400
@@ -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 Fri Apr 04 16:17:28 2008 -0400
+++ b/include/geometry/kinetic-sort.hpp Fri Apr 04 21:14:44 2008 -0400
@@ -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 Fri Apr 04 16:17:28 2008 -0400
+++ b/include/geometry/polynomial.h Fri Apr 04 21:14:44 2008 -0400
@@ -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>
@@ -20,30 +21,34 @@
{
public:
typedef typename SynapsTraits<T>::Polynomial Polynomial;
- typedef RationalFunction<Polynomial> RationalFunction;
+ typedef RationalFunction<Polynomial> RationalFunctionP;
+ typedef RationalFunctionP Function;
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 void solve(const RationalFunctionP& 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 int sign_at(const RationalFunctionP& 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 RationalFunctionP& 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 Fri Apr 04 16:17:28 2008 -0400
+++ b/include/geometry/polynomial.hpp Fri Apr 04 21:14:44 2008 -0400
@@ -1,12 +1,16 @@
+#include <utilities/log.h>
+
template<class T>
void
UPolynomial<T>::
-solve(const RationalFunction& rf, RootStack& stack)
+solve(const RationalFunctionP& rf, RootStack& stack)
{
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
@@ -28,7 +32,45 @@
template<class T>
int
UPolynomial<T>::
-sign_at(const RationalFunction& rf, const RootType& r)
+sign_at(const RationalFunctionP& rf, const RootType& r)
{
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 RationalFunctionP& 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 Fri Apr 04 16:17:28 2008 -0400
+++ b/include/geometry/simulator.h Fri Apr 04 21:14:44 2008 -0400
@@ -2,69 +2,84 @@
#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;
typedef EventComparison_<Event> EventComparison;
typedef EventQueue<Event*, IndirectComparison<EventComparison> >
- EventQueue;
- typedef typename EventQueue::iterator Key;
- typedef typename EventQueue::const_iterator const_Key;
+ EventQueueS;
+ typedef typename EventQueueS::iterator Key;
+ typedef typename EventQueueS::const_iterator const_Key;
- Simulator(Time start = PolynomialKernel::root(0)):
- current_(start),
- reached_infinity_(false) {}
+ Simulator(Time start = FunctionKernel::root(0)):
+ current_(start) {}
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(); }
Time current_time() const { return current_; }
Time audit_time() const;
- bool reached_infinity() const { return reached_infinity_; }
+ bool reached_infinity() const { return queue_.empty() || (*queue_.top())->root_stack().empty(); }
+
+ Event* top() const { return *(queue_.top()); }
+ unsigned size() const { return queue_.size(); }
- std::ostream& print(std::ostream& out) const;
+ std::ostream& operator<<(std::ostream& out) const;
private:
void update(Key i);
private:
Time current_;
- EventQueue queue_;
- bool reached_infinity_;
+ EventQueueS queue_;
};
-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 +95,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 Fri Apr 04 16:17:28 2008 -0400
+++ b/include/geometry/simulator.hpp Fri Apr 04 21:14:44 2008 -0400
@@ -1,82 +1,138 @@
-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); // going to be sign after current time
+ 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");
+ rLog(rlSimulator, "Popping the root because of negative sign (degeneracy)");
+ 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);
+ if (reached_infinity()) return;
+ 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(); }
+ 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
+ EventQueueS 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/conesimplex.h Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/conesimplex.h Fri Apr 04 21:14:44 2008 -0400
@@ -9,6 +9,9 @@
#include <list>
#include <iostream>
+#include "utilities/types.h"
+
+
template<class S>
class ConeSimplex: public S
{
@@ -18,12 +21,18 @@
typedef std::list<Self> Cycle;
public:
+ ConeSimplex(const Self& s):
+ Parent(s), coned_(s.coned_) {}
ConeSimplex(const Parent& parent,
bool coned = false):
Parent(parent), coned_(coned) {}
Cycle boundary() const;
bool coned() const { return coned_; }
+ Dimension dimension() const { return coned_ ? (Parent::dimension() + 1) : Parent::dimension(); }
+
+ bool operator<(const Self& other) const { if (coned_ ^ other.coned_) return !coned_; else return Parent::operator<(other); }
+ bool operator==(const Self& other) const { return !(coned_ ^ other.coned_) && Parent::operator==(other); }
std::ostream& operator<<(std::ostream& out) const;
--- a/include/topology/cycle.h Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/cycle.h Fri Apr 04 21:14:44 2008 -0400
@@ -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/cycle.hpp Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/cycle.hpp Fri Apr 04 21:14:44 2008 -0400
@@ -166,7 +166,8 @@
{
for (const_iterator cur = begin(); cur != end(); ++cur)
{
- out << **cur << ", ";
+ if (cur != begin()) out << ", ";
+ out << **cur;
}
// out << "(last: " << *last << ")"; // For debugging only
return out;
--- a/include/topology/filtration.h Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/filtration.h Fri Apr 04 21:14:44 2008 -0400
@@ -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
@@ -61,7 +65,7 @@
/// Computes RU decomposition of the simplices in [bg, end) range, assuming that everything before bg has been paired
void pair_simplices(Index bg, Index end, bool store_trails = true);
void pair_simplices(bool store_trails = true) { pair_simplices(begin(), end(), store_trails); }
- bool transpose(Index i, bool maintain_lazy = true);
+ bool transpose(Index i, bool maintain_lazy = false);
bool is_paired() const;
Index append(const Simplex& s); ///< Appends s to the filtration
Index insert(Index prior, const Simplex& s); ///< Inserts s after prior
--- a/include/topology/filtration.hpp Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/filtration.hpp Fri Apr 04 21:14:44 2008 -0400
@@ -117,6 +117,15 @@
Index i_orig = i++;
+ rLog(rlFiltrationTranspositions, "Transposing:");
+ rLog(rlFiltrationTranspositions, " %s: (%s | %s)",
+ tostring(*i_orig).c_str(), tostring(i_orig->cycle()).c_str(),
+ tostring(i_orig->trail()).c_str());
+ rLog(rlFiltrationTranspositions, " and");
+ rLog(rlFiltrationTranspositions, " %s: (%s | %s)",
+ tostring(*i).c_str(), tostring(i->cycle()).c_str(),
+ tostring(i->trail()).c_str());
+
AssertMsg(i_orig->pair() != i, "Transposing simplices must not be paired");
bool result = transpose_simplices(i_orig, maintain_lazy);
AssertMsg(i_orig == boost::next(i), "Wrong indices after transposition");
--- a/include/topology/filtrationcontainer.h Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/filtrationcontainer.h Fri Apr 04 21:14:44 2008 -0400
@@ -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 Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/filtrationsimplex.h Fri Apr 04 21:14:44 2008 -0400
@@ -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 Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/simplex.h Fri Apr 04 21:14:44 2008 -0400
@@ -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 Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/vineyard.h Fri Apr 04 21:14:44 2008 -0400
@@ -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
@@ -32,8 +34,9 @@
typedef typename FiltrationSimplex::Simplex Simplex;
typedef Vine<Simplex> VineS;
typedef Knee<Simplex> KneeS;
- typedef std::list<VineS> VineList;
- typedef std::vector<VineList> VineListVector;
+ typedef std::list<VineS> VineList;
+ typedef std::list<VineList> VineListList;
+ typedef std::vector<typename VineListList::iterator> VineListVector;
typedef typename FiltrationSimplex::Cycle Cycle;
typedef typename FiltrationSimplex::Index Index;
@@ -45,6 +48,7 @@
void start_vines(Index bg, Index end);
void switched(Index i, Index j);
+ void record_knee(Index i);
void record_diagram(Index bg, Index end);
void set_evaluator(Evaluator* eval) { evaluator = eval; }
@@ -55,16 +59,18 @@
typename KneeS::SimplexList resolve_cycle(Index i) const;
private:
- void record_knee(Index i);
void start_vine(Index i);
private:
- VineListVector vines;
+ VineListList vines; // stores vine lists
+ VineListVector vines_vector; // stores pointers (iterators) to vine lists
Evaluator* evaluator;
};
/**
* 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
@@ -117,11 +123,17 @@
typedef typename VineRepresentation::const_iterator const_knee_iterator;
Vine() {}
- Vine(const KneeS& k) { add(k); }
+ Vine(const Vine& other):
+ VineRepresentation(other) {}
+ Vine(const VineRepresentation& other):
+ VineRepresentation(other) {}
+ Vine(const KneeS& k) { add(k); }
void add(RealType b, RealType d, RealType t) { push_back(KneeS(b,d,t)); }
void add(const KneeS& k) { push_back(k); }
+ std::ostream& operator<<(std::ostream& out) const { for (const_knee_iterator cur = begin(); cur != end(); ++cur) out << *cur; return out; }
+
using VineRepresentation::begin;
using VineRepresentation::end;
using VineRepresentation::front;
@@ -139,6 +151,9 @@
void serialize(Archive& ar, version_type );
};
+template<class S>
+std::ostream& operator<<(std::ostream& out, const Vine<S>& v) { return v.operator<<(out); }
+
#include "vineyard.hpp"
--- a/include/topology/vineyard.hpp Fri Apr 04 16:17:28 2008 -0400
+++ b/include/topology/vineyard.hpp Fri Apr 04 21:14:44 2008 -0400
@@ -23,7 +23,8 @@
if (dim >= vines.size())
{
AssertMsg(dim == vines.size(), "New dimension has to be contiguous");
- vines.push_back(std::list<VineS>());
+ vines.push_back(VineList());
+ vines_vector.push_back(boost::prior(vines.end()));
}
start_vine(cur);
@@ -60,8 +61,8 @@
AssertMsg(i->sign(), "Can only start vines for positive simplices");
Dimension dim = i->dimension();
- vines[dim].push_back(VineS());
- i->set_vine(&vines[dim].back());
+ vines_vector[dim]->push_back(VineS());
+ i->set_vine(&vines_vector[dim]->back());
i->pair()->set_vine(i->vine());
}
@@ -88,12 +89,12 @@
Vineyard<FS>::
save_edges(const std::string& filename) const
{
- for (unsigned int i = 0; i < vines.size(); ++i)
+ for (unsigned int i = 0; i < vines_vector.size(); ++i)
{
std::ostringstream os; os << i;
std::string fn = filename + os.str() + ".edg";
std::ofstream out(fn.c_str());
- for (typename VineList::const_iterator vi = vines[i].begin(); vi != vines[i].end(); ++vi)
+ for (typename VineList::const_iterator vi = vines_vector[i]->begin(); vi != vines_vector[i]->end(); ++vi)
for (typename VineS::const_iterator ki = vi->begin(), kiprev = ki++; ki != vi->end(); kiprev = ki++)
{
if (kiprev->is_infinite() || ki->is_infinite()) continue;
@@ -122,6 +123,7 @@
rLog(rlVineyard, "Creating knee");
KneeS k(evaluator->value(*i), evaluator->value(*(i->pair())), evaluator->time());
rLog(rlVineyard, "Knee created: %s", tostring(k).c_str());
+ rLog(rlVineyard, "Vine: %s", tostring(*(i->vine())).c_str());
if (!k.is_diagonal() || i->vine()->empty()) // non-diagonal k, or empty vine
{
--- a/include/utilities/consistencylist.h Fri Apr 04 16:17:28 2008 -0400
+++ b/include/utilities/consistencylist.h Fri Apr 04 21:14:44 2008 -0400
@@ -144,7 +144,7 @@
T data;
OrderType consistency;
- std::ostream& operator<<(std::ostream& out) const { return out << data << ": " << consistency; }
+ std::ostream& operator<<(std::ostream& out) const { return out << consistency << ": " << data; }
};
template<class T>
--- a/include/utilities/eventqueue.h Fri Apr 04 16:17:28 2008 -0400
+++ b/include/utilities/eventqueue.h Fri Apr 04 21:14:44 2008 -0400
@@ -8,6 +8,7 @@
#include <iostream>
#include <cassert> // TODO: switch to internal debugging system
#include <string>
+#include <algorithm>
// TODO: change inefficient list-based implementation to something heap-based
// Need a queue that supports deleting arbitrary items (given by iterator),
@@ -34,6 +35,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(); }
@@ -54,7 +58,8 @@
{
QueueRepresentation tmp;
tmp.splice(tmp.end(), queue_, i);
- iterator pos = std::find_if(queue_.begin(), queue_.end(), std::not1(std::bind2nd(EventComparison(), *i)));
+ //iterator pos = std::find_if(queue_.begin(), queue_.end(), std::not1(std::bind2nd(EventComparison(), *i)));
+ iterator pos = std::find_if(queue_.begin(), queue_.end(), std::bind1st(EventComparison(), *i));
queue_.splice(pos, tmp);
}
@@ -64,7 +69,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 Fri Apr 04 16:17:28 2008 -0400
+++ b/include/utilities/log.h Fri Apr 04 21:14:44 2008 -0400
@@ -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 Fri Apr 04 16:17:28 2008 -0400
+++ b/tests/CMakeLists.txt Fri Apr 04 21:14:44 2008 -0400
@@ -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 Fri Apr 04 21:14:44 2008 -0400
@@ -0,0 +1,16 @@
+set (targets
+ euclidean
+ test-eventqueue)
+
+if (use_synaps)
+ set (targets ${targets}
+ polynomial
+ test-kinetic-sort
+ test-linalg)
+endif (use_synaps)
+
+foreach (t ${targets})
+ add_executable (${t} ${t}.cpp)
+ target_link_libraries (${t} ${synaps_libraries} ${libraries})
+endforeach (t ${targets})
+
--- a/tests/geometry/Makefile Fri Apr 04 16:17:28 2008 -0400
+++ /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 Fri Apr 04 16:17:28 2008 -0400
+++ b/tests/geometry/euclidean.cpp Fri Apr 04 21:14:44 2008 -0400
@@ -1,4 +1,4 @@
-#include "euclidean.h"
+#include "geometry/euclidean.h"
#include <vector>
#include <iostream>
#include <cmath>
--- a/tests/geometry/polynomial.cpp Fri Apr 04 16:17:28 2008 -0400
+++ b/tests/geometry/polynomial.cpp Fri Apr 04 21:14:44 2008 -0400
@@ -1,5 +1,5 @@
-#include <euclidean.h>
-#include <polynomial.h>
+#include <geometry/euclidean.h>
+#include <geometry/polynomial.h>
#include <vector>
#include <iostream>
@@ -7,7 +7,7 @@
typedef UPolynomial<ZZ> PolynomialKernel;
typedef PolynomialKernel::Polynomial Polynomial;
-typedef PolynomialKernel::RationalFunction RationalF;
+typedef PolynomialKernel::Function RationalF;
typedef Kernel<RationalF> K;
typedef K::Point Point;
@@ -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 Fri Apr 04 16:17:28 2008 -0400
+++ b/tests/geometry/test-eventqueue.cpp Fri Apr 04 21:14:44 2008 -0400
@@ -1,4 +1,4 @@
-#include <eventqueue.h>
+#include <utilities/eventqueue.h>
#include <functional>
#include <iostream>
--- a/tests/geometry/test-kinetic-sort.cpp Fri Apr 04 16:17:28 2008 -0400
+++ b/tests/geometry/test-kinetic-sort.cpp Fri Apr 04 21:14:44 2008 -0400
@@ -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 Fri Apr 04 16:17:28 2008 -0400
+++ b/tests/geometry/test-linalg.cpp Fri Apr 04 21:14:44 2008 -0400
@@ -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;