--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/consistency/rips-consistency-zigzag.cpp Tue Mar 24 15:17:39 2009 -0700
@@ -0,0 +1,338 @@
+#include <topology/rips.h>
+#include <topology/zigzag-persistence.h>
+#include <utilities/types.h>
+#include <utilities/containers.h>
+
+#include <utilities/log.h>
+#include <utilities/memory.h> // for report_memory()
+#include <utilities/timer.h>
+
+#include <map>
+#include <cmath>
+#include <fstream>
+#include <sstream>
+#include <stack>
+#include <cstdlib>
+
+#include <boost/tuple/tuple.hpp>
+#include <boost/program_options.hpp>
+#include <boost/progress.hpp>
+
+#include "../rips/l2distance.h" // Point, PointContainer, L2DistanceType
+
+#ifdef COUNTERS
+static Counter* cComplexSize = GetCounter("rips/size");
+static Counter* cOperations = GetCounter("rips/operations");
+#endif // COUNTERS
+
+typedef PairwiseDistances<PointContainer, L2Distance> PairDistances;
+typedef PairDistances::DistanceType DistanceType;
+
+typedef PairDistances::IndexType Vertex;
+typedef std::vector<Vertex> SubsampleVector;
+
+typedef Simplex<Vertex> Smplx;
+typedef std::vector<Smplx> SimplexVector;
+typedef std::set<Smplx, Smplx::VertexDimensionComparison> SimplexSet;
+
+typedef std::vector<Vertex> VertexVector;
+typedef std::vector<DistanceType> EpsilonVector;
+typedef std::vector<std::pair<Vertex, Vertex> > EdgeVector;
+
+typedef Rips<PairDistances, Smplx> RipsGenerator;
+typedef RipsGenerator::Evaluator SimplexEvaluator;
+
+struct BirthInfo;
+typedef ZigzagPersistence<BirthInfo> Zigzag;
+typedef Zigzag::SimplexIndex Index;
+typedef Zigzag::Death Death;
+typedef std::map<Smplx, Index,
+ Smplx::VertexDimensionComparison> Complex;
+typedef Zigzag::ZColumn Boundary;
+
+// Information we need to know when a class dies
+struct BirthInfo
+{
+ BirthInfo(Dimension d = 0, unsigned i = 0, bool u = false):
+ dimension(d), index(i), un(u) {}
+
+ bool operator<(const BirthInfo& other) const { if (index == other.index) return (!un && other.un); else return index < other.index; }
+ bool operator>(const BirthInfo& other) const { return other.operator<(*this); }
+
+ Dimension dimension;
+ unsigned index;
+ bool un;
+};
+std::ostream& operator<<(std::ostream& out, const BirthInfo& bi)
+{ out << "(d=" << bi.dimension << ") " << bi.index; if (bi.un) out << " U " << (bi.index + 1); return out; }
+
+// Forward declarations of auxilliary functions
+void report_death(std::ostream& out, const BirthInfo& birth, const BirthInfo& death);
+void make_boundary(const Smplx& s, Complex& c, const Zigzag& zz, Boundary& b);
+void process_command_line_options(int argc,
+ char* argv[],
+ unsigned& ambient_dimension,
+ unsigned& skeleton_dimension,
+ float& epsilon,
+ std::string& points_filename,
+ std::string& subsample_filename,
+ std::string& outfilename);
+void add_simplex(const Smplx& s, const BirthInfo& birth, const BirthInfo& death,
+ Complex& complex, Zigzag& zz, std::ostream& out, Timer& add);
+void remove_simplex(const Smplx& s, const BirthInfo& birth, const BirthInfo& death,
+ Complex& complex, Zigzag& zz, std::ostream& out, Timer& remove);
+
+int main(int argc, char* argv[])
+{
+#ifdef LOGGING
+ rlog::RLogInit(argc, argv);
+
+ stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+#endif
+
+ Timer total, remove, add, coface, generate;
+ total.start();
+
+#if 0
+ SetFrequency(cOperations, 25000);
+ SetTrigger(cOperations, cComplexSize);
+#endif
+
+ unsigned ambient_dimension;
+ unsigned skeleton_dimension;
+ float epsilon;
+ std::string points_filename, subsample_filename, outfilename;
+ process_command_line_options(argc, argv, ambient_dimension, skeleton_dimension, epsilon,
+ points_filename, subsample_filename, outfilename);
+
+ // Read in points
+ PointContainer points;
+ read_points(points_filename, points, ambient_dimension);
+
+ // Read in subsamples
+ std::ifstream subsample_in(subsample_filename.c_str());
+ std::vector<SubsampleVector> subsamples;
+ subsamples.push_back(SubsampleVector()); // Empty subsample to start from
+ while(subsample_in)
+ {
+ std::string line; std::getline(subsample_in, line);
+ std::istringstream line_in(line);
+ subsamples.push_back(SubsampleVector());
+ unsigned sample;
+ while (line_in >> sample)
+ subsamples.back().push_back(sample);
+ }
+
+ std::cout << "Subsample size:" << std::endl;
+ for (unsigned i = 0; i < subsamples.size(); ++i)
+ std::cout << " " << subsamples[i].size() << std::endl;
+
+ // Create output file
+ std::ofstream out(outfilename.c_str());
+
+ // Create pairwise distances
+ PairDistances distances(points);
+
+
+ // Construct zigzag
+ Complex complex;
+ Zigzag zz;
+ RipsGenerator rips(distances);
+ SimplexEvaluator size(distances);
+ SimplexVector subcomplex, across;
+
+ rInfo("Commencing computation");
+ boost::progress_display show_progress(subsamples.size());
+ for (unsigned i = 0; i < subsamples.size() - 1; ++i)
+ {
+ // Take union of subsamples i and i+1
+ SubsampleVector forward_difference, backward_difference;
+ std::set_difference(subsamples[i+1].begin(), subsamples[i+1].end(),
+ subsamples[i].begin(), subsamples[i].end(),
+ std::back_inserter(forward_difference));
+ std::set_difference(subsamples[i].begin(), subsamples[i].end(),
+ subsamples[i+1].begin(), subsamples[i+1].end(),
+ std::back_inserter(backward_difference));
+ rInfo("Forward difference size: %d", forward_difference.size());
+ rInfo("Backward difference size: %d", backward_difference.size());
+
+#if LOGGING
+ rDebug("Forward difference:");
+ for (SubsampleVector::const_iterator cur = forward_difference.begin(); cur != forward_difference.end(); ++cur)
+ rDebug(" %d", *cur);
+
+ rDebug("Backward difference:");
+ for (SubsampleVector::const_iterator cur = backward_difference.begin(); cur != backward_difference.end(); ++cur)
+ rDebug(" %d", *cur);
+#endif
+
+ // Add simplices: take star of the vertices in subsamples[i+1] - subsamples[i]
+ // by first computing the Rips complex on those vertices, and then
+ // taking the star of each individual simplex (perhaps not the most
+ // efficient procedure, but it will do for the first pass)
+ rInfo("Adding");
+ subcomplex.clear(); across.clear();
+ generate.start();
+ rips.generate(skeleton_dimension, epsilon, make_push_back_functor(subcomplex), forward_difference.begin(), forward_difference.end());
+ std::sort(subcomplex.begin(), subcomplex.end(), Smplx::VertexDimensionComparison());
+ generate.stop();
+ rInfo(" Subcomplex size: %d", subcomplex.size());
+ for (SimplexVector::const_iterator cur = subcomplex.begin(); cur != subcomplex.end(); ++cur)
+ {
+ add_simplex(*cur, BirthInfo(cur->dimension(), i, true), BirthInfo(cur->dimension() - 1, i),
+ complex, zz, out, add);
+
+ // Record cofaces with all other vertices outside of the subcomplex
+ coface.start();
+ rips.cofaces(*cur, skeleton_dimension, epsilon, make_push_back_functor(across), subsamples[i].begin(), subsamples[i].end());
+ coface.stop();
+ }
+ std::sort(across.begin(), across.end(), Smplx::VertexDimensionComparison());
+ rInfo(" Subcomplex simplices added");
+ rInfo(" Cross simplices size: %d", across.size());
+
+ // Add simplices that cut across VR(K_i) and VR(K_{i+1})
+ for (SimplexVector::const_iterator cur = across.begin(); cur != across.end(); ++cur)
+ add_simplex(*cur, BirthInfo(cur->dimension(), i, true), BirthInfo(cur->dimension() - 1, i),
+ complex, zz, out, add);
+ rInfo(" Cross simplices added");
+
+
+ // Remove simplices: take star of the vertices in subsamples[i] - subsamples[i+1]
+ rInfo("Removing");
+ subcomplex.clear(); across.clear();
+ generate.start();
+ rips.generate(skeleton_dimension, epsilon, make_push_back_functor(subcomplex), backward_difference.begin(), backward_difference.end());
+ std::sort(subcomplex.begin(), subcomplex.end(), Smplx::VertexDimensionComparison());
+ generate.stop();
+ rInfo(" Subcomplex size: %d", subcomplex.size());
+
+ for (SimplexVector::const_iterator cur = subcomplex.begin(); cur != subcomplex.end(); ++cur)
+ rips.cofaces(*cur, skeleton_dimension, epsilon, make_push_back_functor(across), subsamples[i+1].begin(), subsamples[i+1].end());
+ std::sort(across.begin(), across.end(), Smplx::VertexDimensionComparison());
+ rInfo(" Cross simplices size: %d", across.size());
+
+ for (SimplexVector::const_reverse_iterator cur = across.rbegin(); cur != across.rend(); ++cur)
+ remove_simplex(*cur, BirthInfo(cur->dimension() - 1, i+1), BirthInfo(cur->dimension(), i, true),
+ complex, zz, out, remove);
+ rInfo(" Cross simplices removed");
+
+ for (SimplexVector::const_reverse_iterator cur = subcomplex.rbegin(); cur != subcomplex.rend(); ++cur)
+ remove_simplex(*cur, BirthInfo(cur->dimension() - 1, i+1), BirthInfo(cur->dimension(), i, true),
+ complex, zz, out, remove);
+ rInfo(" Subcomplex simplices removed");
+
+ rInfo("Complex size: %d", complex.size());
+
+ ++show_progress;
+ }
+
+ std::cout << std::endl;
+ total.stop();
+ remove.check("Remove timer ");
+ add.check ("Add timer ");
+ coface.check("Coface timer ");
+ total.check ("Total timer ");
+}
+
+
+void make_boundary(const Smplx& s, Complex& c, const Zigzag& zz, Boundary& b)
+{
+ // rDebug(" Boundary of <%s>", tostring(s).c_str());
+ for (Smplx::BoundaryIterator cur = s.boundary_begin(); cur != s.boundary_end(); ++cur)
+ {
+ b.append(c[*cur], zz.cmp);
+ // rDebug(" %d", c[*cur]->order);
+ }
+}
+
+void report_death(std::ostream& out, const BirthInfo& birth, const BirthInfo& death)
+{
+ if (death > birth)
+ out << birth << " --- " << death << std::endl;
+}
+
+void add_simplex(const Smplx& s, const BirthInfo& birth, const BirthInfo& death,
+ Complex& complex, Zigzag& zz, std::ostream& out, Timer& add)
+{
+ rDebug("Adding simplex: %s", tostring(s).c_str());
+
+ Index idx; Death d; Boundary b;
+ make_boundary(s, complex, zz, b);
+ add.start();
+ boost::tie(idx, d) = zz.add(b, birth);
+ if (d) report_death(out, *d, death);
+ add.stop();
+ complex.insert(std::make_pair(s, idx));
+}
+
+void remove_simplex(const Smplx& s, const BirthInfo& birth, const BirthInfo& death,
+ Complex& complex, Zigzag& zz, std::ostream& out, Timer& remove)
+{
+ rDebug("Removing simplex: %s", tostring(s).c_str());
+
+ Complex::iterator si = complex.find(s);
+ remove.start();
+ Death d = zz.remove(si->second, birth);
+ remove.stop();
+ if (d) report_death(out, *d, death);
+ complex.erase(si);
+}
+
+void process_command_line_options(int argc,
+ char* argv[],
+ unsigned& ambient_dimension,
+ unsigned& skeleton_dimension,
+ float& epsilon,
+ std::string& points_filename,
+ std::string& subsample_filename,
+ std::string& outfilename)
+{
+ namespace po = boost::program_options;
+
+ po::options_description hidden("Hidden options");
+ hidden.add_options()
+ ("points-file", po::value<std::string>(&points_filename), "Point set whose Rips consistency zigzag we want to compute")
+ ("subsample-file", po::value<std::string>(&subsample_filename),"Subsample list")
+ ("output-file", po::value<std::string>(&outfilename), "Location to save persistence pairs");
+
+ po::options_description visible("Allowed options", 100);
+ visible.add_options()
+ ("help,h", "produce help message")
+ ("ambient-dimsnion,a", po::value<unsigned>(&ambient_dimension)->default_value(3), "The ambient dimension of the point set")
+ ("skeleton-dimsnion,s", po::value<unsigned>(&skeleton_dimension)->default_value(2), "Dimension of the Rips complex we want to compute")
+ ("epsilon,e", po::value<float>(&epsilon)->default_value(0), "epsilon to use when computing the Rips complex");
+#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("points-file", 1);
+ pos.add("subsample-file", 1);
+ pos.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(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()) );
+ /**
+ * Interesting channels
+ * "info", "debug", "topology/persistence"
+ */
+#endif
+
+ if (vm.count("help") || !vm.count("points-file") || !vm.count("subsample-file") || !vm.count("output-file"))
+ {
+ std::cout << "Usage: " << argv[0] << " [options] points-file subsample-file output-file" << std::endl;
+ std::cout << visible << std::endl;
+ std::abort();
+ }
+}
--- a/include/topology/rips.h Mon Mar 23 11:40:29 2009 -0700
+++ b/include/topology/rips.h Tue Mar 24 15:17:39 2009 -0700
@@ -53,8 +53,14 @@
void edge_cofaces(IndexType u, IndexType v, Dimension k, DistanceType max, const Functor& f,
Iterator candidates_begin, Iterator candidates_end) const;
+ // Calls functor f on all the simplices of the Rips complex that contain the given Simplex s
+ // (unlike the previous methods it does not call the functor on the Simplex s itself)
+ template<class Functor, class Iterator>
+ void cofaces(const Simplex& s, Dimension k, DistanceType max, const Functor& f,
+ Iterator candidates_begin, Iterator candidates_end) const;
- // No Iterator argument means IndexType and distances().begin() - distances().end()
+
+ /* No Iterator argument means Iterator = IndexType and the range is [distances().begin(), distances().end()) */
template<class Functor>
void generate(Dimension k, DistanceType max, const Functor& f) const
{ generate(k, max, f, boost::make_counting_iterator(distances().begin()), boost::make_counting_iterator(distances().end())); }
@@ -67,6 +73,10 @@
void edge_cofaces(IndexType u, IndexType v, Dimension k, DistanceType max, const Functor& f) const
{ edge_cofaces(u, v, k, max, f, boost::make_counting_iterator(distances().begin()), boost::make_counting_iterator(distances().end())); }
+ template<class Functor>
+ void cofaces(const Simplex& s, Dimension k, DistanceType max, const Functor& f) const
+ { cofaces(s, k, max, f, boost::make_counting_iterator(distances().begin()), boost::make_counting_iterator(distances().end())); }
+
const Distances& distances() const { return distances_; }
DistanceType max_distance() const;
@@ -83,7 +93,8 @@
typename VertexContainer::const_iterator excluded,
Dimension max_dim,
const NeighborTest& neighbor,
- const Functor& functor) const;
+ const Functor& functor,
+ bool check_initial = true) const;
private:
const Distances& distances_;
--- a/include/topology/rips.hpp Mon Mar 23 11:40:29 2009 -0700
+++ b/include/topology/rips.hpp Tue Mar 24 15:17:39 2009 -0700
@@ -4,6 +4,7 @@
#include <iostream>
#include <utilities/log.h>
#include <utilities/counter.h>
+#include <utilities/indirect.h>
#include <boost/iterator/counting_iterator.hpp>
#include <functional>
@@ -77,6 +78,42 @@
bron_kerbosch(current, candidates, boost::prior(candidates.begin()), k, neighbor, f);
}
+template<class D, class S>
+template<class Functor, class Iterator>
+void
+Rips<D,S>::
+cofaces(const Simplex& s, Dimension k, DistanceType max, const Functor& f, Iterator bg, Iterator end) const
+{
+ rLog(rlRipsDebug, "In cofaces(%s)", tostring(s).c_str());
+
+ WithinDistance neighbor(distances(), max);
+
+ // current = s.vertices()
+ VertexContainer current(s.vertices().begin(), s.vertices().end());
+
+ // candidates = everything - s.vertices() that is a neighbor() of every vertex in the simplex
+ VertexContainer candidates;
+ typedef difference_iterator<Iterator,
+ typename VertexContainer::const_iterator,
+ std::less<Vertex> > DifferenceIterator;
+ for (DifferenceIterator cur = DifferenceIterator(bg, end, s.vertices().begin(), s.vertices().end());
+ cur != DifferenceIterator(end, end, s.vertices().end(), s.vertices().end());
+ ++cur)
+ {
+ bool nghbr = true;
+ for (typename VertexContainer::const_iterator v = s.vertices().begin(); v != s.vertices().end(); ++v)
+ if (!neighbor(*v, *cur)) { nghbr = false; break; }
+
+ if (nghbr)
+ {
+ candidates.push_back(*cur);
+ rLog(rlRipsDebug, " added candidate: %d", *cur);
+ }
+ }
+
+ bron_kerbosch(current, candidates, boost::prior(candidates.begin()), k, neighbor, f, false);
+}
+
template<class D, class S>
template<class Functor, class NeighborTest>
@@ -87,11 +124,12 @@
typename VertexContainer::const_iterator excluded,
Dimension max_dim,
const NeighborTest& neighbor,
- const Functor& functor) const
+ const Functor& functor,
+ bool check_initial) const
{
rLog(rlRipsDebug, "Entered bron_kerbosch");
- if (!current.empty())
+ if (check_initial && !current.empty())
{
Simplex s(current);
rLog(rlRipsDebug, "Reporting simplex: %s", tostring(s).c_str());