Merged in Chris' changes dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Tue, 04 Aug 2009 17:12:47 -0700
branchdev
changeset 158 ce2aa05994f0
parent 149 3d15aca95dfb (current diff)
parent 157 700cbac5b23c (diff)
child 159 92318b22ffa5
Merged in Chris' changes
--- a/examples/cohomology/CMakeLists.txt	Tue Aug 04 14:28:17 2009 -0700
+++ b/examples/cohomology/CMakeLists.txt	Tue Aug 04 17:12:47 2009 -0700
@@ -1,6 +1,7 @@
 set                         (targets                        
                              rips-cohomology
                              rips-pairwise-cohomology
+                             rips-weighted-cohomology
                              triangle-cohomology
                             )
                              
--- a/examples/cohomology/output.h	Tue Aug 04 14:28:17 2009 -0700
+++ b/examples/cohomology/output.h	Tue Aug 04 17:12:47 2009 -0700
@@ -13,14 +13,16 @@
     return cmp(s1, s2) || cmp(s2, s1);
 }
 
-unsigned index(const SimplexVector& v, const Smplx& s, const Generator::Comparison& cmp)
+template<class Comparison>
+unsigned index(const SimplexVector& v, const Smplx& s, const Comparison& cmp)
 {
     SimplexVector::const_iterator it = std::lower_bound(v.begin(), v.end(), s, cmp);
     while (neq(*it, s)) ++it;
     return it - v.begin();
 }
 
-void output_boundary_matrix(std::ostream& out, const SimplexVector& v, const Generator::Comparison& cmp)
+template<class Comparison>
+void output_boundary_matrix(std::ostream& out, const SimplexVector& v, const Comparison& cmp)
 {
     unsigned i = 0;
     for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
@@ -49,7 +51,7 @@
     }
 }
 
-void output_cocycle(std::string cocycle_prefix, unsigned i, const SimplexVector& v, const Persistence::Cocycle& c, ZpField::Element prime, const Generator::Comparison& cmp)
+void output_cocycle(std::string cocycle_prefix, unsigned i, const SimplexVector& v, const Persistence::Cocycle& c, ZpField::Element prime)
 {
     std::ostringstream istr; istr << '-' << i;
     std::string filename = cocycle_prefix + istr.str() + ".ccl";
@@ -57,8 +59,8 @@
     out << "# Cocycle born at " << c.birth.get<1>() << std::endl;
     for (Persistence::ZColumn::const_iterator zcur = c.zcolumn.begin(); zcur != c.zcolumn.end(); ++zcur)
     {
-        const Smplx& s = **(zcur->si);
+        //const Smplx& s = **(zcur->si);
         out << (zcur->coefficient <= prime/2 ? zcur->coefficient : zcur->coefficient - prime) << " ";
-        out << index(v, s, cmp) << "\n";
+        out << zcur->si->getValue() << "\n";
     }
 }
--- a/examples/cohomology/rips-pairwise-cohomology.cpp	Tue Aug 04 14:28:17 2009 -0700
+++ b/examples/cohomology/rips-pairwise-cohomology.cpp	Tue Aug 04 17:12:47 2009 -0700
@@ -5,6 +5,7 @@
 #include <geometry/distances.h>
 
 #include <utilities/containers.h>           // for BackInsertFunctor
+#include <utilities/property-maps.h>
 #include <utilities/timer.h>
 #include <utilities/log.h>
 
@@ -13,19 +14,22 @@
 #include <boost/tuple/tuple.hpp>
 #include <boost/program_options.hpp>
 
+#include "wrappers.h"
+
 typedef     PairwiseDistances<PointContainer, L2Distance>           PairDistances;
 typedef     PairDistances::DistanceType                             DistanceType;
 typedef     PairDistances::IndexType                                Vertex;
  
-typedef     Rips<PairDistances>                                     Generator;
+typedef     boost::tuple<Dimension, DistanceType>                   BirthInfo;
+typedef     CohomologyPersistence<BirthInfo, Wrapper<unsigned> >    Persistence;
+typedef     Persistence::SimplexIndex                               Index;
+typedef     Persistence::Death                                      Death;
+
+typedef     Rips<PairDistances, Simplex<Vertex, Index> >            Generator;
 typedef     Generator::Simplex                                      Smplx;
 typedef     std::vector<Smplx>                                      SimplexVector;
 typedef     SimplexVector::const_iterator                           SV_const_iterator;
 
