--- a/examples/CMakeLists.txt Fri Apr 26 22:31:58 2013 -0700
+++ b/examples/CMakeLists.txt Mon Jun 10 16:58:01 2013 +0200
@@ -1,14 +1,15 @@
add_executable (bottleneck-distance bottleneck-distance.cpp)
target_link_libraries (bottleneck-distance ${libraries} ${Boost_PROGRAM_OPTIONS_LIBRARY})
-add_subdirectory (alphashapes)
-#add_subdirectory (ar-vineyard)
-add_subdirectory (cech-complex)
-add_subdirectory (consistency)
-add_subdirectory (cohomology)
-add_subdirectory (fitness)
-add_subdirectory (pl-functions)
-add_subdirectory (triangle)
-add_subdirectory (poincare)
-add_subdirectory (rips)
-add_subdirectory (filtration)
+add_subdirectory (alphashapes)
+#add_subdirectory (ar-vineyard)
+add_subdirectory (cech-complex)
+add_subdirectory (consistency)
+add_subdirectory (cohomology)
+add_subdirectory (fitness)
+add_subdirectory (pl-functions)
+add_subdirectory (triangle)
+add_subdirectory (poincare)
+add_subdirectory (rips)
+add_subdirectory (filtration)
+add_subdirectory (homology-zigzags)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/homology-zigzags/CMakeLists.txt Mon Jun 10 16:58:01 2013 +0200
@@ -0,0 +1,11 @@
+set (targets
+ rips-pairwise
+ M-ZZ
+ dM-ZZ
+ iR-ZZ
+ oR-ZZ)
+
+foreach (t ${targets})
+ add_executable (${t} ${t}.cpp)
+ target_link_libraries (${t} ${libraries} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+endforeach (t ${targets})
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/homology-zigzags/M-ZZ.cpp Mon Jun 10 16:58:01 2013 +0200
@@ -0,0 +1,430 @@
+/***************************************************************************
+
+ M-ZZ: Morozov zigzag implementation
+ Copyright (C) 2009-2012 Dmitriy Morozov
+
+ Adapted from its original version ("rips-zigzag.cpp" in
+ the Dionysus library) by Steve Oudot (2012).
+
+ Changelog (2012-11-22):
+
+ - the barcode is now output on a log scale and in a format that is
+ compatible with the Matlab layer of the PLEX 2.5 library,
+
+ - the barcode representation is now by closed intervals instead of
+ half-open intervals.
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+***************************************************************************/
+
+#include <topology/rips.h>
+#include <topology/zigzag-persistence.h>
+#include <utilities/types.h>
+#include <utilities/containers.h>
+
+#include <geometry/l2distance.h> // Point, PointContainer, L2DistanceType, read_points
+#include <geometry/distances.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 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 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()):
+ distance(dist), dimension(dim) {}
+ DistanceType distance;
+ Dimension dimension;
+};
+
+// Forward declarations of auxilliary functions
+// Note: min_value is used only for log-scale, so set it up to zero by default
+void write_intervals(std::ostream& out, std::list<std::pair<double,double> >* intervals, int skeleton_dimension, bool logscale, double min_value=0);
+void report_death(std::list<std::pair<double, double> >* intervals, Death d, DistanceType epsilon, DistanceType birthEpsilon, 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& skeleton_dimension,
+ float& multiplier,
+ bool& logscale,
+ 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 skeleton_dimension;
+ float multiplier;
+ std::string infilename, outfilename;
+ bool logscale = false;
+ process_command_line_options(argc, argv, skeleton_dimension, multiplier, logscale, infilename, outfilename);
+
+ // Read in points
+ PointContainer points;
+ read_points(infilename, points);
+
+ std::cout << "Number of points: " << points.size() << std::endl;
+
+ // Create output file
+ std::ofstream out(outfilename.c_str());
+
+ // Create pairwise distances
+ PairDistances distances(points);
+
+ // Create intervals DS
+ std::list<std::pair<double, double> > intervals [skeleton_dimension];
+ // for (int i=0; i<skeleton_dimension; i++)
+ // intervals[i] = new std::list<std::pair<double, double> > ();
+
+ // 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) <= multiplier*epsilons[j-1])
+ edges.push_back(std::make_pair(u,v));
+ }
+ std::sort(edges.begin(), edges.end(), RipsGenerator::ComparePair(distances));
+ std::cout << "Total participating edges: " << edges.size() << std::endl;
+ 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 (initial complex is just a disjoint union of 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(),
+ 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 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;
+ // std::cout << "Current stage " << stage << " "
+ // << vertices[i] << " " << epsilons[i-1] << ": "
+ // << multiplier*epsilons[i-1] << std::endl;
+
+ /* 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())
+ {
+ 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(" Recorded new cofaces to add");
+
+ // 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,
+ 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(intervals, d, epsilons[i-1], epsilons[i], 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 != (SimplexSet::const_reverse_iterator)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(intervals, d, epsilons[i-1], epsilons[i], 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(intervals, d, epsilons[0], 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 ");
+
+ std::cout << "Writing intervals...";
+ // Note (hack): use epsilons[vertices.size()-2]/2 as minimal value for the log-scale intervals to avoid intervals starting at -infinity and thus a scaling effect
+ write_intervals(out, intervals, skeleton_dimension, logscale, epsilons[vertices.size()-2]/2);
+ std::cout << " done." << std::endl;
+}
+
+
+const double LOG2 = std::log(2.0);
+double log2(double x) {
+ return std::log(x) / LOG2;
+}
+
+void write_intervals(std::ostream& out, std::list<std::pair<double,double> >* intervals, int skeleton_dimension, bool logscale, double min_value) {
+ out << "I = { ";
+ for (int d = 0; d<skeleton_dimension; d++) {
+ out << "[ ";
+ for (std::list<std::pair<double,double> >::iterator pit = intervals[d].begin(); pit != intervals[d].end(); pit++)
+ if (logscale)
+ out << "[" << log2(std::max(pit->first, min_value)) << ";" << log2(std::max(pit->second, min_value)) << "] ";
+ else
+ out << "[" << pit->first << ";" << pit->second << "] ";
+
+ // add dummy interval if intervals list is empty
+ if (intervals[d].empty())
+ out << "[0;0] ";
+ out << "] ,";
+ }
+ out << "} ";
+}
+
+void report_death(std::list<std::pair<double,double> >* intervals, Death d, DistanceType epsilon, DistanceType birthEpsilon, Dimension skeleton_dimension)
+{
+ // std::cerr << " d = " << d;
+ // if (d)
+ // std::cerr << " d->dim = " << d->dimension
+ // << " d->dist = " << d->distance;
+ // std::cerr << " epsilon = " << epsilon << std::endl;
+
+ if (d && ((d->distance - epsilon) != 0) && (d->dimension < skeleton_dimension))
+ intervals[d->dimension].push_back(std::pair<double,double>(d->distance, birthEpsilon));
+}
+
+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);
+ }
+}
+
+std::ostream& operator<<(std::ostream& out, const BirthInfo& bi)
+{ return (out << bi.distance); }
+
+void process_command_line_options(int argc,
+ char* argv[],
+ unsigned& skeleton_dimension,
+ float& multiplier,
+ bool& logscale,
+ 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")
+ ("dim,d", po::value<unsigned>(&skeleton_dimension)->default_value(3), "maximal dimension of the Rips complex")
+ ("rho,r", po::value<float>(&multiplier)->default_value(3), "multiplier rho used in the Rips parameter")
+ ( "logscale,l" , po::value( &logscale )->zero_tokens(), "outputs interval bounds on log_2 scale" ) ;
+#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();
+ }
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/homology-zigzags/README Mon Jun 10 16:58:01 2013 +0200
@@ -0,0 +1,43 @@
+********************************************************************************
+* Rips-ZZ: Rips Zigzags for Homology Inference *
+********************************************************************************
+
+This is an extension to the Rips package of the Dionysus library. It
+adapts the code of the Morozov zigzag (M-ZZ) and image Rips zigzag
+(iR-ZZ) to the context of the following paper:
+
+Title: Zigzag Zoology: Rips Zigzags for Homology Inference
+Authors: Steve Y. Oudot and Donald R. Sheehy
+
+It also provides two variants of the M-ZZ: the discretized Morozov
+zigzag (dM-ZZ) and the oscillating Rips zigzag (oR-ZZ).
+
+
+Source files:
+
+ - M-ZZ.cpp: implementation of the Morozov zigzag
+ - dM-ZZ.cpp: implementation of the discretized Morozov zigzag
+ - oR-ZZ.cpp: implementation of the oscillating Rips zigzag
+ - iR-ZZ.cpp: implementation of the image Rips zigzag
+ - rips-pairwise.cpp: computation of the standard Rips filtration
+ and its persistent homology.
+
+Execution:
+
+ - the list of arguments required by a program can be obtained by
+ running that program without arguments.
+
+ - every program takes in a point cloud file, in xyz... format (no
+ need to indicate the ambient dimension at the beginning of the
+ file. A sample point cloud is provided in the file
+ input/spiral_3d_10k.xyz
+
+ - every program also asks for the name of the output file in which
+ the barcode will be written. The output format is such that the
+ file can be executed in Matlab. This creates a cells structures
+ that can be read by the px_homologyplot function from the PLEX 2.5
+ library (a copy of which is provided). This function will plot the
+ barcode dimension per dimension.
+
+
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/homology-zigzags/dM-ZZ.cpp Mon Jun 10 16:58:01 2013 +0200
@@ -0,0 +1,449 @@
+/***************************************************************************
+
+ dM-ZZ: discretized Morozov zigzag implementation
+ Copyright (C) 2012 Steve Oudot
+
+ Adapted from the Morozov zigzag implementation provided in the
+ Rips package of the Dionysus library (see the file "M-ZZ.cpp").
+ The Dionysus library is (C) 2006-2012 Dmitriy Morozov.
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+***************************************************************************/
+
+#include <topology/rips.h>
+#include <topology/zigzag-persistence.h>
+#include <utilities/types.h>
+#include <utilities/containers.h>
+
+#include <geometry/l2distance.h> // Point, PointContainer, L2DistanceType, read_points
+#include <geometry/distances.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 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 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()):
+ distance(dist), dimension(dim) {}
+ DistanceType distance;
+ Dimension dimension;
+};
+
+// Forward declarations of auxilliary functions
+// Note: min_value is used only for log-scale, so set it up to zero by default
+void write_intervals(std::ostream& out, std::list<std::pair<double,double> >* intervals, int skeleton_dimension, bool logscale, double min_value=0);
+// void report_death(std::ostream& out, Death d, DistanceType epsilon, Dimension skeleton_dimension);
+void report_death(std::list<std::pair<double, double> >* intervals, Death d, DistanceType epsilon, DistanceType birthEpsilon, 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& skeleton_dimension,
+ float& multiplier,
+ float& discretization_factor,
+ bool& logscale,
+ 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 skeleton_dimension;
+ float multiplier;
+ float discretization_factor;
+ std::string infilename, outfilename;
+ bool logscale;
+ process_command_line_options(argc, argv, skeleton_dimension, multiplier, discretization_factor, logscale, infilename, outfilename);
+
+ // Read in points
+ PointContainer points;
+ read_points(infilename, points);
+
+ std::cout << "Number of points: " << points.size() << std::endl;
+
+ // Create output file
+ std::ofstream out(outfilename.c_str());
+
+ // Create pairwise distances
+ PairDistances distances(points);
+
+ // Create intervals DS
+ std::list<std::pair<double, double> > intervals [skeleton_dimension];
+ // for (int i=0; i<skeleton_dimension; i++)
+ // intervals[i] = new std::list<std::pair<double, double> > ();
+
+ // 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);
+ }
+
+ // Discretize sequence
+ unsigned anchor[vertices.size()];
+ anchor[0] = 0;
+ unsigned ref = 0;
+ // std::cout << "Anchors: 0 (" << epsilons[0] << ")";
+ for (unsigned i=0; i<vertices.size(); i++)
+ if (epsilons[i] < epsilons[anchor[ref]]*discretization_factor) {
+ anchor[i] = i;
+ ref = i;
+ // std::cout << " " << i << " (" << epsilons[i] << ")";
+ }
+ else
+ anchor[i] = ref;
+ // std::cout << std::endl;
+
+
+ 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) <= multiplier*epsilons[anchor[j-1]])
+ edges.push_back(std::make_pair(u,v));
+ }
+ std::sort(edges.begin(), edges.end(), RipsGenerator::ComparePair(distances));
+ std::cout << "Total participating edges: " << edges.size() << std::endl;
+ 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 (initial complex is just a disjoint union of 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(),
+ 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 ce = 0; // index of the current one past last edge in the complex
+ for (unsigned stage = vertices.size()-1; stage > 0; stage = anchor[stage-1])
+ {
+ // std::cout << "Current stage " << stage << ": "
+ // << anchor[stage-1] << " " << epsilons[anchor[stage-1]] << " "
+ // << multiplier*epsilons[anchor[stage-1]] << std::endl;
+
+ /* Increase epsilon */
+ // Record the cofaces of all the simplices that need to be removed and r`einserted
+ SimplexSet cofaces;
+ rDebug(" Cofaces size: %d", cofaces.size());
+
+ // 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) <= multiplier*epsilons[anchor[stage-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[anchor[stage-1]],
+ make_insert_functor(cofaces),
+ vertices.begin(),
+ vertices.begin() + stage + 1);
+ ec.stop();
+ }
+ rDebug(" Recorded new cofaces to add");
+
+ // 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,
+ BirthInfo(epsilons[anchor[stage-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(intervals, d, epsilons[anchor[stage-1]], epsilons[stage], skeleton_dimension);
+ }
+ rInfo("Increased epsilon; complex size: %d", complex.size());
+ report_memory();
+
+ /* Remove the vertices */
+ for (unsigned i = stage; i>anchor[stage-1]; i--) {
+ // std::cout << "removing vertex " << i << std::endl;
+ cofaces.clear();
+ rDebug(" Cofaces size: %d", cofaces.size());
+ vc.start();
+ rips.vertex_cofaces(vertices[i],
+ skeleton_dimension,
+ multiplier*epsilons[anchor[stage-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 != (SimplexSet::const_reverse_iterator)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[anchor[stage-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(intervals, d, epsilons[anchor[stage-1]], epsilons[stage], 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
+ // std::cout << "removing vertex 0" << std::endl;
+ if (complex.size() != 1) {
+ std::cerr << "Error: Only one vertex must remain" << std::endl;
+ abort();
+ }
+ 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(intervals, d, epsilons[0], 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 ");
+
+ std::cout << "Writing intervals...";
+ // Note (hack): use epsilons[vertices.size()-2]/2 as minimal value for the log-scale intervals to avoid intervals starting at -infinity and thus a scaling effect
+ write_intervals(out, intervals, skeleton_dimension, logscale, epsilons[vertices.size()-2]/2);
+ std::cout << " done." << std::endl;
+}
+
+
+const double LOG2 = std::log(2.0);
+double log2(double x) {
+ return std::log(x) / LOG2;
+}
+
+void write_intervals(std::ostream& out, std::list<std::pair<double,double> >* intervals, int skeleton_dimension, bool logscale, double min_value) {
+ out << "I = { ";
+ for (int d = 0; d<skeleton_dimension; d++) {
+ out << "[ ";
+ for (std::list<std::pair<double,double> >::iterator pit = intervals[d].begin(); pit != intervals[d].end(); pit++)
+ if (logscale)
+ out << "[" << log2(std::max(pit->first, min_value)) << ";" << log2(std::max(pit->second, min_value)) << "] ";
+ else
+ out << "[" << pit->first << ";" << pit->second << "] ";
+
+ // add dummy interval if intervals list is empty
+ if (intervals[d].empty())
+ out << "[0;0] ";
+ out << "] ,";
+ }
+ out << "} ";
+}
+
+void report_death(std::list<std::pair<double,double> >* intervals, Death d, DistanceType epsilon, DistanceType birthEpsilon, Dimension skeleton_dimension)
+{
+ // std::cerr << " d = " << d;
+ // if (d)
+ // std::cerr << " d->dim = " << d->dimension
+ // << " d->dist = " << d->distance;
+ // std::cerr << " epsilon = " << epsilon << std::endl;
+
+ if (d && ((d->distance - epsilon) != 0) && (d->dimension < skeleton_dimension))
+ intervals[d->dimension].push_back(std::pair<double,double>(d->distance, birthEpsilon));
+}
+
+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);
+ }
+}
+
+std::ostream& operator<<(std::ostream& out, const BirthInfo& bi)
+{ return (out << bi.distance); }
+
+void process_command_line_options(int argc,
+ char* argv[],
+ unsigned& skeleton_dimension,
+ float& multiplier,
+ float& discretization_factor,
+ bool& logscale,
+ 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")
+ ("dim,d", po::value<unsigned>(&skeleton_dimension)->default_value(3), "maximal dimension of the Rips complex")
+ ("rho,r", po::value<float>(&multiplier)->default_value(3), "multiplier rho used in the Rips parameter")
+ ("zeta,z", po::value<float>(&discretization_factor)->default_value(1), "geometric scale gap in the discretization (cst<=1)")
+ ( "logscale,l" , po::value( &logscale )->zero_tokens(), "outputs interval bounds on log_2 scale" ) ;
+#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();
+ }
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/homology-zigzags/iR-ZZ.cpp Mon Jun 10 16:58:01 2013 +0200
@@ -0,0 +1,523 @@
+/***************************************************************************
+
+ iR-ZZ: image Rips zigzag implementation
+ Copyright (C) 2009-2012 Dmitriy Morozov
+
+ Adapted from its original version ("rips-image-zigzag.cpp" in
+ the Dionysus library) by Steve Oudot (2012).
+
+ Changelog (2012-11-22):
+
+ - the barcode is now output on a log scale and in a format that is
+ compatible with the Matlab layer of the PLEX 2.5 library,
+
+ - the barcode representation is now by closed intervals instead of
+ half-open intervals.
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+***************************************************************************/
+
+#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 <geometry/l2distance.h> // for L2Distance and read_points()
+#include <geometry/distances.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 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
+// Note: min_value is used only for log-scale, so set it up to zero by default
+void write_intervals(std::ostream& out, std::list<std::pair<double,double> >* intervals, int skeleton_dimension, bool logscale, double min_value=0);
+void report_death(std::list<std::pair<double, double> >* intervals, Death d, DistanceType epsilon, DistanceType birthEpsilon, Dimension skeleton_dimension);
+// 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& skeleton_dimension,
+ float& from_multiplier,
+ float& to_multiplier,
+ bool& logscale,
+ 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 skeleton_dimension;
+ float from_multiplier, to_multiplier;
+ std::string infilename, outfilename;
+ bool logscale = false;
+ process_command_line_options(argc, argv, skeleton_dimension, from_multiplier, to_multiplier, logscale, infilename, outfilename);
+
+ // Read in points
+ PointContainer points;
+ read_points(infilename, points);
+
+ // Create output file
+ std::ofstream out(outfilename.c_str());
+
+ // Create pairwise distances
+ PairDistances distances(points);
+
+ // Create intervals DS
+ std::list<std::pair<double, double> > intervals [skeleton_dimension];
+ // for (int i=0; i<skeleton_dimension; i++)
+ // intervals[i] = new std::list<std::pair<double, double> > ();
+
+ // 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 != (SimplexSet::const_reverse_iterator)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(intervals, d, epsilons[i-1], epsilons[i], 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(intervals, d, epsilons[i-1], epsilons[i], 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 != (SimplexSet::const_reverse_iterator)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(intervals, d, epsilons[i-1], epsilons[i], 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(intervals, d, epsilons[0], 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 ");
+
+ std::cout << "Writing intervals...";
+ // Note (hack): use epsilons[vertices.size()-2]/2 as minimal value for the log-scale intervals to avoid intervals starting at -infinity and thus a scaling effect
+ write_intervals(out, intervals, skeleton_dimension, logscale, epsilons[vertices.size()-2]/2);
+ std::cout << " done." << std::endl;
+}
+
+
+void write_intervals(std::ostream& out, std::list<std::pair<double,double> >* intervals, int skeleton_dimension, bool logscale, double min_value) {
+ out << "I = { ";
+ for (int d = 0; d<skeleton_dimension; d++) {
+ out << "[ ";
+ for (std::list<std::pair<double,double> >::iterator pit = intervals[d].begin(); pit != intervals[d].end(); pit++)
+ if (logscale)
+ out << "[" << log2(std::max(pit->first, min_value)) << ";" << log2(std::max(pit->second, min_value)) << "] ";
+ else
+ out << "[" << pit->first << ";" << pit->second << "] ";
+
+ // add dummy interval if intervals list is empty
+ if (intervals[d].empty())
+ out << "[0;0] ";
+ out << "] ,";
+ }
+ out << "} ";
+}
+
+void report_death(std::list<std::pair<double,double> >* intervals, Death d, DistanceType epsilon, DistanceType birthEpsilon, Dimension skeleton_dimension)
+{
+ // std::cerr << " d = " << d;
+ // if (d)
+ // std::cerr << " d->dim = " << d->dimension
+ // << " d->dist = " << d->distance;
+ // std::cerr << " epsilon = " << epsilon << std::endl;
+
+ if (d && ((d->distance - epsilon) != 0) && (d->dimension < skeleton_dimension))
+ intervals[d->dimension].push_back(std::pair<double,double>(d->distance, birthEpsilon));
+}
+
+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& skeleton_dimension,
+ float& from_multiplier,
+ float& to_multiplier,
+ bool& logscale,
+ 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")
+ ("dim,d", po::value<unsigned>(&skeleton_dimension)->default_value(3), "maximal dimension of the Rips complex")
+ ("eta,e", po::value<float>(&from_multiplier)->default_value(3), "multiplier eta used in the Rips parameter")
+ ("rho,r", po::value<float>(&to_multiplier)->default_value(3), "multiplier rho used in the Rips parameter")
+ ( "logscale,l" , po::value<bool>(&logscale)->zero_tokens(), "outputs interval bounds on log_2 scale" ) ;
+#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();
+ }
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/homology-zigzags/oR-ZZ.cpp Mon Jun 10 16:58:01 2013 +0200
@@ -0,0 +1,500 @@
+/***************************************************************************
+
+ oR-ZZ: oscillating Rips zigzag implementation
+ Copyright (C) 2012 Steve Oudot
+
+ Adapted from the Morozov zigzag implementation provided in the
+ Rips package of the Dionysus library (see the file "M-ZZ.cpp").
+ The Dionysus library is (C) 2006-2012 Dmitriy Morozov.
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+***************************************************************************/
+
+#include <topology/rips.h>
+#include <topology/zigzag-persistence.h>
+#include <utilities/types.h>
+#include <utilities/containers.h>
+
+#include <geometry/l2distance.h> // Point, PointContainer, L2DistanceType, read_points
+#include <geometry/distances.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 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 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()):
+ distance(dist), dimension(dim) {}
+ DistanceType distance;
+ Dimension dimension;
+};
+
+// Forward declarations of auxilliary functions
+void write_intervals(std::ostream& out, std::list<std::pair<double,double> >* intervals, int skeleton_dimension, bool logscale, double min_value=0);
+// void report_death(std::ostream& out, Death d, DistanceType epsilon, Dimension skeleton_dimension);
+void report_death(std::list<std::pair<double, double> >* intervals, Death d, DistanceType epsilon, DistanceType birthEpsilon, 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& skeleton_dimension,
+ float& eta,
+ float& rho,
+ bool& logscale,
+ 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 skeleton_dimension;
+ float eta;
+ float rho;
+ bool logscale = false;
+ std::string infilename, outfilename;
+ process_command_line_options(argc, argv, skeleton_dimension, eta, rho, logscale, infilename, outfilename);
+
+ // Read in points
+ PointContainer points;
+ read_points(infilename, points);
+
+ std::cout << "Number of points: " << points.size() << std::endl;
+
+ // Create output file
+ std::ofstream out(outfilename.c_str());
+
+ // Create pairwise distances
+ PairDistances distances(points);
+
+ // Create intervals DS
+ std::list<std::pair<double, double> > intervals [skeleton_dimension];
+
+ // 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) <= rho*epsilons[j-1])
+ edges.push_back(std::make_pair(u,v));
+ }
+ std::sort(edges.begin(), edges.end(), RipsGenerator::ComparePair(distances));
+ std::cout << "Total participating edges: " << edges.size() << std::endl;
+ 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 (initial complex is just a disjoint union of 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(),
+ 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");
+ // std::cerr << std::setprecision(15);
+ bool erased[vertices.size()];
+ for (unsigned i = 0; i<vertices.size(); i++)
+ erased[i] = false;
+ 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;
+
+ /* Increase epsilon */
+ // Record the cofaces of all the simplices that need to be removed and reinserted
+ SimplexSet cofaces, large_cofaces, tmp_cofaces;
+ rDebug(" Cofaces size: %d", cofaces.size());
+ // std::cerr << "Cofaces sizes: " << cofaces.size() << " "
+ // << large_cofaces.size() << " "
+ // << tmp_cofaces.size() << std::endl;
+
+ // Add anything else that needs to be inserted into the complex
+ unsigned cee = ce;
+ bool ce_set = false;
+
+ // if (stage > 0) {
+ // std::cerr << "Stage " << stage << " :";
+ // Vertex u,v;
+ // boost::tie(u,v) = edges[cee-1];
+ // std::cerr << distances(u,v) << " <= " << multiplier*epsilons[i];
+ // boost::tie(u,v) = edges[cee];
+ // std::cerr << " < " << distances(u,v) << " <= "
+ // << epsilons[i-1] << std::endl
+ // << " vertices[i] = " << vertices[i] << std::endl
+ // << " vertices[i+1] = " << vertices[i+1] << std::endl;
+ // }
+
+ // if (stage > 0)
+ // std::cerr << "Stage " << stage << " :" << std::endl;
+
+ while (cee < edges.size())
+ {
+ Vertex u,v;
+ boost::tie(u,v) = edges[cee];
+ // if (stage > 0 && (u == vertices[i+1] || v == vertices[i+1]))
+ // std::cerr << "ATTENTION: [" << u << "," << v << "]" << std::endl;
+ // skip already erased edges (may appear since set of participating is larger than in the Morozov zigzag)
+
+ if (distances(u,v) > eta*epsilons[i-1] && !ce_set) {
+ ce = cee;
+ ce_set = true;
+ // std::cerr << "ce = " << distances(u,v) << " > " << eta*epsilons[i-1] << std::endl;
+ }
+ if (distances(u,v) > rho*epsilons[i-1]) {
+ // std::cerr << "cee = " << distances(u,v) << " > " << rho*epsilons[i-1] << std::endl;
+ break;
+ }
+ rDebug(" Recording cofaces of edges[%d]=(%d, %d) with size=%f", (cee-1), u, v, distances(u,v));
+ ec.start();
+ tmp_cofaces.clear();
+
+ // Ignore edges with already removed vertices
+ if (!erased[u] && !erased[v])
+ rips.edge_cofaces(u, v,
+ skeleton_dimension,
+ rho*epsilons[i-1],
+ make_insert_functor(tmp_cofaces),
+ vertices.begin(),
+ vertices.begin() + i + 1);
+ // insert computed cofaces to cofaces sets
+ cofaces.insert(tmp_cofaces.begin(), tmp_cofaces.end());
+ if (distances(u,v) > eta*epsilons[i-1])
+ large_cofaces.insert(tmp_cofaces.begin(), tmp_cofaces.end());
+ ec.stop();
+
+ ++cee;
+ if (cee == edges.size() && !ce_set)
+ ce = cee;
+ }
+ rDebug(" Recorded new cofaces to add");
+
+ // 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));
+ // std::cerr << " Adding " << *cur << ", its size " << size(*cur) << std::endl;
+ 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(intervals, d, epsilons[i-1], epsilons[i], skeleton_dimension);
+ }
+ rInfo("Increased epsilon to %f; complex size: %d", rho*epsilons[i-1], complex.size());
+ // std::cerr << "Increased epsilon to " << rho*epsilons[i-1] << ", complex size " << complex.size() << std::endl;
+ report_memory();
+
+
+ /* Remove edges of length between eta*epsilons[i-1] and rho*epsilons[i-1] and their cofaces */
+ for (SimplexSet::const_reverse_iterator cur = large_cofaces.rbegin(); cur != (SimplexSet::const_reverse_iterator)large_cofaces.rend(); ++cur)
+ {
+ rDebug(" Removing: %s", tostring(*cur).c_str());
+ // std::cerr << " Removing " << *cur << std::endl;
+ 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(intervals, d, epsilons[i-1], epsilons[i], skeleton_dimension);
+ }
+ rInfo("Decreased epsilon to %f; complex size: %d", eta*epsilons[i-1], complex.size());
+ for (Complex::const_iterator cur = complex.begin(); cur != complex.end(); ++cur)
+ rDebug(" %s", tostring(cur->first).c_str());
+ report_memory();
+ // std::cerr << "Decreased epsilon to " << eta*epsilons[i-1] << ", complex size " << complex.size() << std::endl;
+
+
+ /* Remove the vertex */
+ cofaces.clear();
+ rDebug(" Cofaces size: %d", cofaces.size());
+ vc.start();
+ rips.vertex_cofaces(vertices[i],
+ skeleton_dimension,
+ eta*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 != (SimplexSet::const_reverse_iterator)cofaces.rend(); ++cur)
+ {
+ rDebug(" Removing: %s", tostring(*cur).c_str());
+ // std::cerr << " Removing " << *cur << std::endl;
+ 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(intervals, d, epsilons[i-1], epsilons[i], 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());
+ // std::cerr << " Removed vertex " << vertices[i] << ", complex size " << complex.size() << std::endl;
+
+ // Book-keep set of erased vertices
+ erased[vertices[i]] = true;
+
+ 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(intervals, d, epsilons[0], 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 ");
+
+ std::cout << "Writing intervals...";
+ // Note (hack): use epsilons[vertices.size()-2]/2 as minimal value for the log-scale intervals to avoid intervals starting at -infinity and thus a scaling effect
+ write_intervals(out, intervals, skeleton_dimension, logscale, epsilons[vertices.size()-2]/2);
+ std::cout << " done." << std::endl;
+}
+
+
+void write_intervals(std::ostream& out, std::list<std::pair<double,double> >* intervals, int skeleton_dimension, bool logscale, double min_value) {
+ out << "I = { ";
+ for (int d = 0; d<skeleton_dimension; d++) {
+ out << "[ ";
+ for (std::list<std::pair<double,double> >::iterator pit = intervals[d].begin(); pit != intervals[d].end(); pit++)
+ if (logscale)
+ out << "[" << log2(std::max(pit->first, min_value)) << ";" << log2(std::max(pit->second, min_value)) << "] ";
+ else
+ out << "[" << pit->first << ";" << pit->second << "] ";
+
+ // add dummy interval if intervals list is empty
+ if (intervals[d].empty())
+ out << "[0;0] ";
+ out << "] ,";
+ }
+ out << "} ";
+}
+
+void report_death(std::list<std::pair<double,double> >* intervals, Death d, DistanceType epsilon, DistanceType birthEpsilon, Dimension skeleton_dimension)
+{
+ // std::cerr << " d = " << d;
+ // if (d)
+ // std::cerr << " d->dim = " << d->dimension
+ // << " d->dist = " << d->distance;
+ // std::cerr << " epsilon = " << epsilon << std::endl;
+
+ if (d && ((d->distance - epsilon) != 0) && (d->dimension < skeleton_dimension))
+ intervals[d->dimension].push_back(std::pair<double,double>(d->distance, birthEpsilon));
+}
+
+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);
+ }
+}
+
+std::ostream& operator<<(std::ostream& out, const BirthInfo& bi)
+{ return (out << bi.distance); }
+
+
+void process_command_line_options(int argc,
+ char* argv[],
+ unsigned& skeleton_dimension,
+ float& eta,
+ float& rho,
+ bool& logscale,
+ 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")
+ ("dim,d", po::value<unsigned>(&skeleton_dimension)->default_value(3), "maximal dimension of the Rips complex")
+ ("eta,e", po::value<float>(&eta)->default_value(3), "multiplier eta used in the Rips parameter")
+ ("rho,r", po::value<float>(&rho)->default_value(3), "multiplier rho used in the Rips parameter")
+ ( "logscale,l" , po::value<bool>(&logscale)->zero_tokens(), "outputs interval bounds on log_2 scale" ) ;
+#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();
+ }
+}
+
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/homology-zigzags/rips-pairwise.cpp Mon Jun 10 16:58:01 2013 +0200
@@ -0,0 +1,220 @@
+/***************************************************************************
+
+ rips-pairwise: Rips filtration + persistent homology
+ Copyright (C) 2009-2012 Dmitriy Morozov
+
+ Adapted from its original version ("rips-pairwise.cpp" in
+ the Dionysus library) by Steve Oudot (2012).
+
+ Changelog (2012-11-22):
+
+ - the barcode is now output on a log scale and in a format that is
+ compatible with the Matlab layer of the PLEX 2.5 library,
+
+ - the barcode representation is now by closed intervals instead of
+ half-open intervals.
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+***************************************************************************/
+
+#include <topology/rips.h>
+#include <topology/filtration.h>
+#include <topology/static-persistence.h>
+#include <topology/dynamic-persistence.h>
+#include <topology/persistence-diagram.h>
+
+#include <geometry/l2distance.h>
+#include <geometry/distances.h>
+
+#include <utilities/containers.h> // for BackInsertFunctor
+#include <utilities/timer.h>
+
+#include <vector>
+
+#include <boost/program_options.hpp>
+
+
+typedef PairwiseDistances<PointContainer, L2Distance> PairDistances;
+typedef PairDistances::DistanceType DistanceType;
+typedef PairDistances::IndexType Vertex;
+
+typedef Rips<PairDistances> Generator;
+typedef Generator::Simplex Smplx;
+typedef Filtration<Smplx> Fltr;
+// typedef StaticPersistence<> Persistence;
+typedef DynamicPersistenceChains<> Persistence;
+typedef PersistenceDiagram<> PDgm;
+
+// Forward declarations of auxilliary functions
+void write_intervals(std::ostream& out, std::list<std::pair<double,double> >* intervals, int skeleton_dimension);
+void program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, std::string& diagram_name);
+
+int main(int argc, char* argv[])
+{
+ Dimension skeleton;
+ DistanceType max_distance;
+ std::string infilename, diagram_name;
+
+ program_options(argc, argv, infilename, skeleton, max_distance, diagram_name);
+ std::ofstream diagram_out(diagram_name.c_str());
+ std::cout << "Diagram: " << diagram_name << std::endl;
+
+ PointContainer points;
+ read_points(infilename, points);
+
+ PairDistances distances(points);
+ Generator rips(distances);
+ Generator::Evaluator size(distances);
+ Fltr f;
+
+ // Generate 2-skeleton of the Rips complex for epsilon = 50
+ rips.generate(skeleton, max_distance, make_push_back_functor(f));
+ std::cout << "# Generated complex of size: " << f.size() << std::endl;
+
+ // Generate filtration with respect to distance and compute its persistence
+ f.sort(Generator::Comparison(distances));
+
+ Timer persistence_timer; persistence_timer.start();
+ Persistence p(f);
+ p.pair_simplices();
+ persistence_timer.stop();
+
+#if 1
+ // Create intervals DS
+ std::list<std::pair<double, double> > intervals [skeleton];
+
+ // Output cycles
+ Persistence::SimplexMap<Fltr> m = p.make_simplex_map(f);
+ for (Persistence::iterator cur = p.begin(); cur != p.end(); ++cur)
+ {
+ const Persistence::Cycle& cycle = cur->cycle;
+
+ if (!cur->sign()) // only negative simplices have non-empty cycles
+ {
+ Persistence::OrderIndex birth = cur->pair; // the cycle that cur killed was born when we added birth (another simplex)
+
+ const Smplx& b = m[birth];
+ const Smplx& d = m[cur];
+
+ // if (b.dimension() != 1) continue;
+ // std::cout << "Pair: (" << size(b) << ", " << size(d) << ")" << std::endl;
+ if (b.dimension() >= skeleton) continue;
+ // diagram_out << b.dimension() << " " << size(b) << " " << size(d) << std::endl;
+ intervals[b.dimension()].push_back(std::pair<double,double>(size(b), size(d)));
+ } else if (cur->unpaired()) // positive could be unpaired
+ {
+ const Smplx& b = m[cur];
+ // if (b.dimension() != 1) continue;
+
+ // std::cout << "Unpaired birth: " << size(b) << std::endl;
+ // cycle = cur->chain; // TODO
+ if (b.dimension() >= skeleton) continue;
+ // diagram_out << b.dimension() << " " << size(b) << " inf" << std::endl;
+ intervals[b.dimension()].push_back(std::pair<double,double>(size(b), -1));
+ }
+
+ // Iterate over the cycle
+ // for (Persistence::Cycle::const_iterator si = cycle.begin();
+ // si != cycle.end(); ++si)
+ // {
+ // const Smplx& s = m[*si];
+ // //std::cout << s.dimension() << std::endl;
+ // const Smplx::VertexContainer& vertices = s.vertices(); // std::vector<Vertex> where Vertex = Distances::IndexType
+ // AssertMsg(vertices.size() == s.dimension() + 1, "dimension of a simplex is one less than the number of its vertices");
+ // std::cout << vertices[0] << " " << vertices[1] << std::endl;
+ // }
+ }
+#endif
+
+ persistence_timer.check("# Persistence timer");
+
+ std::cout << "Writing intervals...";
+ // Note (hack): use epsilons[vertices.size()-2]/2 as minimal value for the log-scale intervals to avoid intervals starting at -infinity and thus a scaling effect
+ write_intervals(diagram_out, intervals, skeleton);
+ std::cout << " done." << std::endl;
+
+}
+
+
+const double LOG2 = std::log(2.0);
+double log2(double x) {
+ return std::log(x) / LOG2;
+}
+
+void write_intervals(std::ostream& out, std::list<std::pair<double,double> >* intervals, int skeleton_dimension) {
+ out << "I = { ";
+ for (int d = 0; d<skeleton_dimension; d++) {
+ out << "[ ";
+ for (std::list<std::pair<double,double> >::iterator pit = intervals[d].begin(); pit != intervals[d].end(); pit++) {
+ if (pit->first == 0)
+ out << "[-Inf;";
+ else
+ out << "[" << log2(pit->first) << ";";
+ if (pit->second >= 0)
+ out << log2(pit->second) << "] ";
+ else
+ out << "Inf] ";
+ }
+ // add dummy interval if intervals list is empty
+ if (intervals[d].empty())
+ out << "[0;0] ";
+ out << "] ,";
+ }
+ out << "} " << std::endl;
+}
+
+
+void program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, std::string& diagram_name)
+{
+ 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");
+
+ po::options_description visible("Allowed options", 100);
+ visible.add_options()
+ ("help,h", "produce help message")
+ ("skeleton-dimsnion,s", po::value<Dimension>(&skeleton)->default_value(2), "Dimension of the Rips complex we want to compute")
+ ("max-distance,m", po::value<DistanceType>(&max_distance)->default_value(Infinity), "Maximum value for the Rips complex construction")
+ ("diagram,d", po::value<std::string>(&diagram_name), "Filename where to output the persistence diagram");
+#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);
+
+ 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-file"))
+ {
+ std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl;
+ std::cout << visible << std::endl;
+ std::abort();
+ }
+}