--- 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(); + } +}