-typedef     boost::tuple<Dimension, DistanceType>                   BirthInfo;
-typedef     CohomologyPersistence<BirthInfo, SV_const_iterator>     Persistence;
-typedef     Persistence::SimplexIndex                               Index;
-typedef     Persistence::Death                                      Death;
 typedef     std::map<Smplx, Index, Smplx::VertexComparison>         Complex;
 
 #include "output.h"         // for output_*()
@@ -63,33 +67,41 @@
     Generator::Evaluator    size(distances);
     Generator::Comparison   cmp(distances);
     SimplexVector           v;
-    Complex                 c;
     
     Timer rips_timer; rips_timer.start();
     rips.generate(skeleton, max_distance, make_push_back_functor(v));
-    std::sort(v.begin(), v.end(), cmp);
+    std::sort(v.begin(), v.end(), Smplx::VertexComparison());
+
+    std::vector<unsigned> index_in_v(v.size());
+    for (unsigned idx = 0; idx < v.size(); ++idx)
+        index_in_v[idx] = idx;
+    std::sort(index_in_v.begin(), index_in_v.end(), IndirectIndexComparison<SimplexVector, Generator::Comparison>(v, cmp));
+
+    BinarySearchMap<Smplx, SimplexVector::iterator, Smplx::VertexComparison> map_of_v(v.begin(), v.end());
+
     rips_timer.stop();
     std::cout << "Simplex vector generated, size: " << v.size() << std::endl;
 
-    output_boundary_matrix(bdry_out, v, cmp);
+    output_boundary_matrix(bdry_out, v, Smplx::VertexComparison());
     output_vertex_indices(vertices_out, v);
 
     Timer persistence_timer; persistence_timer.start();
     ZpField                 zp(prime);
     Persistence             p(zp);
