Added homology zigzags by Steve Oudot and Don Sheehy dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Mon, 10 Jun 2013 16:58:01 +0200
branchdev
changeset 274 73d69c01ed6c
parent 273 815d6a978c6c
child 275 5f5f4cc70333
Added homology zigzags by Steve Oudot and Don Sheehy
examples/CMakeLists.txt
examples/homology-zigzags/CMakeLists.txt
examples/homology-zigzags/M-ZZ.cpp
examples/homology-zigzags/README
examples/homology-zigzags/dM-ZZ.cpp
examples/homology-zigzags/iR-ZZ.cpp
examples/homology-zigzags/oR-ZZ.cpp
examples/homology-zigzags/rips-pairwise.cpp
--- 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();
+    }
+}