Fixed a serious bug in DynamicPersistenceTrails + LSVineyard sorts filtration by attachment + miscellanea:
* added program options to pdbdistance-vineyard
* pdbdistance-vineyard takes frame list in a file
* more logging
* verification of pairing after vertex transpositions (commented out)
--- a/examples/pl-functions/CMakeLists.txt Sun Dec 20 10:43:00 2009 -0800
+++ b/examples/pl-functions/CMakeLists.txt Fri Dec 25 06:56:27 2009 -0800
@@ -15,7 +15,7 @@
foreach (t ${targets})
add_executable (${t} ${t}.cpp)
- target_link_libraries (${t} ${libraries} ${CGAL_LIBRARY} ${dsrpdb_LIBRARY})
+ target_link_libraries (${t} ${libraries} ${CGAL_LIBRARY} ${dsrpdb_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
endforeach (t ${targets})
# Extra special for program_options
--- a/examples/pl-functions/pdbdistance-vineyard.cpp Sun Dec 20 10:43:00 2009 -0800
+++ b/examples/pl-functions/pdbdistance-vineyard.cpp Fri Dec 25 06:56:27 2009 -0800
@@ -8,50 +8,32 @@
#include <string>
#include <sstream>
-std::string frame_filename(const std::string& prefix, int frame, int subframe)
-{
- std::ostringstream os;
- os << prefix << frame << "_" << subframe << ".pdb";
- return os.str();
-}
+#include <boost/program_options.hpp>
+
+void program_options(int argc, char* argv[], std::string& input_fn_list,
+ std::string& output_prefix,
+ bool& all_atoms,
+ bool& save_vines);
int main(int argc, char** argv)
{
#ifdef LOGGING
rlog::RLogInit(argc, argv);
-
- stdoutLog.subscribeTo( RLOG_CHANNEL("topology/filtration") );
- stdoutLog.subscribeTo( RLOG_CHANNEL("topology/cycle") );
- stdoutLog.subscribeTo( RLOG_CHANNEL("topology/vineyard") );
- stdoutLog.subscribeTo( RLOG_CHANNEL("topology/lowerstar") );
#endif
- if (argc < 5)
- {
- std::cout << "Usage: pdbdistance FILENAME LASTFRAME LASTSUBFRAME OUTFILENAME [CAs_ONLY]" << std::endl;
- std::cout << " FILENAME - prefix of the filenames of the PDB frames" << std::endl;
- std::cout << " LASTFRAME - the last frame number" << std::endl;
- std::cout << " LASTSUBFRAME - the last subframe number" << std::endl;
- std::cout << " OUTFILENAME - filename prefix for the resulting vineyards" << std::endl;
- std::cout << " CAs_ONLY - only use alpha carbons [1 = true, 0 = false, default: 1]" << std::endl;
- std::cout << std::endl;
- std::cout << "Computes a vineyard of the pairwise distance function for a sequence of PDB frames." << std::endl;
- std::cout << "Frames are in files FILENAME#1_#2.pdb, where #1 goes from 0 to LASTFRAME, " << std::endl;
- std::cout << "and #2 goes from 0 to LASTSUBFRAME." << std::endl;
- exit(0);
- }
- std::string infilename = argv[1];
- int lastframe; std::istringstream(argv[2]) >> lastframe;
- int lastsubframe; std::istringstream(argv[3]) >> lastsubframe;
- std::string outfilename = argv[4];
- bool cas_only = true;
- if (argc > 5)
- std::istringstream(argv[5]) >> cas_only;
+ std::string input_fn_list, output_fn;
+ bool all_atoms = false, save_vines = false;
+ program_options(argc, argv, input_fn_list, output_fn, all_atoms, save_vines);
// Compute initial filtration
- int f = 0; int sf = 0;
- std::ifstream in(frame_filename(infilename, f, sf++).c_str());
- PDBDistanceGrid ginit(in, cas_only);
+ std::vector<std::string> frame_fns;
+ std::ifstream in_fns(input_fn_list.c_str());
+ std::string fn;
+ while(std::getline(in_fns, fn))
+ frame_fns.push_back(fn);
+
+ std::ifstream in(frame_fns[0].c_str());
+ PDBDistanceGrid ginit(in, !all_atoms);
in.close();
typedef LSVineyard<Grid2D::CoordinateIndex, Grid2D> Grid2DVineyard;
@@ -68,23 +50,75 @@
std::cout << "Pairing computed" << std::endl;
// Process frames computing the vineyard
- while (f <= lastframe)
+ for (size_t i = 1; i < frame_fns.size(); ++i)
{
- std::string fn = frame_filename(infilename, f, sf++);
- std::cout << "Processing " << fn << std::endl;
- in.open(fn.c_str());
- v.compute_vineyard(PDBDistanceGrid(in, cas_only));
+ std::cout << "Processing " << frame_fns[i] << std::endl;
+ in.open(frame_fns[i].c_str());
+ v.compute_vineyard(PDBDistanceGrid(in, !all_atoms));
in.close();
- if (sf == lastsubframe) { sf = 0; ++f; }
}
std::cout << "Vineyard computed" << std::endl;
- v.vineyard().save_edges(outfilename);
+ if (save_vines)
+ v.vineyard().save_vines(output_fn);
+ else
+ v.vineyard().save_edges(output_fn);
#if 0
- std::ofstream ofs(outfilename.c_str(), std::ios::binary);
+ std::ofstream ofs(output_fn.c_str(), std::ios::binary);
boost::archive::binary_oarchive oa(ofs);
oa << make_nvp("Filtration", pgf);
ofs.close();
#endif
}
+
+void program_options(int argc, char* argv[], std::string& input_fn_list,
+ std::string& output_prefix,
+ bool& all_atoms,
+ bool& save_vines)
+{
+ namespace po = boost::program_options;
+
+ po::options_description hidden("Hidden options");
+ hidden.add_options()
+ ("input-fn-list", po::value<std::string>(&input_fn_list), "prefix of the input frames")
+ ("output-prefix", po::value<std::string>(&output_prefix), "output prefix");
+
+ po::options_description visible("Allowed options", 100);
+ visible.add_options()
+ ("all-atoms,a", po::bool_switch(&all_atoms), "process all atoms (not only alpha carbons)")
+ ("save-vines,v", po::bool_switch(&save_vines), "save vines instead of edges")
+ ("help,h", "produce help message");
+#if LOGGING
+ std::vector<std::string> log_channels;
+ visible.add_options()
+ ("log,l", po::value< std::vector<std::string> >(&log_channels), "log channels to turn on (info, debug, etc)");
+#endif
+
+ po::positional_options_description pos;
+ pos.add("input-fn-list", 1);
+ pos.add("output-prefix", 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(pos).run(), vm);
+ po::notify(vm);
+
+#if LOGGING
+ for (std::vector<std::string>::const_iterator cur = log_channels.begin(); cur != log_channels.end(); ++cur)
+ stderrLog.subscribeTo( RLOG_CHANNEL(cur->c_str()) );
+#endif
+
+ if (vm.count("help") || !vm.count("input-fn-list") || !vm.count("output-prefix"))
+ {
+ std::cout << "Usage: " << argv[0] << " [options] input-fn-list output-prefix" << std::endl;
+ std::cout << visible << std::endl;
+ std::cout << std::endl;
+ std::cout << "Computes a vineyard of the pairwise distance function for a sequence of PDB frames." << std::endl;
+ std::cout << "Frames are listed in input-fn-list file." << std::endl;
+ std::abort();
+ }
+}
+
--- a/include/topology/dynamic-persistence.h Sun Dec 20 10:43:00 2009 -0800
+++ b/include/topology/dynamic-persistence.h Fri Dec 25 06:56:27 2009 -0800
@@ -74,6 +74,7 @@
typedef typename Parent::ContainerTraits Traits;
typedef typename Parent::Order Order;
+ typedef typename Traits::ConsistentContainer Consistency;
typedef typename Parent::OrderComparison OrderComparison;
typedef typename Parent::OrderIndex OrderIndex;
typedef ConsistencyComparison_ ConsistencyComparison;
@@ -107,6 +108,10 @@
using Parent::iterator_to;
using Parent::index;
using Parent::size;
+ using Parent::order_comparison;
+
+ template<class Iter>
+ void rearrange(Iter i);
// Struct: TranspositionVisitor
//
@@ -131,6 +136,9 @@
using Parent::set_pair;
using Parent::swap_cycle;
+ Consistency& consistent_order() { return order().get<consistency>(); }
+ const Consistency& consistent_order() const { return order().get<consistency>(); }
+
bool trail_remove_if_contains
(iterator i, OrderIndex j) { TrailRemover rm(j); order().modify(i, rm); return rm.result; }
void cycle_add(iterator i, const Cycle& z) { order().modify(i, boost::bind(&Element::template cycle_add<ConsistencyComparison>, bl::_1, boost::ref(z), ccmp_)); } // i->cycle_add(z, ccmp_)
@@ -219,6 +227,7 @@
typedef typename Parent::ContainerTraits Traits;
typedef typename Parent::Order Order;
+
typedef typename Parent::OrderComparison OrderComparison;
typedef typename Parent::OrderIndex OrderIndex;
typedef ConsistencyComparison_ ConsistencyComparison;
@@ -276,7 +285,7 @@
using Parent::order;
using Parent::set_pair;
using Parent::swap_cycle;
-
+
void cycle_add(iterator i, const Cycle& z) { order().modify(i, boost::bind(&Element::template cycle_add<ConsistencyComparison>, bl::_1, boost::ref(z), ccmp_)); } // i->cycle_add(z, ccmp_)
void chain_add(iterator i, const Chain& c) { order().modify(i, boost::bind(&Element::template chain_add<ConsistencyComparison>, bl::_1, boost::ref(c), ccmp_)); } // i->chain_add(c, ccmp_)
--- a/include/topology/dynamic-persistence.hpp Sun Dec 20 10:43:00 2009 -0800
+++ b/include/topology/dynamic-persistence.hpp Fri Dec 25 06:56:27 2009 -0800
@@ -26,7 +26,7 @@
template<class Filtration>
DynamicPersistenceTrails<D,CT,OT,E,Cmp,CCmp>::
DynamicPersistenceTrails(const Filtration& f):
- Parent(f), ccmp_(order().get<consistency>())
+ Parent(f), ccmp_(consistent_order())
{}
template<class D, class CT, class OT, class E, class Cmp, class CCmp>
@@ -113,6 +113,7 @@
if (!(l->cycle.contains(index(i_prev))))
{
// Case 1.2
+ rLog(rlTranspositions, "k is in l: %d", (bool) l->trail.contains(index(k))); // if true, a special update would be needed to maintain lazy decomposition
swap(i_prev, i);
rLog(rlTranspositions, "Case 1.2");
Count(cTranspositionCase12);
@@ -120,7 +121,7 @@
} else
{
// Case 1.1
- if (std::not2(ccmp_)(index(k),index(l)))
+ if (std::not2(order_comparison())(index(k),index(l)))
{
// Case 1.1.1
swap(i_prev, i);
@@ -160,7 +161,7 @@
trail_add(i_prev, i->trail); // Add row i to i_prev
cycle_add(i, i_prev->cycle); // Add column i_prev to i
swap(i_prev, i);
- if (std::not2(ccmp_)(index(low_ii), index(low_i)))
+ if (std::not2(order_comparison())(index(low_ii), index(low_i)))
{
// Case 2.1.2
cycle_add(i_prev, i->cycle); // Add column i to i_prev (after transposition)
@@ -215,6 +216,26 @@
return false; // to avoid compiler complaints; we should never reach this point
}
+
+template<class D, class CT, class OT, class E, class Cmp, class CCmp>
+template<class Iter>
+void
+DynamicPersistenceTrails<D,CT,OT,E,Cmp,CCmp>::
+rearrange(Iter i)
+{
+ order().rearrange(i);
+ consistent_order().rearrange(i);
+
+ // Resort the cycles
+ Cycle z;
+ for(iterator i = begin(); i != end(); ++i)
+ {
+ Parent::swap_cycle(i, z);
+ z.sort(ccmp_);
+ Parent::swap_cycle(i, z);
+ }
+}
+
template<class D, class CT, class OT, class E, class Cmp, class CCmp>
void
DynamicPersistenceTrails<D,CT,OT,E,Cmp,CCmp>::
--- a/include/topology/filtration.h Sun Dec 20 10:43:00 2009 -0800
+++ b/include/topology/filtration.h Fri Dec 25 06:56:27 2009 -0800
@@ -71,6 +71,8 @@
void push_back(const Simplex& s) { container_.get<order>().push_back(s); }
void transpose(Index i) { container_.get<order>().relocate(i, i+1); }
void clear() { container_.get<order>().clear(); }
+ template<class Iter>
+ void rearrange(Iter i) { container_.get<order>().rearrange(i); }
Index begin() const { return container_.get<order>().begin(); }
Index end() const { return container_.get<order>().end(); }
--- a/include/topology/lsvineyard.h Sun Dec 20 10:43:00 2009 -0800
+++ b/include/topology/lsvineyard.h Fri Dec 25 06:56:27 2009 -0800
@@ -19,7 +19,8 @@
#include <CGAL/Kinetic/Sort.h>
#include <CGAL/Kinetic/Sort_visitor_base.h>
-#include <list>
+#include <boost/tuple/tuple.hpp>
+namespace b = boost;
template<class Vertex_, class VertexEvaluator_, class Simplex_ = Simplex<Vertex_>, class Filtration_ = Filtration<Simplex_> >
@@ -36,10 +37,37 @@
typedef Filtration_ LSFiltration;
typedef typename LSFiltration::Index LSFIndex;
- class KineticVertexType;
+ 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;
+ class KineticVertexType
+ {
+ public:
+ KineticVertexType(const Vertex& v):
+ vertex_(v) {}
+
+ Key kinetic_key() const { return key_; }
+ void set_kinetic_key(Key k) { key_ = k; }
+
+ Vertex vertex() const { return vertex_; }
+ void set_vertex(Vertex v) { vertex_ = v; }
+
+ LSFIndex simplex_index() const { return simplex_index_; }
+ void set_simplex_index(LSFIndex i) { simplex_index_ = i; }
+
+ private:
+ Key key_;
+ Vertex vertex_;
+ LSFIndex simplex_index_;
+ };
class KineticVertexComparison;
- typedef std::list<KineticVertexType> VertexContainer;
+ typedef typename OrderContainer<KineticVertexType>::Container
+ VertexContainer;
typedef typename VertexContainer::iterator VertexIndex;
+ typedef std::map<Key, VertexIndex> KeyVertexMap;
struct AttachmentData: public VineData
{
@@ -58,7 +86,9 @@
class KineticEvaluator;
class DimensionFromIterator;
+ typedef std::map<Vertex, LSFIndex> VertexLSFIndexMap;
typedef ThroughEvaluatorComparison<VertexEvaluator> VertexComparison;
+ class VertexAttachmentComparison;
typedef MaxVertexComparison<Simplex, VertexComparison> SimplexComparison;
class TranspositionVisitor;
@@ -66,14 +96,6 @@
typedef Vineyard<Index, iterator, Evaluator> Vnrd;
- 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 std::map<Key, VertexIndex> KeyVertexMap;
-
public:
template<class VertexIterator>
LSVineyard(VertexIterator begin, VertexIterator end,
@@ -95,6 +117,7 @@
VertexValue simplex_value(const Simplex& s) const { return vertex_value(*std::max_element(s.vertices().begin(), s.vertices().end(), vcmp_)); }
const Simplex& pfmap(iterator i) const { return pfmap_[i]; }
const Simplex& pfmap(Index i) const { return pfmap_[i]; }
+ VertexIndex filtration_attachment(LSFIndex i) const { return (persistence().begin() + (i - filtration().begin()))->attachment; }
Index index(iterator i) const { return persistence_.index(i); }
@@ -107,6 +130,13 @@
void set_attachment(iterator i, VertexIndex vi) { persistence_.modifier()(i, boost::bind(&AttachmentData::set_attachment, bl::_1, vi)); }
void transpose_filtration(iterator i) { filtration_.transpose(filtration_.begin() + (i - persistence_.begin())); }
+ bool verify_pairing() const;
+
+ typedef b::tuple< b::reference_wrapper<const Simplex>,
+ b::reference_wrapper<const typename Persistence::Element> >
+ SimplexPersistenceElementTuple;
+ struct AttachmentCmp;
+
private:
VertexContainer vertices_;
@@ -143,27 +173,27 @@
//BOOST_CLASS_EXPORT(LSVineyard)
-template<class V, class VE, class S, class C>
-class LSVineyard<V,VE,S,C>::KineticVertexType
-{
- public:
- KineticVertexType(const Vertex& v):
- vertex_(v) {}
+// template<class V, class VE, class S, class C>
+// class LSVineyard<V,VE,S,C>::KineticVertexType
+// {
+// public:
+// KineticVertexType(const Vertex& v):
+// vertex_(v) {}
- Key kinetic_key() const { return key_; }
- void set_kinetic_key(Key k) { key_ = k; }
+// Key kinetic_key() const { return key_; }
+// void set_kinetic_key(Key k) { key_ = k; }
- Vertex vertex() const { return vertex_; }
- void set_vertex(Vertex v) { vertex_ = v; }
+// Vertex vertex() const { return vertex_; }
+// void set_vertex(Vertex v) { vertex_ = v; }
- iterator simplex_index() const { return simplex_index_; }
- void set_simplex_index(iterator i) { simplex_index_ = i; }
+// LSFIndex simplex_index() const { return simplex_index_; }
+// void set_simplex_index(iterator i) { simplex_index_ = i; }
- private:
- Key key_;
- Vertex vertex_;
- iterator simplex_index_;
-};
+// private:
+// Key key_;
+// Vertex vertex_;
+// LSFIndex simplex_index_;
+// };
template<class V, class VE, class S, class C>
std::ostream&
--- a/include/topology/lsvineyard.hpp Sun Dec 20 10:43:00 2009 -0800
+++ b/include/topology/lsvineyard.hpp Fri Dec 25 06:56:27 2009 -0800
@@ -1,5 +1,14 @@
#include <utilities/log.h>
+#include <boost/function.hpp>
+#include <boost/bind.hpp>
+#include <boost/lambda/lambda.hpp>
+namespace bl = boost::lambda;
+
+#include <boost/iterator/zip_iterator.hpp>
+#include <boost/iterator/transform_iterator.hpp>
+#include <boost/foreach.hpp>
+
#ifdef LOGGING
static rlog::RLogChannel* rlLSVineyard = DEF_CHANNEL("lsvineyard/info", rlog::Log_Debug);
static rlog::RLogChannel* rlLSVineyardDebug = DEF_CHANNEL("lsvineyard/debug", rlog::Log_Debug);
@@ -15,44 +24,78 @@
template<class VertexIterator>
LSVineyard<V,VE,S,F>::
LSVineyard(VertexIterator begin, VertexIterator end,
- LSFiltration& filtration,
+ LSFiltration& fltr,
const VertexEvaluator& veval):
- vertices_(begin, end), filtration_(filtration),
+ filtration_(fltr),
+ vertices_(begin, end),
persistence_(filtration_),
veval_(veval), vcmp_(veval_), scmp_(vcmp_),
pfmap_(persistence_.make_simplex_map(filtration_)),
time_count_(0)
{
- // TODO: LSVineyard really should sort the filtration_ itself
- // filtration_.sort(scmp_);
- // persistence_.init(filtration_);
-
- rLog(rlLSVineyardDebug, "Initializing LSVineyard");
- persistence_.pair_simplices();
- rLog(rlLSVineyardDebug, "Simplices paired");
-
vertices_.sort(KineticVertexComparison(vcmp_)); // sort vertices w.r.t. vcmp_
#if LOGGING
rLog(rlLSVineyardDebug, "Vertex order:");
- for (VertexIndex cur = vertices_.begin(); cur != vertices_.end(); ++cur)
+ for (typename VertexContainer::iterator cur = vertices_.begin(); cur != vertices_.end(); ++cur)
rLog(rlLSVineyardDebug, " %d", cur->vertex());
#endif
- // Initialize simplex_index() and attachment
- VertexIndex vi = boost::prior(vertices_.begin());
- for (iterator i = persistence_.begin(); i != persistence_.end(); ++i)
+ // Record Vertex -> LSFIndex map
+ VertexLSFIndexMap vimap;
+ for (LSFIndex i = filtration().begin(); i != filtration().end(); ++i)
{
- const Simplex& s = pfmap_[i];
+ const Simplex& s = *i;
+ rLog(rlLSVineyardDebug, "Simplex: %s", tostring(*i).c_str());
if (s.dimension() == 0)
- {
- ++vi;
- AssertMsg(s.vertices().front() == vi->vertex(), "In constructor, simplices and vertices must match.");
- vi->set_simplex_index(i);
- }
- set_attachment(i, vi);
- rLog(rlLSVineyardDebug, "%s attached to %d", tostring(pfmap(i)).c_str(), vi->vertex());
+ vimap[s.vertices().front()] = i;
+ }
+
+ // Assign vertex attachments and simplex_index
+ OffsetMap<LSFIndex, iterator> fpmap(filtration().begin(), persistence().begin());
+ for (typename VertexContainer::iterator vi = vertices_.begin(); vi != vertices_.end(); ++vi)
+ {
+ LSFIndex i = vimap[vi->vertex()];
+ const Simplex& s = *i;
+ AssertMsg(s.vertices().front() == vi->vertex(), "In constructor, simplices and vertices must match.");
+ vertices_.modify(vi, b::bind(&KineticVertexType::set_simplex_index, bl::_1, i)); // vi->set_simplex_index(i)
+ set_attachment(fpmap[i], vi);
+ rLog(rlLSVineyardDebug, "%s attached to %d", tostring(*i).c_str(), vi->vertex());
}
+ // Assign attachments for all the simplices
+ VertexAttachmentComparison vacmp(vimap, *this);
+ for (LSFIndex i = filtration().begin(); i != filtration().end(); ++i)
+ set_attachment(fpmap[i], fpmap[vimap[*std::max_element(i->vertices().begin(), i->vertices().end(), vacmp)]]->attachment);
+
+ // Order filtration_ and persistence_ based on attachment
+ rLog(rlLSVineyardDebug, "Ordering the simplices");
+
+ std::vector<SimplexPersistenceElementTuple> fporder
+ (b::make_zip_iterator(b::make_tuple(filtration().begin(), persistence().begin())),
+ b::make_zip_iterator(b::make_tuple(filtration().end(), persistence().end())));
+ std::sort(fporder.begin(), fporder.end(), AttachmentCmp());
+
+ // Rearrage filtration
+ std::vector< b::reference_wrapper<const Simplex> > sv;
+ BOOST_FOREACH(const SimplexPersistenceElementTuple& t, fporder) sv.push_back(b::get<0>(t));
+ filtration_.rearrange (sv.begin());
+
+ // Rearrange persistence
+ std::vector< b::reference_wrapper<const typename Persistence::Element> > pev;
+ BOOST_FOREACH(const SimplexPersistenceElementTuple& t, fporder) pev.push_back(b::get<1>(t));
+ persistence_.rearrange(pev.begin());
+
+#if LOGGING
+ rLog(rlLSVineyardDebug, "Simplices:");
+ for(iterator i = persistence().begin(); i != persistence().end(); ++i)
+ rLog(rlLSVineyardDebug, " %s attached to %d", tostring(pfmap(i)).c_str(), i->attachment->vertex());
+#endif
+
+ // Pair simplices
+ rLog(rlLSVineyardDebug, "Initializing LSVineyard");
+ persistence_.pair_simplices();
+ rLog(rlLSVineyardDebug, "Simplices paired");
+
evaluator_ = new StaticEvaluator(*this, time_count_);
vineyard_.set_evaluator(evaluator_);
vineyard_.start_vines(persistence_.begin(), persistence_.end());
@@ -90,7 +133,7 @@
rLog(rlLSVineyardDebug, "Vertex %d: %f -> %f", cur->vertex(), val0, val1);
F x = cf(F::NT(val0), F::NT(val1 - val0)); // x = val0 + (val1 - val0)*t = (1-t)*val0 + t*val1
Point p(x);
- cur->set_kinetic_key(apt->insert(p));
+ vertices_.modify(cur, b::bind(&KineticVertexType::set_kinetic_key, bl::_1, apt->insert(p))); // cur->set_kinetic_key(apt->insert(p))
kinetic_map_[cur->kinetic_key()] = cur;
rLog(rlLSVineyardDebug, "Scheduling: %d %s", cur->vertex(), tostring(x).c_str());
}
@@ -111,7 +154,7 @@
veval_ = veval;
change_evaluator(new StaticEvaluator(*this, ++time_count_));
- vineyard_.record_diagram(persistence_.begin(), persistence_.end());
+ vineyard_.record_diagram(persistence().begin(), persistence().end());
}
template<class V, class VE, class S, class F>
@@ -138,54 +181,56 @@
vineyard_.set_evaluator(evaluator_);
}
-#include <boost/function.hpp>
-#include <boost/bind.hpp>
-#include <boost/lambda/lambda.hpp>
-namespace bl = boost::lambda;
-
template<class V, class VE, class S, class F>
bool
LSVineyard<V,VE,S,F>::
transpose_vertices(VertexIndex vi)
{
Count(cVertexTransposition);
- rLog(rlLSVineyard, "Transposing vertices (%d, %d)", vi->vertex(), boost::next(vi)->vertex());
+ rLog(rlLSVineyard, "Transposing vertices (%d:%d, %d:%d)", vi->vertex(), (vi - vertices_.begin()),
+ b::next(vi)->vertex(), (b::next(vi) - vertices_.begin()));
DimensionFromIterator dim(pfmap_);
TranspositionVisitor visitor(*this);
- iterator i = vi->simplex_index();
- iterator i_prev = boost::prior(i);
- iterator i_next = boost::next(vi)->simplex_index();
- iterator i_next_prev = boost::prior(i_next); // transpositions are done in terms of the first index in the pair
- iterator j = boost::next(i_next);
+ OffsetMap<LSFIndex, iterator> fpmap(filtration().begin(), persistence().begin());
+ iterator i = fpmap[vi->simplex_index()];
+ iterator i_prev = b::prior(i);
+ iterator i_next = fpmap[b::next(vi)->simplex_index()];
+ iterator i_next_prev = b::prior(i_next); // transpositions are done in terms of the first index in the pair
+ iterator j = b::next(i_next);
- VertexIndex vi_next = boost::next(vi);
+ VertexIndex vi_next = b::next(vi);
const Vertex& v = vi->vertex();
bool result = false; // has a switch in pairing occurred
- // First, move the vertex --- this can be sped up if we devise special "vertex transpose" operation
+ // First move the vertex --- this can be sped up if we devise special "vertex transpose" operation
rLog(rlLSVineyardDebug, "Starting to move the vertex");
while (i_next_prev != i_prev)
{
rLog(rlLSVineyardDebug, " Transposing %s %s", tostring(pfmap(i_next_prev)).c_str(),
- tostring(pfmap(boost::next(i_next_prev))).c_str());
+ tostring(pfmap(b::next(i_next_prev))).c_str());
result |= persistence_.transpose(i_next_prev, dim, visitor);
- i_next_prev = boost::prior(i_next);
+ AssertMsg((i_next_prev <= persistence().iterator_to(i_next_prev->pair)) == i_next_prev->sign(), "Pairing must respect order");
+ AssertMsg((i_next <= persistence().iterator_to(i_next->pair)) == i_next->sign(), "Pairing must respect order");
+ i_next_prev = b::prior(i_next);
}
rLog(rlLSVineyardDebug, "Done moving the vertex");
// Second, move the simplices attached to it
rLog(rlLSVineyardDebug, "Moving attached simplices");
+ // rLog(rlLSVineyardDebug, " Considering %s", tostring(pfmap(j)).c_str());
+ // rLog(rlLSVineyardDebug, " attachment %d", j->attachment->vertex());
while (j != persistence_.end() && j->attachment == vi_next)
{
rLog(rlLSVineyardDebug, " Considering %s", tostring(pfmap(j)).c_str());
if (pfmap(j).contains(v)) // j becomes attached to v and does not move
{
Count(cAttachment);
- rLog(rlLSVineyardDebug, " Attachment changed for ", tostring(pfmap(j)).c_str());
+ rLog(rlLSVineyardDebug, " Attachment changed for %s to %d", tostring(pfmap(j)).c_str(), vi->vertex());
set_attachment(j, vi);
+ AssertMsg(fpmap[vi->simplex_index()] < j, "The simplex must be attached to a preceding vertex");
++j;
continue;
}
@@ -195,18 +240,58 @@
{
rLog(rlLSVineyardDebug, " Moving: %s, %s",
tostring(pfmap(j_prev)).c_str(),
- tostring(pfmap(boost::next(j_prev))).c_str());
+ tostring(pfmap(b::next(j_prev))).c_str());
AssertMsg(j_prev->attachment == vi, "Simplex preceding the one being moved must be attached to v");
result |= persistence_.transpose(j_prev, dim, visitor);
+ AssertMsg((j_prev <= persistence().iterator_to(j_prev->pair)) == j_prev->sign(), "Pairing must respect order");
+ AssertMsg((j <= persistence().iterator_to(j->pair)) == j->sign(), "Pairing must respect order");
--j_prev;
}
}
rLog(rlLSVineyard, "Done moving attached simplices");
- vertices_.splice(vi, vertices_, vi_next); // swap vi and vi_next
+ vertices_.relocate(vi, vi_next); // swap vi and vi_next
+#if LSVINEYARD_CONSISTENCY
+ AssertMsg(verify_pairing(), "Pairing must be correct after vertex transposition");
+#endif
+
return result;
}
+template<class V, class VE, class S, class F>
+bool
+LSVineyard<V,VE,S,F>::
+verify_pairing() const
+{
+ rLog(rlLSVineyardDebug, "Verifying pairing");
+ StaticPersistence<> p(filtration());
+ p.pair_simplices();
+ iterator i = persistence().begin();
+ StaticPersistence<>::iterator ip = p.begin();
+ StaticPersistence<>::SimplexMap<LSFiltration> m = p.make_simplex_map(filtration());
+
+ while (ip != p.end())
+ {
+ if (&pfmap(i) != &m[ip])
+ {
+ rError("DP: %s %s", tostring(pfmap(i)).c_str(), tostring(pfmap(i->pair)).c_str());
+ rError("SP: %s %s", tostring(m[ip]).c_str(), tostring(m[ip->pair]).c_str());
+ rError("The order must match");
+ return false;
+ }
+ if (&pfmap(i->pair) != &m[ip->pair])
+ {
+ rError("DP: %s %s", tostring(pfmap(i)).c_str(), tostring(pfmap(i->pair)).c_str());
+ rError("SP: %s %s", tostring(m[ip]).c_str(), tostring(m[ip->pair]).c_str());
+ rError("The pairing must match");
+ return false;
+ }
+ ++i; ++ip;
+ }
+
+ return true;
+}
+
/* Evaluators */
template<class V, class VE, class S, class C>
@@ -252,3 +337,35 @@
ActivePointsTable::Handle apt_;
RealType time_offset_;
};
+
+
+template<class V, class VE, class S, class C>
+class LSVineyard<V,VE,S,C>::VertexAttachmentComparison:
+ public std::binary_function<Vertex, Vertex, bool>
+{
+ public:
+ VertexAttachmentComparison(const VertexLSFIndexMap& vimap,
+ const LSVineyard& vnrd):
+ vimap_(vimap), vnrd_(vnrd) {}
+ bool operator()(Vertex v1, Vertex v2) const
+ { return vnrd_.filtration_attachment(vimap_.find(v1)->second) < vnrd_.filtration_attachment(vimap_.find(v2)->second); }
+
+ private:
+ const VertexLSFIndexMap& vimap_;
+ const LSVineyard& vnrd_;
+};
+
+
+template<class V, class VE, class S, class C>
+struct LSVineyard<V,VE,S,C>::AttachmentCmp:
+ public std::binary_function<const SimplexPersistenceElementTuple&, const SimplexPersistenceElementTuple&, bool>
+{
+ bool operator()(const SimplexPersistenceElementTuple& t1, const SimplexPersistenceElementTuple& t2) const
+ {
+ if (b::get<1>(t1).get().attachment == b::get<1>(t2).get().attachment)
+ return b::get<0>(t1).get().dimension() < b::get<0>(t2).get().dimension();
+ else
+ return b::get<1>(t1).get().attachment < b::get<1>(t2).get().attachment;
+ }
+};
+