-    for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
+    for (unsigned j = 0; j < index_in_v.size(); ++j)
     {
+        SimplexVector::const_iterator cur = v.begin() + index_in_v[j];
         std::vector<Index>      boundary;
         for (Smplx::BoundaryIterator bcur  = cur->boundary_begin(); bcur != cur->boundary_end(); ++bcur)
-            boundary.push_back(c[*bcur]);
+            boundary.push_back(map_of_v[*bcur]->data());
         
         Index idx; Death d;
         bool store = cur->dimension() < skeleton;
-        boost::tie(idx, d)      = p.add(boundary.begin(), boundary.end(), boost::make_tuple(cur->dimension(), size(*cur)), store, cur);
+        boost::tie(idx, d)      = p.add(boundary.begin(), boundary.end(), boost::make_tuple(cur->dimension(), size(*cur)), store, index_in_v[j]);
         
         // c[*cur] = idx;
         if (store)
-            c.insert(std::make_pair(*cur, idx));
+            map_of_v[*cur]->data() = idx;
 
         if (d && (size(*cur) - d->get<1>()) > 0)
         {
@@ -109,7 +121,7 @@
     for (Persistence::Cocycles::const_iterator cur = p.begin(); cur != p.end(); ++cur)
     {
         if (cur->birth.get<0>() != 1) continue;
-        output_cocycle(cocycle_prefix, i, v, *cur, prime, cmp);
+        output_cocycle(cocycle_prefix, i, v, *cur, prime);
         // std::cout << "Cocycle of dimension: " << cur->birth.get<0>() << " born at " << cur->birth.get<1>() << std::endl;
         ++i;
     }
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/cohomology/rips-weighted-cohomology.cpp	Tue Aug 04 17:12:47 2009 -0700
@@ -0,0 +1,184 @@
+#include <topology/cohomology-persistence.h>
+#include <topology/weighted-rips.h>
+
+#include <geometry/weighted-l2distance.h>
+#include <geometry/distances.h>
+
+#include <utilities/containers.h>           // for BackInsertFunctor
+#include <utilities/property-maps.h>
+#include <utilities/timer.h>
+#include <utilities/log.h>
+
+#include <string>
+
+#include <boost/tuple/tuple.hpp>
+#include <boost/program_options.hpp>
+#include <boost/progress.hpp>
+
+#include "wrappers.h"
+
+typedef     PairwiseDistances<PointContainer, WeightedL2Distance>   PairDistances;
+typedef     PairDistances::DistanceType                             DistanceType;
+typedef     PairDistances::IndexType                                Vertex;
+ 
+typedef     boost::tuple<Dimension, DistanceType>                   BirthInfo;
+typedef     CohomologyPersistence<BirthInfo, Wrapper<unsigned> >    Persistence;
+typedef     Persistence::SimplexIndex                               Index;
+typedef     Persistence::Death                                      Death;
+
+typedef     WeightedRips<PairDistances, Simplex<Vertex, Index> >    Generator;
+typedef     Generator::Simplex                                      Smplx;
+typedef     std::vector<Smplx>                                      SimplexVector;
+typedef     SimplexVector::const_iterator                           SV_const_iterator;
+
+#include "output.h"         // for output_*()
+
+void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime, std::string& boundary_name, std::string& cocycle_prefix, std::string& vertices_name, std::string& diagram_name);
+
+int main(int argc, char* argv[])
+{
+#ifdef LOGGING
+    rlog::RLogInit(argc, argv);
+
+    stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+#endif
+
+    Dimension               skeleton;
+    DistanceType            max_distance;
+    ZpField::Element        prime;
+    std::string             infilename, boundary_name, cocycle_prefix, vertices_name, diagram_name;
+
+    program_options(argc, argv, infilename, skeleton, max_distance, prime, boundary_name, cocycle_prefix, vertices_name, diagram_name);
+    std::ofstream           bdry_out(boundary_name.c_str());
+    std::ofstream           vertices_out(vertices_name.c_str());
+    std::ofstream           diagram_out(diagram_name.c_str());
+    std::cout << "Boundary matrix: " << boundary_name << std::endl;
+    std::cout << "Cocycles:        " << cocycle_prefix << "*.ccl" << std::endl;
+    std::cout << "Vertices:        " << vertices_name << std::endl;
+    std::cout << "Diagram:         " << diagram_name << std::endl;
+
+    Timer total_timer; total_timer.start();
+    PointContainer          points;
+    read_weighted_points(infilename, points);
+
+    PairDistances           distances(points);
+    Generator               rips(distances);
+    Generator::Evaluator    size(distances);
+    Generator::Comparison   cmp(distances);
+    SimplexVector           v;
+
+    Timer rips_timer; rips_timer.start();
+    rips.generate(skeleton, max_distance, make_push_back_functor(v));
+
+    /* Keep simplices sorted lexicographically (so that we can binary search through them) */
+    std::sort(v.begin(), v.end(), Smplx::VertexComparison());
+
+    /* We also need the simplices sorted by value though for the filtration:
+       index_in_v[j] refers to the simplex v[index_in_v[j]]                      */
+    std::vector<unsigned> index_in_v(v.size());
+    for (unsigned idx = 0; idx < v.size(); ++idx)
+        index_in_v[idx] = idx;
+    std::sort(index_in_v.begin(), index_in_v.end(), IndirectIndexComparison<SimplexVector, Generator::Comparison>(v, cmp));
+
+    /* Set up map access to the lexicographically sorted simplices */
+    BinarySearchMap<Smplx, SimplexVector::iterator, Smplx::VertexComparison> map_of_v(v.begin(), v.end());
+
+    rips_timer.stop();
+    std::cout << "Simplex vector generated, size: " << v.size() << std::endl;
+
+    /*output_boundary_matrix(bdry_out, v, Smplx::VertexComparison());
+    output_vertex_indices(vertices_out, v);*/
+
+    Timer persistence_timer; persistence_timer.start();
+    ZpField                 zp(prime);
+    Persistence             p(zp);
+    boost::progress_display show_progress(v.size());
+    for (unsigned j = 0; j < index_in_v.size(); ++j)
+    {
+        SimplexVector::const_iterator cur = v.begin() + index_in_v[j];
+        std::vector<Index>      boundary;
+        for (Smplx::BoundaryIterator bcur  = cur->boundary_begin(); bcur != cur->boundary_end(); ++bcur)
+            boundary.push_back(map_of_v[*bcur]->data());
+        
+        Index idx; Death d;
+        bool store = cur->dimension() < skeleton;
+        boost::tie(idx, d)      = p.add(boundary.begin(), boundary.end(), boost::make_tuple(cur->dimension(), size(*cur)), store, index_in_v[j]);
+        
+        if (store)
+            map_of_v[*cur]->data() = idx;
+
+        if (d && (size(*cur) - d->get<1>()) > 0)
+        {
+            AssertMsg(d->get<0>() == cur->dimension() - 1, "Dimensions must match");
+            diagram_out << (cur->dimension() - 1) << " " << d->get<1>() << " " << size(*cur) << std::endl;
+        }
+        ++show_progress;
+    }
+    // output infinte persistence pairs 
+    for (Persistence::CocycleIndex cur = p.begin(); cur != p.end(); ++cur)
+        diagram_out << cur->birth.get<0>() << " " << cur->birth.get<1>() << " inf" << std::endl;
+    persistence_timer.stop();
+
+
+    // p.show_cocycles();
+    // Output alive cocycles of dimension 1
+    unsigned i = 0;
+    for (Persistence::Cocycles::const_iterator cur = p.begin(); cur != p.end(); ++cur)
+    {
+        if (cur->birth.get<0>() != 1) continue;
+        output_cocycle(cocycle_prefix, i, v, *cur, prime);
+        // std::cout << "Cocycle of dimension: " << cur->birth.get<0>() << " born at " << cur->birth.get<1>() << std::endl;
+        ++i;
+    }
+    total_timer.stop();
+    rips_timer.check("Rips timer");
+    persistence_timer.check("Persistence timer");
+    total_timer.check("Total timer");
+}
+
+void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime, std::string& boundary_name, std::string& cocycle_prefix, std::string& vertices_name, 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")
+        ("prime,p",             po::value<ZpField::Element>(&prime)->default_value(11),             "Prime p for the field F_p")
+        ("max-distance,m",      po::value<DistanceType>(&max_distance)->default_value(Infinity),    "Maximum value for the Rips complex construction")
+        ("boundary,b",          po::value<std::string>(&boundary_name),                             "Filename where to output the boundary matrix")
+        ("cocycle,c",           po::value<std::string>(&cocycle_prefix),                            "Prefix of the filename where to output the 1-dimensional cocycles")
+        ("vertices,v",          po::value<std::string>(&vertices_name),                             "Filename where to output the simplex-vertex mapping")
+        ("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();
+    }
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/cohomology/wrappers.h	Tue Aug 04 17:12:47 2009 -0700
@@ -0,0 +1,59 @@
+/**
+ * Wrapper class
+ * 
+ * At points we need to wrap primitive data types in classes, so that we can
+ * pass them as template arguments to classes that like to inherit from their
+ * template arguments--this is done in CohomologyPersistence, for example.
+ */
+
+template<typename Primitive>
+class    Wrapper
+{
+
+    public:
+
+        Wrapper () {}
+        Wrapper (Primitive v)                       { value = v;    }
+
+        void       setValue  (const Primitive &v)   { value = v;    }
+        Primitive &getValue  ()                     { return value; }
+
+        /* provide seemless integration */
+        Wrapper   &operator =(const Primitive &v)   { setValue(v);     return *this; }
+        operator  Primitive()                       { return getValue;               }
+
+    protected:
+
+        Primitive value;
+
+};
+
+/**
+ * IndirectIndexComparison class
+ * 
+ * This class serves as a comparison function for arrays that are being sorted
+ * even though they only contain *indices* to the real data. Therefore, a reference
+ * to the original data as well as the data comparison function needs to be passed
+ * to this class for it to be functional.
+ */
+
+template<class DataContainer, class DataComparison>
+class IndirectIndexComparison
+{
+
+    public:
+
+        IndirectIndexComparison(const DataContainer &dstor, const DataComparison &dcmp) :
+            container(dstor), comparison(dcmp) { }
+
+        bool operator()(const unsigned &idx_1, const unsigned &idx_2) const
+        {
+            return comparison(container[idx_1], container[idx_2]);
+        }
+
+    private:
+
+        const DataContainer  &container;
+        const DataComparison &comparison;
+
+};
--- a/examples/rips/CMakeLists.txt	Tue Aug 04 14:28:17 2009 -0700
+++ b/examples/rips/CMakeLists.txt	Tue Aug 04 17:12:47 2009 -0700
@@ -1,6 +1,7 @@
 set                         (targets                        
                              rips
                              rips-pairwise
+                             rips-weighted
                              rips-image-zigzag
                              rips-zigzag)
                              
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/rips/rips-weighted.cpp	Tue Aug 04 17:12:47 2009 -0700
@@ -0,0 +1,146 @@
+#include <topology/weighted-rips.h>
+#include <topology/filtration.h>
+#include <topology/static-persistence.h>
+#include <topology/dynamic-persistence.h>
+#include <topology/persistence-diagram.h>
+
+//#define RIPS_CLOSURE_CECH_SKELETON
+#ifndef RIPS_CLOSURE_CECH_SKELETON
+#include <geometry/weighted-l2distance.h>
+#else
+#include <geometry/weighted-cechdistance.h>
+#endif
+
+#include <geometry/distances.h>
+
+#include <utilities/containers.h>           // for BackInsertFunctor
+#include <utilities/timer.h>
+
+#include <vector>
+
+#include <boost/program_options.hpp>
+
+
+typedef         PairwiseDistances<PointContainer, WeightedL2Distance>   PairDistances;
+typedef         PairDistances::DistanceType                             DistanceType;
+typedef         PairDistances::IndexType                                Vertex;
+
+typedef         WeightedRips<PairDistances>                             Generator;
+typedef         Generator::Simplex                                      Smplx;
+typedef         std::vector<Smplx>                                      SimplexVector;
+typedef         Filtration<SimplexVector, unsigned>                     Fltr;
+typedef         StaticPersistence<>                                     Persistence;
+//typedef         DynamicPersistenceChains<>                              Persistence;
+typedef         PersistenceDiagram<>                                    PDgm;
+
+void            program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance);
+
+int main(int argc, char* argv[])
+{
+#ifdef LOGGING
+    rlog::RLogInit(argc, argv);
+
+    stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+#endif
+
+    Dimension               skeleton;
+    DistanceType            max_distance;
+    std::string             infilename;
+
+    program_options(argc, argv, infilename, skeleton, max_distance);
+
+    PointContainer          points;
+    read_weighted_points(infilename, points);
+
+    PairDistances           distances(points);
+    Generator               rips(distances);
+    Generator::Evaluator    size(distances);
+    Generator::Comparison   cmp (distances);
+    SimplexVector           complex;
+    
+    // Generate skeleton of the weighted Rips complex for epsilon = 50
+    rips.generate(skeleton, max_distance, make_push_back_functor(complex));
+
+    std::sort(complex.begin(), complex.end(), Smplx::VertexComparison());       // unnecessary
+    std::cout << "# Generated complex of size: " << complex.size() << std::endl;
+
+    // Generate filtration with respect to distance and compute its persistence
+    Fltr f(complex.begin(), complex.end(), cmp);
+
+    Timer persistence_timer; persistence_timer.start();
+    Persistence p(f);
+    p.pair_simplices();
+    persistence_timer.stop();
+
+    // Output cycles
+    for (Persistence::OrderIndex cur = p.begin(); cur != p.end(); ++cur)
+    {
+        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 = f.simplex(f.begin() + (birth - p.begin()));        // eventually this will be encapsulated behind an interface
+            const Smplx& d = f.simplex(f.begin() + (cur   - p.begin()));
+            
+            if (b.dimension() >= skeleton) continue;
+            std::cout << b.dimension() << " " << size(b) << " " << size(d) << std::endl;
+        } else if (cur->pair == cur)    // positive could be unpaired
+        {
+            const Smplx& b = f.simplex(f.begin() + (cur - p.begin()));
+            if (b.dimension() >= skeleton) continue;
+            
+            std::cout << b.dimension() << " " << size(b) << " inf" << std::endl;
+            //cycle = cur->chain;
+        }
+    }
+    
+    persistence_timer.check("# Persistence timer");
+}
+
+void        program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance)
+{
+    namespace po = boost::program_options;
+
+    po::options_description     hidden("Hidden options");
+    hidden.add_options()
+        ("input-file",          po::value<std::string>(&infilename),        "Point set whose weighed 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 weighted Rips complex we want to compute")
+        ("max-distance,m",      po::value<DistanceType>(&max_distance)->default_value(Infinity),    "Maximum value for the weighted Rips complex construction");
+#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()) );
+    /**
+     * Interesting channels
+     * "info", "debug", "topology/persistence"
+     */
+#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();
+    }
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/weighted-cechdistance.h	Tue Aug 04 17:12:47 2009 -0700
@@ -0,0 +1,55 @@
+#ifndef __L2_DISTANCE_H__
+#define __L2_DISTANCE_H__
+
+#include <utilities/types.h>
+
+#include <vector>
+#include <fstream>
+#include <functional>
+#include <cmath>
+#include <string>
+#include <sstream>
+
+
+typedef     std::vector<double>                                     Point;
+typedef     std::vector<Point>                                      PointContainer;
+
+struct WeightedL2Distance:
+    public std::binary_function<const Point&, const Point&, double>
+{
+    result_type     operator()(const Point& p1, const Point& p2) const
+    {
+        AssertMsg(p1.size() == p2.size(), "Points must be in the same dimension (in L2Distance): dim1=%d, dim2=%d", p1.size()-1, p2.size()-1);
+
+        /* the distance of a point to itself is the radius at which the "power distance" ball around it appears:
+           d(p,p) := sqrt(-w_p) */
+        if (p1 == p2)
+            return sqrt(-p1[p1.size()-1]); 
+
+        /* otherwise, the distance is the radius at which the power balls around p, q intersect:
+           d(p,q) := sqrt( ( (w_p-w_q)^2 / 4||p-q||^2 ) - (w_p+w_q) / 2 + ||p-q||^2 / 4 ) */
+        result_type sq_l2dist = 0;
+        result_type wp = p1[p1.size()-1];
+        result_type wq = p2[p2.size()-1];
+        for (size_t i = 0; i < p1.size()-1; ++i)
+            sq_l2dist += (p1[i] - p2[i])*(p1[i] - p2[i]);
+        return sqrt( (wp-wq)*(wp-wq) / (4 * sq_l2dist) - (wp+wq) / 2 + sq_l2dist / 4 );
+    }
+};
+
+void    read_weighted_points(const std::string& infilename, PointContainer& points)
+{
+    std::ifstream in(infilename.c_str());
+    std::string   line;
+    while(std::getline(in, line))
+    {
+        if (line[0] == '#') continue;               // comment line in the file
+        std::stringstream linestream(line);
+        double x;
+        points.push_back(Point());
+        while (linestream >> x)
+            points.back().push_back(x);
+    }
+}
+
+#endif // __L2_DISTANCE_H__
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/weighted-l2distance.h	Tue Aug 04 17:12:47 2009 -0700
@@ -0,0 +1,53 @@
+#ifndef __L2_DISTANCE_H__
+#define __L2_DISTANCE_H__
+
+#include <utilities/types.h>
+
+#include <vector>
+#include <fstream>
+#include <functional>
+#include <cmath>
+#include <string>
+#include <sstream>
+
+
+typedef     std::vector<double>                                     Point;
+typedef     std::vector<Point>                                      PointContainer;
+
+struct WeightedL2Distance:
+    public std::binary_function<const Point&, const Point&, double>
+{
+    result_type     operator()(const Point& p1, const Point& p2) const
+    {
+        AssertMsg(p1.size() == p2.size(), "Points must be in the same dimension (in L2Distance): dim1=%d, dim2=%d", p1.size()-1, p2.size()-1);
+
+        /* the distance of a point to itself is the radius at which the "power distance" ball around it appears:
+           d(p,p) := sqrt(-w_p) */
+        if (p1 == p2)
+            return sqrt(-p1[p1.size()-1]);
+
+        /* otherwise, the distance is the radius at which the power balls around p, q intersect:
+           d(p,q) := sqrt( ||p-q||^2 - w_p - w_q ) */
+        result_type sq_l2dist = 0;
+        for (size_t i = 0; i < p1.size()-1; ++i)
+            sq_l2dist += (p1[i] - p2[i])*(p1[i] - p2[i]);
+        return sqrt(sq_l2dist - p1[p1.size()-1] - p2[p2.size()-1]);
+    }
+};
+
+void    read_weighted_points(const std::string& infilename, PointContainer& points)
+{
+    std::ifstream in(infilename.c_str());
+    std::string   line;
+    while(std::getline(in, line))
+    {
+        if (line[0] == '#') continue;               // comment line in the file
+        std::stringstream linestream(line);
+        double x;
+        points.push_back(Point());
+        while (linestream >> x)
+            points.back().push_back(x);
+    }
+}
+
+#endif // __L2_DISTANCE_H__
--- a/include/topology/rips.h	Tue Aug 04 14:28:17 2009 -0700
+++ b/include/topology/rips.h	Tue Aug 04 17:12:47 2009 -0700
@@ -84,7 +84,7 @@
         DistanceType        distance(const Simplex& s1, const Simplex& s2) const;
 
 
