Initial commit of KineticSort-based ARVineyard
Compiles, but needs to be runtime debugged.
--- a/examples/ar-vineyard/CMakeLists.txt Mon Feb 25 04:33:13 2008 -0500
+++ b/examples/ar-vineyard/CMakeLists.txt Mon Feb 25 04:37:49 2008 -0500
@@ -1,8 +1,12 @@
set (targets
ar-vineyard)
-
+
+find_library (boost_signals_LIBRARY NAMES boost_signals
+ PATHS ${Boost_LIBRARY_DIR})
+
add_definitions (${cgal_cxxflags})
foreach (t ${targets})
add_executable (${t} ${t}.cpp ${external_sources})
- target_link_libraries (${t} ${libraries} ${cgal_libraries} ${synaps_libraries})
+ target_link_libraries (${t} ${libraries} ${cgal_libraries}
+ ${synaps_libraries} ${boost_signals_LIBRARY})
endforeach (t)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/ar-vineyard/ar-function-kernel.h Mon Feb 25 04:37:49 2008 -0500
@@ -0,0 +1,79 @@
+#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_ << " " << form2_); }
+
+ 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 bool 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 Mon Feb 25 04:37:49 2008 -0500
@@ -0,0 +1,132 @@
+#include <utilities/log.h>
+#include <cmath>
+
+#if LOGGING
+static rlog::RLogChannel* rlARFunctionKernel = DEF_CHANNEL("ar/function-kernel", 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");
+
+ //if (f1 == Function::phi && f2 == Function::phi) return;
+ //if (f1 == Function::rho && f2 == Function::rho) return;
+
+ if (f1 == Function::phi && f2 == Function::rho)
+ stack.push(root(s2->alpha() - s1->phi_const()));
+
+ if (f1 == Function::phi && f2 == Function::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)
+ {
+ 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)
+ {
+ // 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());
+ }
+}
+
+int
+ARFunctionKernel::
+sign_at(const Function& f, RootType r)
+{
+ FieldType v = value_at(f,r);
+ if (v > 0) return true;
+ else return false;
+}
+
+bool
+ARFunctionKernel::
+sign_at_negative_infinity(const Function& f)
+{
+ FunctionForm f1 = f.form(), f2 = f.form2();
+ const ARSimplex3D *s1 = f.simplex(),
+ *s2 = f.simplex2();
+ bool multiplier = true;
+ if (f1 < f2) { std::swap(f1, f2); std::swap(s1, s2); multiplier = false; }
+
+ 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 true; // multiplier must be 1
+ else return false;
+ }
+
+ if (f1 == Function::phi && f2 == Function::rho)
+ return !multiplier;
+
+ if (f1 == Function::rho && f2 == Function::rho)
+ {
+ if (s1->alpha() > s2->alpha()) return true; // multiplier must be 1
+ else return false;
+ }
+
+ AssertMsg(false, "The case analysis should be exhaustive");
+ 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();
+ 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)
+ return root(s1->phi_const() - s2->phi_const())*multiplier;
+
+ if (f1 == Function::phi && f2 == Function::rho)
+ return (v + root(s1->phi_const() - s2->alpha()))*multiplier;
+
+ if (f1 == Function::rho && f2 == Function::rho)
+ return root(s1->alpha() - s2->alpha())*multiplier;
+
+ AssertMsg(false, "The case analysis should be exhaustive");
+ return 0;
+}
--- a/examples/ar-vineyard/ar-simplex3d.h Mon Feb 25 04:33:13 2008 -0500
+++ b/examples/ar-vineyard/ar-simplex3d.h Mon Feb 25 04:37:49 2008 -0500
@@ -28,7 +28,6 @@
typedef Delaunay::Facet Facet;
typedef Delaunay::Cell Cell;
typedef Delaunay::Cell_handle Cell_handle;
-typedef K::FT RealValue;
typedef Delaunay::Finite_vertices_iterator Vertex_iterator;
typedef Delaunay::Finite_edges_iterator Edge_iterator;
@@ -40,6 +39,7 @@
class ARSimplex3D: public SimplexWithVertices<Vertex_handle>
{
public:
+ typedef K::FT RealValue;
typedef std::map<ARSimplex3D, RealValue> SimplexPhiMap;
typedef SimplexWithVertices<Vertex_handle> Parent;
typedef Parent::VertexContainer VertexSet;
@@ -89,6 +89,8 @@
RealValue phi_const_; // see LHI paper, Appendices A and B
bool attached_;
+ // in paper's notation: s_ = v^2; v_ = d^2
+
};
typedef std::vector<ARSimplex3D> ARSimplex3DVector;
--- a/examples/ar-vineyard/ar-simplex3d.hpp Mon Feb 25 04:33:13 2008 -0500
+++ b/examples/ar-vineyard/ar-simplex3d.hpp Mon Feb 25 04:37:49 2008 -0500
@@ -1,3 +1,9 @@
+#include <utilities/log.h>
+
+#if LOGGING
+static rlog::RLogChannel* rlARSimplex3D = DEF_CHANNEL("ar/simplex3d", rlog::Log_Debug);
+#endif
+
ARSimplex3D::
ARSimplex3D(const ARSimplex3D& s): Parent(s)
{
@@ -226,19 +232,19 @@
ARSimplex3D::SimplexPhiMap simplices;
for(Cell_iterator cur = Dt.finite_cells_begin(); cur != Dt.finite_cells_end(); ++cur)
update_simplex_phi_map(ARSimplex3D(*cur, z), simplices);
- std::cout << "Cells inserted" << std::endl;
+ rLog(rlARSimplex3D, "Cells inserted");
for(Vertex_iterator cur = Dt.finite_vertices_begin(); cur != Dt.finite_vertices_end(); ++cur)
simplices[ARSimplex3D(*cur, z)] = 0; // only one tetrahedron can have non-negative phi_const value
// (namely the one containing z); all other simplices will have a
// negative phi_const value, so 0 is safe
- std::cout << "Vertices inserted" << std::endl;
+ rLog(rlARSimplex3D, "Vertices inserted");
for(Facet_iterator cur = Dt.finite_facets_begin(); cur != Dt.finite_facets_end(); ++cur)
update_simplex_phi_map(ARSimplex3D(*cur, z, simplices), simplices);
- std::cout << "Facets inserted" << std::endl;
+ rLog(rlARSimplex3D, "Facets inserted");
for(Edge_iterator cur = Dt.finite_edges_begin(); cur != Dt.finite_edges_end(); ++cur)
update_simplex_phi_map(ARSimplex3D(*cur, z, simplices, Dt.incident_facets(*cur)), simplices);
- std::cout << "Edges inserted" << std::endl;
+ rLog(rlARSimplex3D, "Edges inserted");
// Sort simplices by their alpha values
alpha_order.resize(simplices.size()); ARSimplex3DVector::iterator out = alpha_order.begin();
--- a/examples/ar-vineyard/ar-vineyard.cpp Mon Feb 25 04:33:13 2008 -0500
+++ b/examples/ar-vineyard/ar-vineyard.cpp Mon Feb 25 04:37:49 2008 -0500
@@ -11,7 +11,9 @@
#ifdef LOGGING
rlog::RLogInit(argc, argv);
- stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
+ stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+ stdoutLog.subscribeTo( RLOG_CHANNEL("ar-vineyard") );
+
//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/filtration") );
//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/cycle") );
//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/vineyard") );
@@ -52,7 +54,7 @@
arv.compute_pairing();
// Compute vineyard
- arv.compute_vineyard(true);
+ arv.compute_vineyard();
std::cout << "Vineyard computed" << std::endl;
arv.vineyard()->save_edges(outfilename);
}
--- a/examples/ar-vineyard/ar-vineyard.h Mon Feb 25 04:33:13 2008 -0500
+++ b/examples/ar-vineyard/ar-vineyard.h Mon Feb 25 04:37:49 2008 -0500
@@ -6,81 +6,106 @@
#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 <list>
#include "ar-simplex3d.h"
+#include "ar-function-kernel.h"
-class ARConeSimplex: public ConeSimplex<ARSimplex3D>
+template <class Simulator_>
+class ARConeSimplex3D: public ConeSimplex<ARSimplex3D>
{
public:
typedef ConeSimplex<ARSimplex3D> Parent;
typedef ARSimplex3D ARSimplex3D;
- typedef Filtration<ARConeSimplex> Filtration;
- typedef Filtration::Index Index;
- /// \name Polynomial Kernel types
+ /// \name Simulator types
/// @{
- typedef double FieldType;
- typedef UPolynomial<FieldType> PolyKernel;
- typedef PolyKernel::Polynomial Polynomial;
- typedef Simulator<PolyKernel> Simulator;
+ typedef Simulator_ Simulator;
+ typedef typename Simulator::FunctionKernel FunctionKernel;
+ typedef typename FunctionKernel::Function Function;
+ /// @}
- typedef KineticSort<Index, SimplexTrajectoryExtractor, Simulator>
- SimplexSort;
- typedef SimplexSort::iterator SimplexSortIterator;
- typedef SimplexSortIterator Key;
+ /// \name ThresholdSort types
+ /// @{
+ typedef std::list<Function> ThresholdList;
+ typedef typename ThresholdList::iterator ThresholdListIterator;
+
+ struct ThresholdTrajectoryExtractor
+ { Function operator()(ThresholdListIterator i) const { return *i; } };
+
+ typedef KineticSort<ThresholdListIterator,
+ ThresholdTrajectoryExtractor, Simulator> ThresholdSort;
/// @}
- /// \name Kinetic Sort types
- /// @{
- typedef std::list<Polynomial> ThresholdList;
-
- struct ThresholdTrajectoryExtractor
- { Polynomial operator()(ThresholdList::iterator i) const { return *i; } }
- struct SimplexTrajectoryExtractor
- { Polynomial operator()(Index i) const { i->thresholds().front(); }
+ typedef boost::signal<void (Simulator*)> NewFrontSignal;
+
+ public:
+ ARConeSimplex3D(const ARSimplex3D& s, bool coned = false);
+ ARConeSimplex3D(const ARConeSimplex3D& other): // need explicit copy-constructor because of the signal
+ Parent(other, other.coned()),
+ thresholds_(other.thresholds_) {}
- typedef KineticSort<ThresholdList::iterator, ThresholdTrajectoryExtractor, Simulator>
- ThresholdSort;
- /// @}
+ const ThresholdList& thresholds() const { return thresholds_; }
+ NewFrontSignal& new_front_signal() { return new_front_signal_; }
- ARConeSimplex(const ARSimplex3D& s, bool coned = false):
- Parent(s, coned) {}
+ void schedule_thresholds(Simulator* simulator);
- Key kinetic_key() const { return key_; }
- void set_kinetic_key(Key k) { key_ = k; }
- ThresholdList& thresholds() { return thresholds_; }
-
- void schedule_thresholds(SimplexSort* sort, Simulator* simulator);
+ // 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_;
+ NewFrontSignal new_front_signal_;
- void swap_thresholds(SimplexSort* sort, ThresholdList::iterator i, Simulator* simulator);
+ void swap_thresholds(ThresholdListIterator i, Simulator* simulator);
};
-
+/**
+ * Encapsulated filtration, and provides compute_vineyard() functionality.
+ */
class ARVineyard
{
public:
typedef ARVineyard Self;
+
+ /// \name FunctionKernel and Simulator types
+ /// @{
+ typedef ARFunctionKernel FunctionKernel;
+ typedef FunctionKernel::Function Function;
+ typedef Simulator<FunctionKernel> Simulator;
+ /// @}
+
+ /// \name Filtration types
+ /// @{
+ typedef ARConeSimplex3D<Simulator> ARConeSimplex3D;
+ typedef Filtration<ARConeSimplex3D> Filtration;
+ typedef Filtration::Simplex Simplex;
+ typedef Filtration::Index Index;
+ typedef Filtration::Vineyard Vineyard;
+ typedef Vineyard::Evaluator Evaluator;
+ /// @}
- typedef ARConeSimplex::Filtration ARFiltration;
- typedef ARFiltration::Simplex Simplex;
- typedef ARFiltration::Index Index;
- typedef ARFiltration::Vineyard Vineyard;
- typedef Vineyard::Evaluator Evaluator;
+ /// \name SimplexSort types
+ /// @{
+ struct SimplexTrajectoryExtractor
+ { Function operator()(Index i) const { return i->thresholds().front(); }};
+
+ typedef KineticSort<Index, SimplexTrajectoryExtractor, Simulator> SimplexSort;
+ typedef SimplexSort::iterator SimplexSortIterator;
- typedef ARConeSimplex::Simulator Simulator;
- typedef ARConeSimplex::SimplexSort SimplexSort;
-
+ class ThresholdChangeSlot; // used to notify of change in front threshold
+ /// @}
typedef std::list<Point> PointList;
@@ -92,9 +117,9 @@
~ARVineyard();
void compute_pairing();
- void compute_vineyard(bool explicit_events = false);
+ void compute_vineyard();
- const ARFiltration* filtration() const { return filtration_; }
+ const Filtration* filtration() const { return filtration_; }
const Vineyard* vineyard() const { return vineyard_; }
public:
@@ -105,7 +130,7 @@
void change_evaluator(Evaluator* eval);
private:
- ARFiltration* filtration_;
+ Filtration* filtration_;
Vineyard* vineyard_;
Evaluator* evaluator_;
@@ -129,35 +154,38 @@
//BOOST_CLASS_EXPORT(ARVineyard)
+class ARVineyard::ThresholdChangeSlot
+{
+ public:
+ ThresholdChangeSlot(SimplexSortIterator iter, SimplexSort* sort):
+ iter_(iter), sort_(sort) { iter_->element->new_front_signal().connect(*this); }
+ void operator()(Simulator* simulator) { sort_->update_trajectory(iter_, simulator); }
+
+ private:
+ SimplexSortIterator iter_;
+ SimplexSort* sort_;
+};
class ARVineyard::StaticEvaluator: public Evaluator
{
public:
- StaticEvaluator(RealType t): time_(t) {}
+ StaticEvaluator() {}
- virtual RealType time() const { return time_; }
+ virtual RealType time() const { return 0; }
virtual RealType value(const Simplex& s) const { return s.value(); }
-
- private:
- RealType time_;
};
class ARVineyard::KineticEvaluator: public Evaluator
{
public:
- KineticEvaluator(Simulator::Handle sp,
- ActivePointsTable::Handle apt,
- RealType time_offset):
- sp_(sp), apt_(apt) {}
+ KineticEvaluator(Simulator* simulator):
+ simulator_(simulator) {}
- virtual RealType time() const { return CGAL::to_double(get_time()); }
- virtual RealType value(const Simplex& s) const { return CGAL::to_double(apt_->at(s.kinetic_key()).x()(get_time())); }
+ virtual RealType time() const { return simulator_->current_time(); }
+ virtual RealType value(const Simplex& s) const { return FunctionKernel::value_at(s.thresholds().front(), time()); }
private:
- Simulator::Time get_time() const { return sp_->current_time(); }
-
- Simulator::Handle sp_;
- ActivePointsTable::Handle apt_;
+ Simulator* simulator_;
};
--- a/examples/ar-vineyard/ar-vineyard.hpp Mon Feb 25 04:33:13 2008 -0500
+++ b/examples/ar-vineyard/ar-vineyard.hpp Mon Feb 25 04:37:49 2008 -0500
@@ -1,45 +1,65 @@
#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);
+#endif
+
+
+template <class Simulator_>
+ARConeSimplex3D<Simulator_>::
+ARConeSimplex3D(const ARSimplex3D& s, bool coned): Parent(s, coned)
+{
+ if (!coned) thresholds_.push_back(Function(Function::rho, this));
+ else
+ {
+ thresholds_.push_back(Function(Function::rho, this));
+ thresholds_.push_back(Function(Function::phi, this));
+ }
+}
+template <class Simulator_>
void
-ARConeSimplex::
-swap_thresholds(SimplexSort* sort, ThresholdList::iterator i, Simulator* simulator)
+ARConeSimplex3D<Simulator_>::
+swap_thresholds(ThresholdListIterator i, Simulator* simulator)
{
typename ThresholdList::iterator n = boost::next(i);
- tl->splice(i, *tl, n);
- if (n == tl->begin())
- sort->update_trajectory(kinetic_key(), simulator);
+ thresholds_.splice(i, thresholds_, n);
+ if (n == thresholds_.begin())
+ new_front_signal_(simulator);
}
+template <class Simulator_>
void
-ARConeSimplex::
-schedule_thresholds(SimplexSort* sort, Simulator* simulator)
+ARConeSimplex3D<Simulator_>::
+schedule_thresholds(Simulator* simulator)
{
thresholds_sort_.initialize(thresholds_.begin(), thresholds_.end(),
- boost::bind(&ARConeSimplex::swap_thresholds, this, sort, _1, _2), simulator);
+ boost::bind(&ARConeSimplex3D::swap_thresholds, this, _1, _2), simulator);
}
ARVineyard::
-ARVineyard(const PointList& points, const Point& z): simplex_sort_(0), z_(z)
+ARVineyard(const PointList& points, const Point& z): z_(z)
{
for (PointList::const_iterator cur = points.begin(); cur != points.end(); ++cur)
dt_.insert(*cur);
- std::cout << "Delaunay triangulation computed" << std::endl;
+ rLog(rlARVineyard, "Delaunay triangulation computed");
ARSimplex3DVector alpha_ordering;
fill_alpha_order(dt_, z_, alpha_ordering);
- std::cout << "Delaunay simplices: " << alpha_ordering.size() << std::endl;
+ rLog(rlARVineyard, "Delaunay simplices: %i", alpha_ordering.size());
- evaluator_ = new StaticEvaluator(0);
+ evaluator_ = new StaticEvaluator;
vineyard_ = new Vineyard(evaluator_);
- filtration_ = new ARFiltration(vineyard_);
+ filtration_ = new Filtration(vineyard_);
for (ARSimplex3DVector::const_iterator cur = alpha_ordering.begin(); cur != alpha_ordering.end(); ++cur)
{
- filtration_->append(ARConeSimplex(*cur)); // Delaunay simplex
- filtration_->append(ARConeSimplex(*cur, true)); // Coned off delaunay simplex
+ filtration_->append(ARConeSimplex3D(*cur)); // Delaunay simplex
+ filtration_->append(ARConeSimplex3D(*cur, true)); // Coned off delaunay simplex
}
}
@@ -58,96 +78,48 @@
filtration_->fill_simplex_index_map();
filtration_->pair_simplices(filtration_->begin(), filtration_->end());
vineyard_->start_vines(filtration_->begin(), filtration_->end());
- std::cout << "Simplices paired" << std::endl;
+ rLog(rlARVineyard, "Simplices paired");
}
void
ARVineyard::
-compute_vineyard(bool explicit_events)
+compute_vineyard()
{
AssertMsg(filtration_->is_paired(), "Simplices must be paired for a vineyard to be computed");
Simulator simulator;
SimplexSort simplex_sort;
- // Set thresholds
+ // Schedule thresholds
for (Index cur = filtration_->begin(); cur != filtration_->end(); ++cur)
- {
- cur->thresholds.push_back(ARConeSimplex::Polynomial(CGAL::to_double(cur->alpha())));
-
- if (!cur->coned()) continue; // non-coned simplices stay put, so we are done
-
- 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());
+ cur->schedule_thresholds(&simulator);
- 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;
-
- cur->set_kinetic_key();
- }
-
-
- // Once thresholds are set (and sorted), we can initialize the simplex_sort
- simplex_sort.initialize(filtration_.begin(), filtration_.end(),
+ // Once thresholds are scheduled, we can initialize the simplex_sort
+ simplex_sort.initialize(filtration_->begin(), filtration_->end(),
boost::bind(&ARVineyard::swap, this, _1, _2), &simulator);
-
+ // 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));
- // Process all the events (compute the vineyard in the process)
- // FIXME: the time should not be 1, but something like twice the radius of
- // the pointset as seen from z
- change_evaluator(new KineticEvaluator(sp, apt, 0));
- if (explicit_events)
- {
- while (sp->next_event_time() < 1)
- {
- std::cout << "Next event time: " << sp->next_event_time() << std::endl;
- sp->set_current_event_number(sp->current_event_number() + 1);
- std::cout << "Processed event" << std::endl;
- }
- } else
- sp->set_current_time(1.0);
- std::cout << "Processed " << sp->current_event_number() << " events" << std::endl;
+ // Simulate
+ change_evaluator(new KineticEvaluator(&simulator));
+ while(!simulator.reached_infinity())
+ {
+ rLog(rlARVineyardComputing, "Current time before: %lf", simulator.current_time());
+ simulator.process();
+ }
- //change_evaluator(new StaticEvaluator(1));
vineyard_->record_diagram(filtration_->begin(), filtration_->end());
}
void
ARVineyard::
-swap(Key a, Key b)
+swap(Index i, Simulator* simulator)
{
- Index ao = kinetic_map_[a], bo = kinetic_map_[b];
- AssertMsg(filtration_->get_trails_cmp()(ao, bo), "In swap(a,b), a must precede b");
- filtration_->transpose(ao);
- AssertMsg(filtration_->get_trails_cmp()(bo, ao), "In swap(a,b), b must precede a after the transposition");
+ filtration_->transpose(i);
}
void
@@ -160,30 +132,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/geometry/polynomial.hpp Mon Feb 25 04:33:13 2008 -0500
+++ b/include/geometry/polynomial.hpp Mon Feb 25 04:37:49 2008 -0500
@@ -1,3 +1,5 @@
+#include <utilities/log.h>
+
template<class T>
void
UPolynomial<T>::
@@ -5,6 +7,8 @@
{
typedef SYNAPS::Seq<RootType> RootSeq;
+ 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());
--- a/include/geometry/simulator.h Mon Feb 25 04:33:13 2008 -0500
+++ b/include/geometry/simulator.h Mon Feb 25 04:37:49 2008 -0500
@@ -2,7 +2,6 @@
#define __SIMULATOR_H__
#include "utilities/eventqueue.h"
-#include "polynomial.h"
template<class Comparison> class IndirectComparison;