ar-vinyeard compiles and runs:
- ar-function-kernel uses int for sign rather than bool in sign_at
- switched to dealing with a single functions themselves in value_at
- order of thresholds is handled correctly in ar-vineyard (max is the last one)
--- a/examples/ar-vineyard/ar-function-kernel.h Tue Feb 26 18:25:27 2008 -0500
+++ b/examples/ar-vineyard/ar-function-kernel.h Wed Feb 27 07:56:26 2008 -0500
@@ -35,7 +35,8 @@
ARSimplex3D* simplex() const { return simplex_; }
ARSimplex3D* simplex2() const { return simplex2_; }
- std::ostream& operator<<(std::ostream& out) const { return (out << form_ << " " << form2_); }
+ std::ostream& operator<<(std::ostream& out) const { return (out << "(" << form_ << ", " << simplex_ << "), ("
+ << form2_ << ", " << simplex2_ << ")"); }
private:
FunctionForm form_, form2_; // the function is form_ - form2_
@@ -69,7 +70,7 @@
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 int sign_at_negative_infinity(const Function& f);
static FieldType value_at(const Function& f, RootType v);
};
--- a/examples/ar-vineyard/ar-function-kernel.hpp Tue Feb 26 18:25:27 2008 -0500
+++ b/examples/ar-vineyard/ar-function-kernel.hpp Wed Feb 27 07:56:26 2008 -0500
@@ -3,6 +3,7 @@
#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
@@ -20,15 +21,23 @@
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)
- stack.push(root(s2->alpha() - s1->phi_const()));
+ {
+ 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();
@@ -44,6 +53,7 @@
// 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); }
@@ -57,6 +67,7 @@
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");
@@ -64,6 +75,8 @@
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
@@ -75,34 +88,36 @@
else return false;
}
-bool
+int
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; }
+ 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 true; // multiplier must be 1
- else return false;
+ 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;
+ return -multiplier;
if (f1 == Function::rho && f2 == Function::rho)
{
- if (s1->alpha() > s2->alpha()) return true; // multiplier must be 1
- else return false;
+ 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");
+ AssertMsg(false, "The case analysis should be exhaustive in sign at -infinity");
return false;
}
@@ -113,20 +128,17 @@
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;
+ 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 && f2 == Function::rho)
- return (v + root(s1->phi_const() - s2->alpha()))*multiplier;
+ if (f1 == Function::phi)
+ return v + root(s1->phi_const());
- if (f1 == Function::rho && f2 == Function::rho)
- return root(s1->alpha() - s2->alpha())*multiplier;
+ if (f1 == Function::rho)
+ return root(s1->alpha());
- AssertMsg(false, "The case analysis should be exhaustive");
+ AssertMsg(false, "The case analysis should be exhaustive in value_at");
return 0;
}
--- a/examples/ar-vineyard/ar-simplex3d.hpp Tue Feb 26 18:25:27 2008 -0500
+++ b/examples/ar-vineyard/ar-simplex3d.hpp Wed Feb 27 07:56:26 2008 -0500
@@ -215,7 +215,7 @@
operator<<(std::ostream& out) const
{
for (VertexSet::const_iterator cur = Parent::vertices().begin(); cur != Parent::vertices().end(); ++cur)
- out << **cur << ", ";
+ out << &(**cur) << ", ";
out << "value = " << value();
return out;
--- a/examples/ar-vineyard/ar-vineyard.cpp Tue Feb 26 18:25:27 2008 -0500
+++ b/examples/ar-vineyard/ar-vineyard.cpp Wed Feb 27 07:56:26 2008 -0500
@@ -12,7 +12,9 @@
rlog::RLogInit(argc, argv);
stderrLog.subscribeTo( RLOG_CHANNEL("error") );
- stdoutLog.subscribeTo( RLOG_CHANNEL("ar-vineyard") );
+ //stdoutLog.subscribeTo( RLOG_CHANNEL("ar/vineyard") );
+ //stdoutLog.subscribeTo( RLOG_CHANNEL("ar/function-kernel/value") );
+ //stdoutLog.subscribeTo( RLOG_CHANNEL("geometry/simulator") );
//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/filtration") );
//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/cycle") );
--- a/examples/ar-vineyard/ar-vineyard.h Tue Feb 26 18:25:27 2008 -0500
+++ b/examples/ar-vineyard/ar-vineyard.h Wed Feb 27 07:56:26 2008 -0500
@@ -46,7 +46,7 @@
ThresholdTrajectoryExtractor, Simulator> ThresholdSort;
/// @}
- typedef boost::signal<void (Simulator*)> NewFrontSignal;
+ typedef boost::signal<void (Simulator*)> NewMaxSignal;
public:
ARConeSimplex3D(const ARSimplex3D& s, bool coned = false);
@@ -55,18 +55,20 @@
thresholds_(other.thresholds_) {}
const ThresholdList& thresholds() const { return thresholds_; }
- NewFrontSignal& new_front_signal() { return new_front_signal_; }
+ NewMaxSignal& new_max_signal() { return new_max_signal_; }
+ const Function& max_threshold() const { return thresholds_.back(); }
void schedule_thresholds(Simulator* simulator);
// need explicit operator= because of the signal
ARConeSimplex3D& operator=(const ARConeSimplex3D& other) { Parent::operator=(other); thresholds_ = other.thresholds_; return *this; }
+ bool operator<(const ARConeSimplex3D& other) const { if (coned() ^ other.coned()) return !coned(); else return Parent::operator<(other); }
private:
ThresholdList thresholds_;
ThresholdSort thresholds_sort_;
- NewFrontSignal new_front_signal_;
+ NewMaxSignal new_max_signal_;
void swap_thresholds(ThresholdListIterator i, Simulator* simulator);
};
@@ -99,12 +101,12 @@
/// \name SimplexSort types
/// @{
struct SimplexTrajectoryExtractor
- { Function operator()(Index i) const { return i->thresholds().front(); }};
+ { Function operator()(Index i) const { return i->max_threshold(); } };
typedef KineticSort<Index, SimplexTrajectoryExtractor, Simulator> SimplexSort;
typedef SimplexSort::iterator SimplexSortIterator;
- class ThresholdChangeSlot; // used to notify of change in front threshold
+ class ThresholdChangeSlot; // used to notify of change in max threshold
/// @}
typedef std::list<Point> PointList;
@@ -158,7 +160,7 @@
{
public:
ThresholdChangeSlot(SimplexSortIterator iter, SimplexSort* sort):
- iter_(iter), sort_(sort) { iter_->element->new_front_signal().connect(*this); }
+ iter_(iter), sort_(sort) { iter_->element->new_max_signal().connect(*this); }
void operator()(Simulator* simulator) { sort_->update_trajectory(iter_, simulator); }
private:
@@ -182,7 +184,7 @@
simulator_(simulator) {}
virtual RealType time() const { return simulator_->current_time(); }
- virtual RealType value(const Simplex& s) const { return FunctionKernel::value_at(s.thresholds().front(), time()); }
+ virtual RealType value(const Simplex& s) const { return FunctionKernel::value_at(s.max_threshold(), time()); }
private:
Simulator* simulator_;
--- a/examples/ar-vineyard/ar-vineyard.hpp Tue Feb 26 18:25:27 2008 -0500
+++ b/examples/ar-vineyard/ar-vineyard.hpp Wed Feb 27 07:56:26 2008 -0500
@@ -5,30 +5,27 @@
#ifdef LOGGING
static rlog::RLogChannel* rlARVineyard = DEF_CHANNEL("ar/vineyard", rlog::Log_Debug);
static rlog::RLogChannel* rlARVineyardComputing = DEF_CHANNEL("ar/vineyard/computing", rlog::Log_Debug);
+static rlog::RLogChannel* rlARVineyardSwap = DEF_CHANNEL("ar/vineyard/swap", rlog::Log_Debug);
+static rlog::RLogChannel* rlARVineyardConeSimplexSwap = DEF_CHANNEL("ar/vineyard/conesimplex/swap", rlog::Log_Debug);
#endif
template <class Simulator_>
ARConeSimplex3D<Simulator_>::
ARConeSimplex3D(const ARSimplex3D& s, bool coned): Parent(s, coned)
-{
- 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
ARConeSimplex3D<Simulator_>::
swap_thresholds(ThresholdListIterator i, Simulator* simulator)
{
+ rLog(rlARVineyardConeSimplexSwap, "Transposing %s and %s", tostring(*i).c_str(),
+ tostring(*boost::next(i)).c_str());
typename ThresholdList::iterator n = boost::next(i);
thresholds_.splice(i, thresholds_, n);
- if (n == thresholds_.begin())
- new_front_signal_(simulator);
+ if (boost::next(i) == thresholds_.end())
+ new_max_signal_(simulator);
}
template <class Simulator_>
@@ -36,6 +33,13 @@
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);
}
@@ -97,12 +101,14 @@
// Once thresholds are scheduled, we can initialize the simplex_sort
simplex_sort.initialize(filtration_->begin(), filtration_->end(),
boost::bind(&ARVineyard::swap, this, _1, _2), &simulator);
+ rLog(rlARVineyardComputing, "SimplexSort initialized");
// Connect signals and slots
std::vector<ThresholdChangeSlot> slots;
slots.reserve(filtration_->size());
for (SimplexSortIterator cur = simplex_sort.begin(); cur != simplex_sort.end(); ++cur)
slots.push_back(ThresholdChangeSlot(cur, &simplex_sort));
+ rLog(rlARVineyardComputing, "Signals and slots connected");
// Simulate
change_evaluator(new KineticEvaluator(&simulator));
@@ -119,6 +125,9 @@
ARVineyard::
swap(Index i, Simulator* simulator)
{
+ rLog(rlARVineyardSwap, "Transposing %p and %p: %s and %s",
+ &(*i), &(*boost::next(i)),
+ tostring(*i).c_str(), tostring(*boost::next(i)).c_str());
filtration_->transpose(i);
}
--- a/include/geometry/simulator.hpp Tue Feb 26 18:25:27 2008 -0500
+++ b/include/geometry/simulator.hpp Wed Feb 27 07:56:26 2008 -0500
@@ -81,6 +81,7 @@
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(); }