-    private:
+    protected:
         class               WithinDistance;
 
         template<class Functor, class NeighborTest>
@@ -96,7 +96,7 @@
                                           const Functor&                            functor,
                                           bool                                      check_initial = true) const;
         
-    private:
+    protected:
         const Distances&    distances_;
 };
         
@@ -111,7 +111,7 @@
 
         bool                operator()(Vertex u, Vertex v) const                        { return distances_(u, v) <= max_; }
 
-    private:
+    protected:
         const Distances&    distances_;  
         DistanceType        max_;
 };
@@ -127,7 +127,7 @@
 
         DistanceType        operator()(const Simplex& s) const;
 
-    private:
+    protected:
         const Distances&    distances_;
 };
 
@@ -150,7 +150,7 @@
             return e1 < e2;
         }
 
-    private:
+    protected:
         Evaluator           eval_;
 };
 
--- a/include/topology/simplex.h	Tue Aug 04 14:28:17 2009 -0700
+++ b/include/topology/simplex.h	Tue Aug 04 17:12:47 2009 -0700
@@ -127,12 +127,12 @@
         struct DataEvaluator;
         struct DimensionExtractor;
     
-    private:
+    protected:
         VertexContainer&        vertices()                                          { return vdpair_.first(); }
 
         VerticesDataPair        vdpair_;
 
