Updated Rips zigzags
* Updated to use Bron-Kerbosch
* Compute in the order of increasing epsilon, and decreasing point sizes
* Split into plain zigzag (rips-zigzag) and image zigzag (rips-image-zigzag)
* Added minor enhancements (show_progress and timers)
--- a/examples/rips/CMakeLists.txt Sat Jan 31 23:02:22 2009 -0800
+++ b/examples/rips/CMakeLists.txt Tue Feb 03 11:31:10 2009 -0800
@@ -1,5 +1,6 @@
set (targets
rips
+ rips-image-zigzag
rips-zigzag)
foreach (t ${targets})
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/rips/rips-image-zigzag.cpp Tue Feb 03 11:31:10 2009 -0800
@@ -0,0 +1,478 @@
+#include <topology/rips.h>
+#include <topology/image-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 <stack>
+#include <cstdlib>
+
+#include <boost/tuple/tuple.hpp>
+#include <boost/program_options.hpp>
+#include <boost/progress.hpp>
+
+
+#ifdef COUNTERS
+static Counter* cComplexSize = GetCounter("rips/size");
+static Counter* cOperations = GetCounter("rips/operations");
+#endif // COUNTERS
+
+typedef std::vector<double> Point;
+typedef std::vector<Point> PointContainer;
+struct L2Distance:
+ public std::binary_function<const Point&, const Point&, double>
+{
+ result_type operator()(const Point& p1, const Point& p2) const
+ {
+ AssertMsg(p1.size() == p2.size(), "Points must be in the same dimension (in L2Distance)");
+ result_type sum = 0;
+ for (size_t i = 0; i < p1.size(); ++i)
+ sum += (p1[i] - p2[i])*(p1[i] - p2[i]);
+
+ return sqrt(sum);
+ }
+};
+typedef PairwiseDistances<PointContainer, L2Distance> PairDistances;
+typedef PairDistances::DistanceType DistanceType;
+
+typedef PairDistances::IndexType Vertex;
+typedef Simplex<Vertex> Smplx;
+typedef std::vector<Smplx> SimplexVector;
+typedef std::list<Smplx> SimplexList;
+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 ImageZigzagPersistence<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(DistanceType dist = DistanceType(), Dimension dim = Dimension()):
+ distance(dist), dimension(dim) {}
+ DistanceType distance;
+ Dimension dimension;
+};
+
+// Forward declarations of auxilliary functions
+void report_death(std::ofstream& out, Death d, DistanceType epsilon, Dimension skeleton_dimension);
+void make_boundary(const Smplx& s, Complex& c, const Zigzag& zz, Boundary& b);
+void show_image_betti(Zigzag& zz, Dimension skeleton);
+std::ostream& operator<<(std::ostream& out, const BirthInfo& bi);
+void process_command_line_options(int argc,
+ char* argv[],
+ unsigned& ambient_dimension,
+ unsigned& skeleton_dimension,
+ float& from_multiplier,
+ float& to_multiplier,
+ std::string& infilename,
+ std::string& outfilename);
+
+int main(int argc, char* argv[])
+{
+#ifdef LOGGING
+ rlog::RLogInit(argc, argv);
+
+ stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+#endif
+
+ Timer total, remove, add, ec, vc;
+ total.start();
+
+#if 0
+ SetFrequency(cOperations, 25000);
+ SetTrigger(cOperations, cComplexSize);
+#endif
+
+ unsigned ambient_dimension;
+ unsigned skeleton_dimension;
+ float from_multiplier, to_multiplier;
+ std::string infilename, outfilename;
+ process_command_line_options(argc, argv, ambient_dimension, skeleton_dimension, from_multiplier, to_multiplier, infilename, outfilename);
+
+ // Read in points
+ std::ifstream in(infilename.c_str());
+ PointContainer points;
+ while(in)
+ {
+ points.push_back(Point());
+ for (unsigned i = 0; i < ambient_dimension; ++i)
+ {
+ DistanceType x;
+ in >> x;
+ points.back().push_back(x);
+ }
+ }
+
+ // Create output file
+ std::ofstream out(outfilename.c_str());
+
+ // Create pairwise distances
+ PairDistances distances(points);
+
+ // Order vertices and epsilons (in maxmin fashion)
+ VertexVector vertices;
+ EpsilonVector epsilons;
+ EdgeVector edges;
+
+ {
+ EpsilonVector dist(distances.size(), Infinity);
+
+ vertices.push_back(distances.begin());
+ //epsilons.push_back(Infinity);
+ while (vertices.size() < distances.size())
+ {
+ for (Vertex v = distances.begin(); v != distances.end(); ++v)
+ dist[v] = std::min(dist[v], distances(v, vertices.back()));
+ EpsilonVector::const_iterator max = std::max_element(dist.begin(), dist.end());
+ vertices.push_back(max - dist.begin());
+ epsilons.push_back(*max);
+ }
+ epsilons.push_back(0);
+ }
+
+ rInfo("Point and epsilon ordering:");
+ for (unsigned i = 0; i < vertices.size(); ++i)
+ rInfo(" %4d: %4d - %f", i, vertices[i], epsilons[i]);
+
+ // Generate and sort all the edges
+ for (unsigned i = 0; i != vertices.size(); ++i)
+ for (unsigned j = i+1; j != vertices.size(); ++j)
+ {
+ Vertex u = vertices[i];
+ Vertex v = vertices[j];
+ if (distances(u,v) <= to_multiplier*epsilons[j-1])
+ edges.push_back(std::make_pair(u,v));
+ }
+ std::sort(edges.begin(), edges.end(), RipsGenerator::ComparePair(distances));
+ rInfo("Total participating edges: %d", edges.size());
+ for (EdgeVector::const_iterator cur = edges.begin(); cur != edges.end(); ++cur)
+ rDebug(" (%d, %d) %f", cur->first, cur->second, distances(cur->first, cur->second));
+
+ // Construct zigzag
+ Complex complex;
+ Zigzag zz;
+ RipsGenerator rips(distances);
+ SimplexEvaluator size(distances);
+
+ // Insert vertices
+ for (unsigned i = 0; i != vertices.size(); ++i)
+ {
+ // Add a vertex
+ Smplx sv; sv.add(vertices[i]);
+ rDebug("Adding %s", tostring(sv).c_str());
+ add.start();
+ complex.insert(std::make_pair(sv,
+ zz.add(Boundary(),
+ true, // vertex is always in the subcomplex
+ BirthInfo(0, 0)).first));
+ add.stop();
+ //rDebug("Newly born cycle order: %d", complex[sv]->low->order);
+ CountNum(cComplexSize, 0);
+ Count(cComplexSize);
+ Count(cOperations);
+ }
+
+ rInfo("Commencing computation");
+ boost::progress_display show_progress(vertices.size());
+ unsigned sce = 0, // index of the current one past last edge in the subcomplex
+ ce = 0; // index of the current one past last edge in the complex
+ for (unsigned stage = 0; stage != vertices.size() - 1; ++stage)
+ {
+ unsigned i = vertices.size() - 1 - stage;
+ rInfo("Current stage %d: %d %f: %f -> %f", stage,
+ vertices[i], epsilons[i-1],
+ from_multiplier*epsilons[i-1],
+ to_multiplier*epsilons[i-1]);
+
+ /* Increase epsilon */
+ // Record the cofaces of all the simplices that need to be removed and reinserted
+ SimplexSet cofaces;
+ rDebug(" Cofaces size: %d", cofaces.size());
+ while(sce < ce)
+ {
+ Vertex u,v;
+ boost::tie(u,v) = edges[sce];
+ if (distances(u,v) <= from_multiplier*epsilons[i-1])
+ ++sce;
+ else
+ break;
+
+ // Skip an edge if any one of its vertices has been removed from the complex
+ bool skip_edge = false;
+ for (unsigned j = i+1; j != vertices.size(); ++j)
+ if (u == vertices[j] || v == vertices[j])
+ {
+ // Debug only: eventually remove
+ rDebug(" Skipping edge (%d, %d)", u, v);
+ Smplx s; s.add(u); s.add(v);
+ AssertMsg(complex.find(s) == complex.end(), "Simplex should not be in the complex.");
+ skip_edge = true;
+ break;
+ }
+ if (skip_edge) continue;
+ rDebug(" Generating cofaces for (%d, %d)", u, v);
+
+ ec.start();
+ rips.edge_cofaces(u, v,
+ skeleton_dimension,
+ to_multiplier*epsilons[i],
+ make_insert_functor(cofaces),
+ vertices.begin(),
+ vertices.begin() + i + 1);
+ ec.stop();
+ }
+ rDebug(" Recorded cofaces to remove");
+ rDebug(" Cofaces size: %d", cofaces.size());
+ // Remove all the cofaces
+ for (SimplexSet::const_reverse_iterator cur = cofaces.rbegin(); cur != cofaces.rend(); ++cur)
+ {
+ rDebug(" Removing %s", tostring(*cur).c_str());
+ Complex::iterator si = complex.find(*cur);
+ remove.start();
+ AssertMsg(!si->second->subcomplex, "We should not remove simplices already in the subcomplex when we increase epsilon");
+ Death d = zz.remove(si->second,
+ BirthInfo(epsilons[i-1], cur->dimension() - 1));
+ remove.stop();
+ complex.erase(si);
+ CountNumBy(cComplexSize, cur->dimension(), -1);
+ CountBy(cComplexSize, -1);
+ Count(cOperations);
+ AssertMsg(zz.check_consistency(), "Zigzag representation must be consistent after removing a simplex");
+ report_death(out, d, epsilons[i-1], skeleton_dimension);
+ }
+ rDebug(" Removed cofaces");
+
+ // Add anything else that needs to be inserted into the complex
+ while (ce < edges.size())
+ {
+ Vertex u,v;
+ boost::tie(u,v) = edges[ce];
+ if (distances(u,v) <= to_multiplier*epsilons[i-1])
+ ++ce;
+ else
+ break;
+ rDebug(" Recording cofaces of edges[%d]=(%d, %d) with size=%f", (ce-1), u, v, distances(u,v));
+ ec.start();
+ rips.edge_cofaces(u, v,
+ skeleton_dimension,
+ to_multiplier*epsilons[i-1],
+ make_insert_functor(cofaces),
+ vertices.begin(),
+ vertices.begin() + i + 1);
+ ec.stop();
+ }
+ rDebug(" Recorded new cofaces to add");
+
+ // Progress sce
+ while (sce < ce)
+ {
+ Vertex u,v;
+ boost::tie(u,v) = edges[sce];
+ rDebug(" Progressing sce=%d over (%d, %d) %f", sce, u, v, distances(u,v));
+ if (distances(u,v) <= from_multiplier*epsilons[i-1])
+ ++sce;
+ else
+ break;
+ }
+ rDebug(" Moved subcomplex index forward");
+
+ // Insert all the cofaces
+ rDebug(" Cofaces size: %d", cofaces.size());
+ for (SimplexSet::const_iterator cur = cofaces.begin(); cur != cofaces.end(); ++cur)
+ {
+ Index idx; Death d; Boundary b;
+ rDebug(" Adding %s, its size %f", tostring(*cur).c_str(), size(*cur));
+ make_boundary(*cur, complex, zz, b);
+ add.start();
+ boost::tie(idx, d) = zz.add(b,
+ size(*cur) <= from_multiplier*epsilons[i-1],
+ BirthInfo(epsilons[i-1], cur->dimension()));
+ add.stop();
+ //if (!d) rDebug("Newly born cycle order: %d", complex[*cur]->low->order);
+ CountNum(cComplexSize, cur->dimension());
+ Count(cComplexSize);
+ Count(cOperations);
+ AssertMsg(zz.check_consistency(), "Zigzag representation must be consistent after removing a simplex");
+ complex.insert(std::make_pair(*cur, idx));
+ report_death(out, d, epsilons[i-1], skeleton_dimension);
+ }
+ rInfo("Increased epsilon; complex size: %d", complex.size());
+ show_image_betti(zz, skeleton_dimension);
+ report_memory();
+
+ /* Remove the vertex */
+ cofaces.clear();
+ rDebug(" Cofaces size: %d", cofaces.size());
+ vc.start();
+ rips.vertex_cofaces(vertices[i],
+ skeleton_dimension,
+ to_multiplier*epsilons[i-1],
+ make_insert_functor(cofaces),
+ vertices.begin(),
+ vertices.begin() + i + 1);
+ vc.stop();
+ rDebug(" Computed cofaces of the vertex, their number: %d", cofaces.size());
+ for (SimplexSet::const_reverse_iterator cur = cofaces.rbegin(); cur != cofaces.rend(); ++cur)
+ {
+ rDebug(" Removing: %s", tostring(*cur).c_str());
+ Complex::iterator si = complex.find(*cur);
+ remove.start();
+ Death d = zz.remove(si->second,
+ BirthInfo(epsilons[i-1], cur->dimension() - 1));
+ remove.stop();
+ complex.erase(si);
+ CountNumBy(cComplexSize, cur->dimension(), -1);
+ CountBy(cComplexSize, -1);
+ Count(cOperations);
+ AssertMsg(zz.check_consistency(), "Zigzag representation must be consistent after removing a simplex");
+ report_death(out, d, epsilons[i-1], skeleton_dimension);
+ }
+ rInfo("Removed vertex; complex size: %d", complex.size());
+ for (Complex::const_iterator cur = complex.begin(); cur != complex.end(); ++cur)
+ rDebug(" %s", tostring(cur->first).c_str());
+ show_image_betti(zz, skeleton_dimension);
+ report_memory();
+
+ ++show_progress;
+ }
+
+ // Remove the last vertex
+ AssertMsg(complex.size() == 1, "Only one vertex must remain");
+ remove.start();
+ Death d = zz.remove(complex.begin()->second, BirthInfo(epsilons[0], -1));
+ remove.stop();
+ complex.erase(complex.begin());
+ if (!d) AssertMsg(false, "The vertex must have died");
+ report_death(out, d, epsilons[0], skeleton_dimension);
+ CountNumBy(cComplexSize, 0, -1);
+ CountBy(cComplexSize, -1);
+ Count(cOperations);
+ rInfo("Removed vertex; complex size: %d", complex.size());
+ ++show_progress;
+
+ total.stop();
+
+ remove.check("Remove timer ");
+ add.check ("Add timer ");
+ ec.check ("Edge coface timer ");
+ vc.check ("Vertex coface timer ");
+ total.check ("Total timer ");
+}
+
+
+
+void report_death(std::ofstream& out, Death d, DistanceType epsilon, Dimension skeleton_dimension)
+{
+ if (d && ((d->distance - epsilon) != 0) && (d->dimension < skeleton_dimension))
+ out << d->dimension << " " << d->distance << " " << epsilon << std::endl;
+}
+
+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 (inL=%d)", c[*cur]->order, b.back()->subcomplex);
+ }
+}
+
+bool face_leaving_subcomplex(Complex::reverse_iterator si, const SimplexEvaluator& size, DistanceType after, DistanceType before)
+{
+ const Smplx& s = si->first;
+ for (Smplx::VertexContainer::const_iterator v1 = s.vertices().begin(); v1 != s.vertices().end(); ++v1)
+ for (Smplx::VertexContainer::const_iterator v2 = boost::next(v1); v2 != s.vertices().end(); ++v2)
+ {
+ Smplx e; e.add(*v1); e.add(*v2);
+ if (size(e) > after && size(e) <= before)
+ return true;
+ }
+
+ return false;
+}
+
+void show_image_betti(Zigzag& zz, Dimension skeleton)
+{
+ for (Zigzag::ZIndex cur = zz.image_begin(); cur != zz.image_end(); ++cur)
+ if (cur->low == zz.boundary_end() && cur->birth.dimension < skeleton)
+ rInfo("Class in the image of dimension: %d", cur->birth.dimension);
+}
+
+
+std::ostream& operator<<(std::ostream& out, const BirthInfo& bi)
+{ return (out << bi.distance); }
+
+void process_command_line_options(int argc,
+ char* argv[],
+ unsigned& ambient_dimension,
+ unsigned& skeleton_dimension,
+ float& from_multiplier,
+ float& to_multiplier,
+ std::string& infilename,
+ std::string& outfilename)
+{
+ namespace po = boost::program_options;
+
+ po::options_description hidden("Hidden options");
+ hidden.add_options()
+ ("input-file", po::value<std::string>(&infilename), "Point set whose Rips zigzag we want to compute")
+ ("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")
+ ("from,f", po::value<float>(&from_multiplier)->default_value(4), "From multiplier for the epsilon (distance to next maxmin point) when computing the Rips complex")
+ ("to,t", po::value<float>(&to_multiplier)->default_value(16), "To multiplier for the epsilon (distance to next maxmin point) 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("input-file", 1);
+ pos.add("output-file", 2);
+
+ 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("input-file") || !vm.count("output-file"))
+ {
+ std::cout << "Usage: " << argv[0] << " [options] input-file output-file" << std::endl;
+ std::cout << visible << std::endl;
+ std::abort();
+ }
+}
--- a/examples/rips/rips-zigzag.cpp Sat Jan 31 23:02:22 2009 -0800
+++ b/examples/rips/rips-zigzag.cpp Tue Feb 03 11:31:10 2009 -0800
@@ -1,18 +1,21 @@
#include <topology/rips.h>
-#include <topology/image-zigzag-persistence.h>
+#include <topology/zigzag-persistence.h>
#include <utilities/types.h>
#include <utilities/containers.h>
#include <utilities/log.h>
-#include <utilities/memory.h>
+#include <utilities/memory.h> // for report_memory()
+#include <utilities/timer.h>
#include <map>
#include <cmath>
#include <fstream>
#include <stack>
+#include <cstdlib>
#include <boost/tuple/tuple.hpp>
#include <boost/program_options.hpp>
+#include <boost/progress.hpp>
#ifdef COUNTERS
@@ -27,7 +30,7 @@
{
result_type operator()(const Point& p1, const Point& p2) const
{
- AssertMsg(p1.size() == p2.size(), "Points must be in the same dimension (in L2Distance");
+ AssertMsg(p1.size() == p2.size(), "Points must be in the same dimension (in L2Distance)");
result_type sum = 0;
for (size_t i = 0; i < p1.size(); ++i)
sum += (p1[i] - p2[i])*(p1[i] - p2[i]);
@@ -40,10 +43,26 @@
typedef PairDistances::IndexType Vertex;
typedef Simplex<Vertex> Smplx;
+typedef std::vector<Smplx> SimplexVector;
+typedef std::list<Smplx> SimplexList;
+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(DistanceType dist = DistanceType(), Dimension dim = Dimension()):
@@ -51,111 +70,40 @@
DistanceType distance;
Dimension dimension;
};
-typedef ImageZigzagPersistence<BirthInfo> Zigzag;
-typedef Zigzag::SimplexIndex Index;
-typedef Zigzag::Death Death;
-typedef std::map<Smplx, Index,
- Smplx::VertexDimensionComparison> Complex;
-typedef Zigzag::ZColumn Boundary;
-
-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 (inL=%d)", c[*cur]->order, b.back()->subcomplex);
- }
-}
-
-bool face_leaving_subcomplex(Complex::reverse_iterator si, const SimplexEvaluator& size, DistanceType after, DistanceType before)
-{
- const Smplx& s = si->first;
- for (Smplx::VertexContainer::const_iterator v1 = s.vertices().begin(); v1 != s.vertices().end(); ++v1)
- for (Smplx::VertexContainer::const_iterator v2 = boost::next(v1); v2 != s.vertices().end(); ++v2)
- {
- Smplx e; e.add(*v1); e.add(*v2);
- if (size(e) > after && size(e) <= before)
- return true;
- }
-
- return false;
-}
-
-void show_image_betti(Zigzag& zz, Dimension skeleton)
-{
- for (Zigzag::ZIndex cur = zz.image_begin(); cur != zz.image_end(); ++cur)
- if (cur->low == zz.boundary_end() && cur->birth.dimension< skeleton)
- rInfo("Class in the image of dimension: %d", cur->birth.dimension);
-}
-
-
-std::ostream& operator<<(std::ostream& out, const BirthInfo& bi)
-{ return (out << bi.distance); }
-
-namespace po = boost::program_options;
-
+// Forward declarations of auxilliary functions
+void report_death(std::ofstream& out, Death d, DistanceType epsilon, Dimension skeleton_dimension);
+void make_boundary(const Smplx& s, Complex& c, const Zigzag& zz, Boundary& b);
+std::ostream& operator<<(std::ostream& out, const BirthInfo& bi);
+void process_command_line_options(int argc,
+ char* argv[],
+ unsigned& ambient_dimension,
+ unsigned& skeleton_dimension,
+ float& multiplier,
+ std::string& infilename,
+ std::string& outfilename);
int main(int argc, char* argv[])
{
#ifdef LOGGING
rlog::RLogInit(argc, argv);
- stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
+ stderrLog.subscribeTo( RLOG_CHANNEL("error") );
#endif
-
+
+ Timer total, remove, add, ec, vc;
+ total.start();
+
+#if 0
SetFrequency(cOperations, 25000);
SetTrigger(cOperations, cComplexSize);
+#endif
unsigned ambient_dimension;
unsigned skeleton_dimension;
- float from_multiplier, to_multiplier;
+ float multiplier;
std::string infilename, outfilename;
-
- po::options_description hidden("Hidden options");
- hidden.add_options()
- ("input-file", po::value<std::string>(&infilename), "Point set whose Rips zigzag we want to compute")
- ("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")
- ("from,f", po::value<float>(&from_multiplier)->default_value(4), "From multiplier for the epsilon (distance to next maxmin point) when computing the Rips complex")
- ("to,t", po::value<float>(&to_multiplier)->default_value(16), "To multiplier for the epsilon (distance to next maxmin point) 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("input-file", 1);
- pos.add("output-file", 2);
-
- 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)
- stdoutLog.subscribeTo( RLOG_CHANNEL(cur->c_str()) );
- /* Interesting channels
- * "info", "debug", "topology/persistence"
- */
-#endif
-
- if (vm.count("help") || !vm.count("input-file"))
- {
- std::cout << "Usage: " << argv[0] << " [options] input-file output-file" << std::endl;
- std::cout << visible << std::endl;
- return 1;
- }
+ process_command_line_options(argc, argv, ambient_dimension, skeleton_dimension, multiplier, infilename, outfilename);
// Read in points
std::ifstream in(infilename.c_str());
@@ -177,12 +125,10 @@
// Create pairwise distances
PairDistances distances(points);
- // Order vertices and epsilons
- typedef std::vector<Vertex> VertexVector;
- typedef std::vector<DistanceType> EpsilonVector;
-
+ // Order vertices and epsilons (in maxmin fashion)
VertexVector vertices;
EpsilonVector epsilons;
+ EdgeVector edges;
{
EpsilonVector dist(distances.size(), Infinity);
@@ -199,147 +145,251 @@
}
epsilons.push_back(0);
}
-
+
rInfo("Point and epsilon ordering:");
for (unsigned i = 0; i < vertices.size(); ++i)
- rInfo(" %d %f", vertices[i], epsilons[i]);
+ rInfo(" %4d: %4d - %f", i, vertices[i], epsilons[i]);
+ // Generate and sort all the edges
+ for (unsigned i = 0; i != vertices.size(); ++i)
+ for (unsigned j = i+1; j != vertices.size(); ++j)
+ {
+ Vertex u = vertices[i];
+ Vertex v = vertices[j];
+ if (distances(u,v) <= multiplier*epsilons[j-1])
+ edges.push_back(std::make_pair(u,v));
+ }
+ std::sort(edges.begin(), edges.end(), RipsGenerator::ComparePair(distances));
+ rInfo("Total participating edges: %d", edges.size());
+ for (EdgeVector::const_iterator cur = edges.begin(); cur != edges.end(); ++cur)
+ rDebug(" (%d, %d) %f", cur->first, cur->second, distances(cur->first, cur->second));
// Construct zigzag
Complex complex;
Zigzag zz;
- RipsGenerator aux(distances);
+ RipsGenerator rips(distances);
SimplexEvaluator size(distances);
- // TODO: it probably makes sense to do things in reverse.
- // I.e., we should start from the smallest epsilon, and grow, rather than
- // starting from the largest epsilon and shrinking since the interesting
- // part of the computation is that with small epsilon.
- rInfo("Commencing computation");
+ // Insert vertices
for (unsigned i = 0; i != vertices.size(); ++i)
{
- rInfo("Current stage %d: %d %f", i, vertices[i], epsilons[i]);
-
- // Add a point
+ // Add a vertex
Smplx sv; sv.add(vertices[i]);
- rDebug("Added %s", tostring(sv).c_str());
+ rDebug("Adding %s", tostring(sv).c_str());
+ add.start();
complex.insert(std::make_pair(sv,
zz.add(Boundary(),
- true, // vertex is always in the subcomplex
- BirthInfo(epsilons[i], 0)).first));
+ BirthInfo(0, 0)).first));
+ add.stop();
+ //rDebug("Newly born cycle order: %d", complex[sv]->low->order);
CountNum(cComplexSize, 0);
Count(cComplexSize);
Count(cOperations);
- if (!zz.check_consistency())
- {
- //zz.show_all();
- rError("Zigzag representation must be consistent after adding a vertex");
- }
- for (Complex::iterator si = complex.begin(); si != complex.end(); ++si)
+ }
+
+ rInfo("Commencing computation");
+ boost::progress_display show_progress(vertices.size());
+ unsigned ce = 0; // index of the current one past last edge in the complex
+ for (unsigned stage = 0; stage != vertices.size() - 1; ++stage)
+ {
+ unsigned i = vertices.size() - 1 - stage;
+ rInfo("Current stage %d: %d %f: %f", stage,
+ vertices[i], epsilons[i-1],
+ multiplier*epsilons[i-1]);
+
+ /* Increase epsilon */
+ // Record the cofaces of all the simplices that need to be removed and reinserted
+ SimplexSet cofaces;
+ rDebug(" Cofaces size: %d", cofaces.size());
+
+ // Add anything else that needs to be inserted into the complex
+ while (ce < edges.size())
{
- if (si->first.contains(sv)) continue;
- if (si->first.dimension() + 1 > skeleton_dimension) continue;
-
- rDebug(" Distance between %s and %s: %f",
- tostring(si->first).c_str(),
- tostring(sv).c_str(),
- aux.distance(si->first, sv));
- if (aux.distance(si->first, sv) <= to_multiplier*epsilons[i-1])
- {
- Boundary b;
- Smplx s(si->first); s.join(sv);
-
- //zz.show_all();
- rDebug("Adding %s", tostring(s).c_str());
- make_boundary(s, complex, zz, b);
- rDebug("Made boundary, %d", b.size());
- Index idx; Death d;
- boost::tie(idx, d) = zz.add(b,
- (size(s) <= from_multiplier*epsilons[i-1]),
- BirthInfo(epsilons[i-1], s.dimension()));
- if (!zz.check_consistency())
- {
- //zz.show_all();
- rError("Zigzag representation must be consistent after adding a simplex");
- }
- complex.insert(std::make_pair(s, idx));
- CountNum(cComplexSize, s.dimension());
- Count(cComplexSize);
- Count(cOperations);
-
- // Death
- if (d && ((d->distance - epsilons[i-1]) != 0) && (d->dimension < skeleton_dimension))
- out << d->dimension << " " << d->distance << " " << epsilons[i-1] << std::endl;
- }
+ Vertex u,v;
+ boost::tie(u,v) = edges[ce];
+ if (distances(u,v) <= multiplier*epsilons[i-1])
+ ++ce;
+ else
+ break;
+ rDebug(" Recording cofaces of edges[%d]=(%d, %d) with size=%f", (ce-1), u, v, distances(u,v));
+ ec.start();
+ rips.edge_cofaces(u, v,
+ skeleton_dimension,
+ multiplier*epsilons[i-1],
+ make_insert_functor(cofaces),
+ vertices.begin(),
+ vertices.begin() + i + 1);
+ ec.stop();
}
- rDebug("Complex after addition:");
- for (Complex::const_iterator si = complex.begin(); si != complex.end(); ++si)
- rDebug(" %s", tostring(si->first).c_str());
+ rDebug(" Recorded new cofaces to add");
- rInfo("Inserted point; complex size: %d", complex.size());
- show_image_betti(zz, skeleton_dimension);
- report_memory();
-
- if (i == 0) continue; // want to skip the removal from the image check (involving epsilons[i-1]),
- // and in any case, there is only one vertex at this point
- rDebug("Removing simplices");
- // Shrink epsilon
+ // Insert all the cofaces
+ rDebug(" Cofaces size: %d", cofaces.size());
+ for (SimplexSet::const_iterator cur = cofaces.begin(); cur != cofaces.end(); ++cur)
{
- std::stack<Complex::reverse_iterator> leaving_subcomplex;
- Complex::reverse_iterator si = complex.rbegin();
+ Index idx; Death d; Boundary b;
+ rDebug(" Adding %s, its size %f", tostring(*cur).c_str(), size(*cur));
+ make_boundary(*cur, complex, zz, b);
+ add.start();
+ boost::tie(idx, d) = zz.add(b,
+ BirthInfo(epsilons[i-1], cur->dimension()));
+ add.stop();
+ //if (!d) rDebug("Newly born cycle order: %d", complex[*cur]->low->order);
+ CountNum(cComplexSize, cur->dimension());
+ Count(cComplexSize);
+ Count(cOperations);
+ AssertMsg(zz.check_consistency(), "Zigzag representation must be consistent after removing a simplex");
+ complex.insert(std::make_pair(*cur, idx));
+ report_death(out, d, epsilons[i-1], skeleton_dimension);
+ }
+ rInfo("Increased epsilon; complex size: %d", complex.size());
+ report_memory();
+
+ /* Remove the vertex */
+ cofaces.clear();
+ rDebug(" Cofaces size: %d", cofaces.size());
+ vc.start();
+ rips.vertex_cofaces(vertices[i],
+ skeleton_dimension,
+ multiplier*epsilons[i-1],
+ make_insert_functor(cofaces),
+ vertices.begin(),
+ vertices.begin() + i + 1);
+ vc.stop();
+ rDebug(" Computed cofaces of the vertex, their number: %d", cofaces.size());
+ for (SimplexSet::const_reverse_iterator cur = cofaces.rbegin(); cur != cofaces.rend(); ++cur)
+ {
+ rDebug(" Removing: %s", tostring(*cur).c_str());
+ Complex::iterator si = complex.find(*cur);
+ remove.start();
+ Death d = zz.remove(si->second,
+ BirthInfo(epsilons[i-1], cur->dimension() - 1));
+ remove.stop();
+ complex.erase(si);
+ CountNumBy(cComplexSize, cur->dimension(), -1);
+ CountBy(cComplexSize, -1);
+ Count(cOperations);
+ AssertMsg(zz.check_consistency(), "Zigzag representation must be consistent after removing a simplex");
+ report_death(out, d, epsilons[i-1], skeleton_dimension);
+ }
+ rInfo("Removed vertex; complex size: %d", complex.size());
+ for (Complex::const_iterator cur = complex.begin(); cur != complex.end(); ++cur)
+ rDebug(" %s", tostring(cur->first).c_str());
+ report_memory();
+
+ ++show_progress;
+ }
+
+ // Remove the last vertex
+ AssertMsg(complex.size() == 1, "Only one vertex must remain");
+ remove.start();
+ Death d = zz.remove(complex.begin()->second, BirthInfo(epsilons[0], -1));
+ remove.stop();
+ complex.erase(complex.begin());
+ if (!d) AssertMsg(false, "The vertex must have died");
+ report_death(out, d, epsilons[0], skeleton_dimension);
+ CountNumBy(cComplexSize, 0, -1);
+ CountBy(cComplexSize, -1);
+ Count(cOperations);
+ rInfo("Removed vertex; complex size: %d", complex.size());
+ ++show_progress;
+
+ total.stop();
+
+ remove.check("Remove timer ");
+ add.check ("Add timer ");
+ ec.check ("Edge coface timer ");
+ vc.check ("Vertex coface timer ");
+ total.check ("Total timer ");
+}
+
+
- while(si != complex.rend())
- {
- rDebug(" Size of %s is %f", tostring(si->first).c_str(), size(si->first));
- if (size(si->first) > to_multiplier*epsilons[i])
- {
- //zz.show_all();
- rDebug(" Removing from complex: %s", tostring(si->first).c_str());
- Death d = zz.remove(si->second,
- BirthInfo(epsilons[i], si->first.dimension() - 1));
- AssertMsg(zz.check_consistency(), "Zigzag representation must be consistent after removing a simplex");
- if (d && ((d->distance - epsilons[i]) != 0) && (d->dimension < skeleton_dimension))
- out << d->dimension << " " << d->distance << " " << epsilons[i] << std::endl;
- CountNumBy(cComplexSize, si->first.dimension(), -1);
- complex.erase(boost::prior(si.base()));
- CountBy(cComplexSize, -1);
- Count(cOperations);
- } else if (face_leaving_subcomplex(si, size, from_multiplier*epsilons[i], from_multiplier*epsilons[i-1]))
- {
- // Remove from subcomplex (i.e., remove and reinsert as outside of the subcomplex)
- rDebug(" Removing from subcomplex: %s", tostring(si->first).c_str());
- Death d = zz.remove(si->second,
- BirthInfo(epsilons[i], si->first.dimension() - 1));
- Count(cOperations);
- AssertMsg(zz.check_consistency(), "Zigzag representation must be consistent after removing a simplex");
- if (d && ((d->distance - epsilons[i]) != 0) && (d->dimension < skeleton_dimension))
- out << d->dimension << " " << d->distance << " " << epsilons[i] << std::endl;
- leaving_subcomplex.push(si++);
- } else
- ++si;
- }
- while(!leaving_subcomplex.empty())
- {
- si = leaving_subcomplex.top(); // copying an iterator onto stack is probably Ok
- Boundary b;
- make_boundary(si->first, complex, zz, b);
- Index idx; Death d;
- boost::tie(idx, d) = zz.add(b,
- false, // now it is outside of the subcomplex
- BirthInfo(epsilons[i], si->first.dimension()));
- Count(cOperations);
- si->second = idx;
- if (d && ((d->distance - epsilons[i]) != 0) && (d->dimension < skeleton_dimension))
- out << d->dimension << " " << d->distance << " " << epsilons[i] << std::endl;
- leaving_subcomplex.pop();
- }
- }
- rDebug("Complex after removal:");
- for (Complex::const_iterator si = complex.begin(); si != complex.end(); ++si)
- rDebug(" %s", tostring(si->first).c_str());
-
- rInfo("Shrunk epsilon; complex size: %d", complex.size());
- show_image_betti(zz, skeleton_dimension);
- report_memory();
+void report_death(std::ofstream& out, Death d, DistanceType epsilon, Dimension skeleton_dimension)
+{
+ if (d && ((d->distance - epsilon) != 0) && (d->dimension < skeleton_dimension))
+ out << d->dimension << " " << d->distance << " " << epsilon << std::endl;
+}
+
+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);
}
}
+
+bool face_leaving_subcomplex(Complex::reverse_iterator si, const SimplexEvaluator& size, DistanceType after, DistanceType before)
+{
+ const Smplx& s = si->first;
+ for (Smplx::VertexContainer::const_iterator v1 = s.vertices().begin(); v1 != s.vertices().end(); ++v1)
+ for (Smplx::VertexContainer::const_iterator v2 = boost::next(v1); v2 != s.vertices().end(); ++v2)
+ {
+ Smplx e; e.add(*v1); e.add(*v2);
+ if (size(e) > after && size(e) <= before)
+ return true;
+ }
+
+ return false;
+}
+
+
+std::ostream& operator<<(std::ostream& out, const BirthInfo& bi)
+{ return (out << bi.distance); }
+
+void process_command_line_options(int argc,
+ char* argv[],
+ unsigned& ambient_dimension,
+ unsigned& skeleton_dimension,
+ float& multiplier,
+ std::string& infilename,
+ std::string& outfilename)
+{
+ namespace po = boost::program_options;
+
+ po::options_description hidden("Hidden options");
+ hidden.add_options()
+ ("input-file", po::value<std::string>(&infilename), "Point set whose Rips zigzag we want to compute")
+ ("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")
+ ("multiplier,m", po::value<float>(&multiplier)->default_value(6), "Multiplier for the epsilon (distance to next maxmin point) 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("input-file", 1);
+ pos.add("output-file", 2);
+
+ 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("input-file") || !vm.count("output-file"))
+ {
+ std::cout << "Usage: " << argv[0] << " [options] input-file output-file" << std::endl;
+ std::cout << visible << std::endl;
+ std::abort();
+ }
+}
--- a/include/topology/zigzag-persistence.h Sat Jan 31 23:02:22 2009 -0800
+++ b/include/topology/zigzag-persistence.h Tue Feb 03 11:31:10 2009 -0800
@@ -94,6 +94,9 @@
// Function: remove(s, birth)
Death remove(SimplexIndex s, const BirthID& birth = BirthID()) { ZigzagVisitor zzv; return remove<ZigzagVisitor>(s, birth, zzv); }
+
+ ZIndex begin() { return z_list.begin(); }
+ ZIndex end() { return z_list.end(); }
protected:
// Function: add(s)
--- a/include/topology/zigzag-persistence.hpp Sat Jan 31 23:02:22 2009 -0800
+++ b/include/topology/zigzag-persistence.hpp Tue Feb 03 11:31:10 2009 -0800
@@ -12,6 +12,11 @@
static rlog::RLogChannel* rlZigzagCheckConsistency= DEF_CHANNEL("topology/persistence/zigzag/check", rlog::Log_Debug);
#endif // LOGGING
+#ifdef COUNTERS
+static Counter* cZigzagAdd = GetCounter("zigzag/add");
+static Counter* cZigzagRemove = GetCounter("zigzag/remove");
+static Counter* cZigzagConsistency = GetCounter("zigzag/consistency");
+#endif // COUNTERS
template<class BID, class SD>
template<class Visitor>
@@ -19,6 +24,8 @@
ZigzagPersistence<BID,SD>::
add(ZColumn bdry, const BirthID& birth, Visitor& visitor)
{
+ Count(cZigzagAdd);
+
rLog(rlZigzagAdd, "Entered ZigzagPersistence::add()");
rLog(rlZigzagAdd, " Boundary: %s", bdry.tostring(out).c_str());
rLog(rlZigzagAdd, " Boundary size: %d", bdry.size());
@@ -118,6 +125,8 @@
ZigzagPersistence<BID,SD>::
remove(SimplexIndex s, const BirthID& birth, Visitor& visitor)
{
+ Count(cZigzagRemove);
+
rLog(rlZigzagRemove, "Entered ZigzagPersistence::remove(%d)", s->order);
AssertMsg(check_consistency(), "Must be consistent before removal");
@@ -172,8 +181,9 @@
boost::make_filter_iterator(std::bind2nd(NotEqualBIndex(), j), first_z->b_row.end(), first_z->b_row.end()),
j, &BNode::c_column, &SimplexNode::c_row);
add_chain(j, j, &BNode::c_column, &SimplexNode::c_row);
- // TODO: remove add_chains(first_z->b_row.rbegin(), first_z->b_row.rend(), j, &BNode::c_column, &SimplexNode::c_row);
- AssertMsg(check_consistency(s_list.end(), z_list.begin(), b_list.end()), "Must be consistent after subtracting C[j] in remove::birth");
+ // TODO: that's how it was done before, now it can be removed
+ // add_chains(first_z->b_row.rbegin(), first_z->b_row.rend(), j, &BNode::c_column, &SimplexNode::c_row);
+ //AssertMsg(check_consistency(s_list.end(), z_list.begin(), b_list.end()), "Must be consistent after subtracting C[j] in remove::birth");
// Subtract B[j] from every other column of B that has l
ZIndex l = j->b_column.back();
@@ -185,7 +195,8 @@
j, &BNode::b_column, &ZNode::b_row);
j->b_column.back()->low = b_list.end(); // redundant since l will be deleted (here for check_consistency only)
add_chain(j, j, &BNode::b_column, &ZNode::b_row);
- AssertMsg(check_consistency(s_list.end(), first_z, b_list.end()), "Must be consistent after subtracting B[j] in remove::birth");
+ // TODO: investigate why this works for ordinary zigzag, but fails for the image zigzag
+ //AssertMsg(check_consistency(s_list.end(), first_z, b_list.end()), "Must be consistent after subtracting B[j] in remove::birth");
// Drop j, l, and s
@@ -406,6 +417,9 @@
check_consistency(SimplexIndex s_skip, ZIndex z_skip, BIndex b_skip)
{
#ifdef ZIGZAG_CONSISTENCY
+ #warning "Checking consistency in ZigzagPersistence"
+
+ Count(cZigzagConsistency);
for (SimplexIndex cur = s_list.begin(); cur != s_list.end(); ++cur)
{
if (cur == s_skip) continue;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/utilities/timer.h Tue Feb 03 11:31:10 2009 -0800
@@ -0,0 +1,69 @@
+#ifndef __TIMER_H__
+#define __TIMER_H__
+
+// Adapted with minor changes from http://oldmill.uchicago.edu/~wilder/Code/timer/
+
+#include <ctime>
+#include <iostream>
+#include <iomanip>
+
+class Timer
+{
+ public:
+ Timer(): start_time(0), acc_time(0) {}
+
+ void start();
+ void stop();
+ void check(const char* msg = 0) const;
+
+ private:
+ clock_t start_clock;
+ time_t start_time;
+ double acc_time;
+
+ double elapsed_time() const;
+};
+
+// Return the total time that the timer has been in the "running" state since
+// it was last "started". For "short" time periods (less than an hour), the
+// actual cpu time used is reported instead of the elapsed time.
+inline double
+Timer::
+elapsed_time() const
+{
+ time_t acc_sec = time(0) - start_time;
+ if (acc_sec < 3600)
+ return (clock() - start_clock) / (1.0 * CLOCKS_PER_SEC);
+ else
+ return (1.0 * acc_sec);
+}
+
+inline void
+Timer::
+start()
+{
+ start_clock = clock();
+ start_time = time(0);
+}
+
+inline void
+Timer::
+stop()
+{
+ acc_time += elapsed_time();
+}
+
+// Print out an optional message followed by the current timer timing.
+inline void
+Timer::
+check(const char* msg) const
+{
+ // Print an optional message, something like "Checking timer t";
+ if (msg) std::cout << msg << " : ";
+
+ std::cout << "Elapsed time [" << std::setiosflags(std::ios::fixed)
+ << std::setprecision(2) << acc_time << "] seconds\n";
+}
+
+#endif // __TIMER_H__
+