-    private:
+    protected:
         /* Serialization */
         friend class boost::serialization::access;
         
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/weighted-rips.h	Tue Aug 04 17:12:47 2009 -0700
@@ -0,0 +1,139 @@
+#ifndef __WEIGHTED_RIPS_H__
+#define __WEIGHTED_RIPS_H__
+
+#include <vector>
+#include <string>
+#include "simplex.h"
+#include "rips.h"
+#include <boost/iterator/counting_iterator.hpp>
+
+/**
+ * WeightedRipsSimplex class
+ * 
+ * This class sits as an invisible layer between the Simplex datatype passed
+ * to WeightedRips and the class itself. The need for this layer is the need
+ * to store the ``value'' (max inter-vertex distance) of each simplex in the
+ * Weighted Rips complex--something that the user of the class does not need
+ * to be aware of.
+ */
+
+template<class Simplex_, class DistanceType_>
+class WeightedRipsSimplex : public Simplex_
+{
+    public:
+        typedef          typename Simplex_::Vertex                 Vertex;
+        typedef          typename Simplex_::VertexContainer        VertexContainer;
+        typedef          DistanceType_                             DistanceType;
+
+        WeightedRipsSimplex(Simplex_ s) : Simplex_(s)  { }
+
+        void             setSimplexValue(const DistanceType &sv) { simplexValue = sv; }
+        DistanceType     getSimplexValue() const       { return    simplexValue;      }
+
+    protected:
+        DistanceType     simplexValue;
+};
+
+/**
+ * WeightedRips class
+ *
+ * Class providing basic operations to work with Rips complexes. It implements Bron-Kerbosch algorithm, 
+ * and provides simple wrappers for various functions.
+ *
+ * Distances_ is expected to define types IndexType and DistanceType as well as 
+ *               provide operator()(...) which given two IndexTypes should return 
+ *               the distance between them. There should be methods begin() and end() 
+ *               for iterating over IndexTypes as well as a method size().
+ */
+template<class Distances_, class Simplex_ = Simplex<typename Distances_::IndexType> >
+class WeightedRips : public Rips<Distances_, Simplex_>
+{
+    public:
+
+        /* redeclaring the typedefs because they cannot be inherited at compile-time */
+        typedef             Distances_                                      Distances; 
+        typedef             typename Distances::IndexType                   IndexType;
+        typedef             typename Distances::DistanceType                DistanceType;
+
+        typedef             WeightedRipsSimplex<Simplex_, DistanceType>     Simplex;
+        typedef             typename Simplex::Vertex                        Vertex;             // should be the same as IndexType
+        typedef             typename Simplex::VertexContainer               VertexContainer;
+
+        class               Evaluator;
+        class               Comparison;
+
+    public:
+                            WeightedRips(const Distances& distances):
+                                Rips<Distances_, Simplex_>(distances)                             {}
+
+                            template<class Functor>
+                            void generate(Dimension k, DistanceType max, const Functor& f) const;
+
+};
+
+/**
+ * DistanceDataStackingFunctor class
+ * 
+ * Class providing a functor that is to be called by WeightedRips::generate(). This functor
+ * takes as an argument (to its constructor) the original functor passed by the user to
+ * generate(), and a new ``double'' functor is created. Assuming that the functor acts on
+ * simplices, first the value of the simplex is computed (the radius at which the simplex
+ * appears in the weighted Rips complex), the data field of the simplex is populated with
+ * this value, and then the original functor is called (it has no idea that it was
+ * intercepted).
+ */
+
+template<class Rips_, class Functor_>
+class DistanceDataStackingFunctor
+{
+	public:
+		typedef     typename Rips_::Simplex         Simplex_;
+
+		DistanceDataStackingFunctor(const Rips_ &r, const Functor_ &f):
+			rips(r), original_functor(f) { }
+
+		void operator()(const Simplex_ &s) const
+		{
+			Simplex_         s_new(s);
+			s_new.setSimplexValue (rips.distance(s_new, s_new));
+			original_functor      (s_new);
+		}
+
+	private:
+		const Rips_       &rips;
+		const Functor_    &original_functor;
+};
+
+template<class Distances_, class Simplex_>
+template<class Functor>
+void WeightedRips<Distances_, Simplex_>::generate(Dimension k, DistanceType max, const Functor &f) const
+{
+	Rips<Distances_,Simplex_>::generate(k, max, DistanceDataStackingFunctor<WeightedRips<Distances_, Simplex_>,Functor>(*this, f));
+}
+
+template<class Distances_, class Simplex_>
+class WeightedRips<Distances_, Simplex_>::Evaluator: public Rips<Distances_,Simplex_>::Evaluator
+{
+    public:
+                            Evaluator(const Distances& distances): 
+                                Rips<Distances_, Simplex_>::Evaluator(distances)                       {}
+
+        DistanceType       operator()(const Simplex& s) const { return s.getSimplexValue(); }
+};
+
+template<class Distances_, class Simplex_>
+class WeightedRips<Distances_, Simplex_>::Comparison: public Rips<Distances_,Simplex_>::Comparison
+{
+    public:
+                            Comparison(const Distances& distances):
+                                Rips<Distances_, Simplex_>::Comparison(distances)                            {}
+
+        bool                operator()(const Simplex& s1, const Simplex& s2) const    
+        { 
+                            if (s1.dimension() != s2.dimension())
+                                return s1.dimension() < s2.dimension();
+                            return s1.getSimplexValue() < s2.getSimplexValue();
+        }
+};
+
+#endif // __WEIGHTED_RIPS_H__