Added weighted rips code and libraries dev
authorChristos Mantoulidis <cmad@stanford.edu>
Thu, 02 Jul 2009 15:42:44 -0700
branchdev
changeset 150 eb629d8f00bf
parent 136 beff535b53ff
child 151 104ea146b9bc
Added weighted rips code and libraries
examples/rips/CMakeLists.txt
examples/rips/rips-weighted.cpp
include/geometry/weighted-cechdistance.h
include/geometry/weighted-l2distance.h
include/topology/weighted-rips.h
include/topology/weighted-rips.hpp
--- a/examples/rips/CMakeLists.txt	Thu May 14 14:04:43 2009 -0700
+++ b/examples/rips/CMakeLists.txt	Thu Jul 02 15:42:44 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	Thu Jul 02 15:42:44 2009 -0700
@@ -0,0 +1,128 @@
+#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[])
+{
+    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);
+    SimplexVector           complex;
+    
+    // Generate skeleton of the weighted Rips complex for epsilon = 50
+    rips.generate(skeleton, max_distance, make_push_back_functor(complex));
+
+    for (SimplexVector::iterator it = complex.begin(); it != complex.end(); ++it)
+        it->data() = rips.distance(*it, *it);
+
+    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(), DataDimensionComparison<Smplx>());
+
+    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");
+
+    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 (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	Thu Jul 02 15:42:44 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	Thu Jul 02 15:42:44 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__
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/weighted-rips.h	Thu Jul 02 15:42:44 2009 -0700
@@ -0,0 +1,175 @@
+#ifndef __WEIGHTED_RIPS_H__
+#define __WEIGHTED_RIPS_H__
+
+#include <vector>
+#include <string>
+#include "simplex.h"
+#include <boost/iterator/counting_iterator.hpp>
+
+
+/**
+ * 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, typename Distances_::DistanceType> >
+class WeightedRips
+{
+    public:
+        typedef             Distances_                                      Distances; 
+        typedef             typename Distances::IndexType                   IndexType;
+        typedef             typename Distances::DistanceType                DistanceType;
+
+        typedef             Simplex_                                        Simplex;
+        typedef             typename Simplex::Vertex                        Vertex;             // should be the same as IndexType
+        typedef             typename Simplex::VertexContainer               VertexContainer;
+
+        class               Evaluator;
+        class               Comparison;
+        class               ComparePair;       
+
+    public:
+                            WeightedRips(const Distances& distances):
+                                distances_(distances)                       {}
+
+        // Calls functor f on each simplex in the k-skeleton of the weighted Rips complex
+        template<class Functor, class Iterator>
+        void                generate(Dimension k, DistanceType max, const Functor& f, 
+                                     Iterator candidates_begin, Iterator candidates_end) const;
+        
+        // Calls functor f on all the simplices of the weighted Rips complex that contain the given vertex v
+        template<class Functor, class Iterator>
+        void                vertex_cofaces(IndexType v, Dimension k, DistanceType max, const Functor& f, 
+                                           Iterator candidates_begin, Iterator candidates_end) const;
+
+        // Calls functor f on all the simplices of the weighted Rips complex that contain the given edge [u,v]
+        template<class Functor, class Iterator>
+        void                edge_cofaces(IndexType u, IndexType v, Dimension k, DistanceType max, const Functor& f, 
+                                         Iterator candidates_begin, Iterator candidates_end) const;
+        
+        // Calls functor f on all the simplices of the weighted Rips complex that contain the given Simplex s
+        // (unlike the previous methods it does not call the functor on the Simplex s itself)
+        template<class Functor, class Iterator>
+        void                cofaces(const Simplex& s, Dimension k, DistanceType max, const Functor& f,
+                                    Iterator candidates_begin, Iterator candidates_end) const;
+
+        
+        /* No Iterator argument means Iterator = IndexType and the range is [distances().begin(), distances().end()) */
+        template<class Functor>
+        void                generate(Dimension k, DistanceType max, const Functor& f) const
+        { generate(k, max, f, boost::make_counting_iterator(distances().begin()), boost::make_counting_iterator(distances().end())); }
+        
+        template<class Functor>
+        void                vertex_cofaces(IndexType v, Dimension k, DistanceType max, const Functor& f) const
+        { vertex_cofaces(v, k, max, f, boost::make_counting_iterator(distances().begin()), boost::make_counting_iterator(distances().end())); }
+
+        template<class Functor>
+        void                edge_cofaces(IndexType u, IndexType v, Dimension k, DistanceType max, const Functor& f) const
+        { edge_cofaces(u, v, k, max, f, boost::make_counting_iterator(distances().begin()), boost::make_counting_iterator(distances().end())); }
+
+        template<class Functor>
+        void                cofaces(const Simplex& s, Dimension k, DistanceType max, const Functor& f) const
+        { cofaces(s, k, max, f, boost::make_counting_iterator(distances().begin()), boost::make_counting_iterator(distances().end())); }
+
+        
+        const Distances&    distances() const                               { return distances_; }
+        DistanceType        max_distance() const;
+        
+        DistanceType        distance(const Simplex& s1, const Simplex& s2) const;
+
+
+    private:
+        class               WithinDistance;
+
+        template<class Functor, class NeighborTest>
+        void                bron_kerbosch(VertexContainer&                          current, 
+                                          const VertexContainer&                    candidates, 
+                                          typename VertexContainer::const_iterator  excluded,
+                                          Dimension                                 max_dim,
+                                          const NeighborTest&                       neighbor,
+                                          const Functor&                            functor,
+                                          bool                                      check_initial = true) const;
+        
+    private:
+        const Distances&    distances_;
+};
+        
+
+template<class Distances_, class Simplex_>
+class WeightedRips<Distances_, Simplex_>::WithinDistance: public std::binary_function<Vertex, Vertex, bool>
+{
+    public:
+                            WithinDistance(const Distances_&    distances, 
+                                           DistanceType         max):
+                                distances_(distances), max_(max)                        {}
+
+        bool                operator()(Vertex u, Vertex v) const                        { return distances_(u, v) <= max_; }
+
+    private:
+        const Distances&    distances_;  
+        DistanceType        max_;
+};
+
+template<class Distances_, class Simplex_>
+class WeightedRips<Distances_, Simplex_>::Evaluator: public std::unary_function<const Simplex&, DistanceType>
+{
+    public:
+        typedef             Simplex_                                        Simplex;
+
+                            Evaluator(const Distances& distances): 
+                                distances_(distances)                       {}
+
+        DistanceType        operator()(const Simplex& s) const;
+
+    private:
+        const Distances&    distances_;
+};
+
+/* TODO: I'm pretty sure we don't even need this anymore */
+template<class Distances_, class Simplex_>
+class WeightedRips<Distances_, Simplex_>::Comparison: public std::binary_function<const Simplex&, const Simplex&, bool>
+{
+    public:
+        typedef             Simplex_                                        Simplex;
+
+                            Comparison(const Distances& distances):
+                                eval_(distances)                            {}
+
+        bool                operator()(const Simplex& s1, const Simplex& s2) const    
+        { 
+                            if (s1.dimension() != s2.dimension())
+                                return s1.dimension() < s2.dimension();
+                            return s1.data() < s2.data();
+        }
+
+    private:
+        Evaluator           eval_;
+};
+
+/* TODO: May need to change this function for weighted Rips, I'm not sure */
+template<class Distances_, class Simplex_>
+struct WeightedRips<Distances_, Simplex_>::ComparePair: 
+    public std::binary_function<const std::pair<IndexType, IndexType>&,
+                                const std::pair<IndexType, IndexType>&,
+                                bool>
+{
+                            ComparePair(const Distances& distances): 
+                                distances_(distances)                       {}
+
+        bool                operator()(const std::pair<IndexType, IndexType>& a,
+                                       const std::pair<IndexType, IndexType>& b)            {  return   distances_(a.first, a.second) <
+                                                                                                        distances_(b.first, b.second);  }
+
+        const Distances&    distances_;
+};
+
+
+#include "weighted-rips.hpp"
+
+#endif // __WEIGHTED_RIPS_H__
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/weighted-rips.hpp	Thu Jul 02 15:42:44 2009 -0700
@@ -0,0 +1,193 @@
+#include <algorithm>
+#include <utility>
+#include <boost/utility.hpp>
+#include <iostream>
+#include <utilities/log.h>
+#include <utilities/counter.h>
+#include <utilities/indirect.h>
+#include <boost/iterator/counting_iterator.hpp>
+#include <functional>
+
+#ifdef LOGGING
+static rlog::RLogChannel* rlRips =                  DEF_CHANNEL("rips/info", rlog::Log_Debug);  /* TODO: change rips/ to weighted-rips/ ? */
+static rlog::RLogChannel* rlRipsDebug =             DEF_CHANNEL("rips/debug", rlog::Log_Debug); /* TODO: same here */
+#endif // LOGGING
+
+#ifdef COUNTERS
+static Counter*  cClique =                          GetCounter("rips/clique");                  /* TODO: same here */
+#endif // COUNTERS
+
+template<class D, class S>
+template<class Functor, class Iterator>
+void
+WeightedRips<D,S>::
+generate(Dimension k, DistanceType max, const Functor& f, Iterator bg, Iterator end) const
+{
+    rLog(rlRipsDebug,       "Entered generate with %d indices", distances().size());
+
+    WithinDistance neighbor(distances(), max);
+
+    // current      = empty
+    // candidates   = everything
+    VertexContainer current;
+    VertexContainer candidates(bg, end);
+    bron_kerbosch(current, candidates, boost::prior(candidates.begin()), k, neighbor, f);
+}
+
+template<class D, class S>
+template<class Functor, class Iterator>
+void
+WeightedRips<D,S>::
+vertex_cofaces(IndexType v, Dimension k, DistanceType max, const Functor& f, Iterator bg, Iterator end) const
+{
+    WithinDistance neighbor(distances(), max);
+
+    // current      = [v]
+    // candidates   = everything - [v]
+    VertexContainer current; current.push_back(v);
+    VertexContainer candidates;
+    for (Iterator cur = bg; cur != end; ++cur)
+        if (*cur != v && neighbor(v, *cur))
+            candidates.push_back(*cur);
+    bron_kerbosch(current, candidates, boost::prior(candidates.begin()), k, neighbor, f);
+}
+
+template<class D, class S>
+template<class Functor, class Iterator>
+void
+WeightedRips<D,S>::
+edge_cofaces(IndexType u, IndexType v, Dimension k, DistanceType max, const Functor& f, Iterator bg, Iterator end) const
+{
+    rLog(rlRipsDebug,   "In edge_cofaces(%d, %d)", u, v);
+
+    WithinDistance neighbor(distances(), max);
+    AssertMsg(neighbor(u,v), "The new edge must be in the complex");
+
+    // current      = [u,v]
+    // candidates   = everything - [u,v]
+    VertexContainer current; current.push_back(u); current.push_back(v);
+
+    VertexContainer candidates;
+    for (Iterator cur = bg; cur != end; ++cur)
+        if (*cur != u && *cur != v && neighbor(v,*cur) && neighbor(u,*cur))
+        {
+            candidates.push_back(*cur);
+            rLog(rlRipsDebug,   "  added candidate: %d", *cur);
+        }
+
+    bron_kerbosch(current, candidates, boost::prior(candidates.begin()), k, neighbor, f);
+}
+
+template<class D, class S>
+template<class Functor, class Iterator>
+void
+WeightedRips<D,S>::
+cofaces(const Simplex& s, Dimension k, DistanceType max, const Functor& f, Iterator bg, Iterator end) const
+{
+    rLog(rlRipsDebug,   "In cofaces(%s)", tostring(s).c_str());
+
+    WithinDistance neighbor(distances(), max);
+
+    // current      = s.vertices()
+    VertexContainer current(s.vertices().begin(), s.vertices().end());
+    
+    // candidates   = everything - s.vertices()     that is a neighbor() of every vertex in the simplex
+    VertexContainer candidates;
+    typedef difference_iterator<Iterator, 
+                                typename VertexContainer::const_iterator, 
+                                std::less<Vertex> >                     DifferenceIterator;
+    for (DifferenceIterator cur =  DifferenceIterator(bg, end, s.vertices().begin(), s.vertices().end()); 
+                            cur != DifferenceIterator(end, end, s.vertices().end(), s.vertices().end()); 
+                            ++cur)
+    {
+        bool nghbr = true;
+        for (typename VertexContainer::const_iterator v = s.vertices().begin(); v != s.vertices().end(); ++v)
+            if (!neighbor(*v, *cur))        { nghbr = false; break; }
+
+        if (nghbr)  
+        {
+            candidates.push_back(*cur);
+            rLog(rlRipsDebug,   "  added candidate: %d", *cur);
+        }
+    }
+
+    bron_kerbosch(current, candidates, boost::prior(candidates.begin()), k, neighbor, f, false);
+}
+
+
+template<class D, class S>
+template<class Functor, class NeighborTest>
+void
+WeightedRips<D,S>::
+bron_kerbosch(VertexContainer&                          current,    
+              const VertexContainer&                    candidates,     
+              typename VertexContainer::const_iterator  excluded,
+              Dimension                                 max_dim,    
+              const NeighborTest&                       neighbor,       
+              const Functor&                            functor,
+              bool                                      check_initial) const
+{
+    rLog(rlRipsDebug,       "Entered bron_kerbosch");
+    
+    if (check_initial && !current.empty())
+    {
+        Simplex s(current);
+        rLog(rlRipsDebug,   "Reporting simplex: %s", tostring(s).c_str());
+        functor(s);
+    }
+
+    if (current.size() == max_dim + 1)
+        return;
+
+    rLog(rlRipsDebug,       "Traversing %d vertices", candidates.end() - boost::next(excluded));
+    for (typename VertexContainer::const_iterator cur = boost::next(excluded); cur != candidates.end(); ++cur)
+    {
+        current.push_back(*cur);
+        rLog(rlRipsDebug,   "  current.size() = %d, current.back() = %d", current.size(), current.back());
+
+        VertexContainer new_candidates;
+        for (typename VertexContainer::const_iterator ccur = candidates.begin(); ccur != cur; ++ccur)
+            if (neighbor(*ccur, *cur))
+                new_candidates.push_back(*ccur);
+        size_t ex = new_candidates.size();
+        for (typename VertexContainer::const_iterator ccur = boost::next(cur); ccur != candidates.end(); ++ccur)
+            if (neighbor(*ccur, *cur))
+                new_candidates.push_back(*ccur);
+        excluded  = new_candidates.begin() + (ex - 1);
+
+        bron_kerbosch(current, new_candidates, excluded, max_dim, neighbor, functor);
+        current.pop_back();
+    }
+}
+
+template<class Distances_, class Simplex_>
+typename WeightedRips<Distances_, Simplex_>::DistanceType
+WeightedRips<Distances_, Simplex_>::
+distance(const Simplex& s1, const Simplex& s2) const
+{
+    DistanceType mx = 0;
+    for (typename Simplex::VertexContainer::const_iterator      a = s1.vertices().begin();   a != s1.vertices().end();    ++a)
+        for (typename Simplex::VertexContainer::const_iterator  b = s2.vertices().begin();   b != s2.vertices().end();    ++b)
+            mx = std::max(mx, distances_(*a,*b));
+    return mx;
+}
+
+template<class Distances_, class Simplex_>
+typename WeightedRips<Distances_, Simplex_>::DistanceType
+WeightedRips<Distances_, Simplex_>::
+max_distance() const
+{
+    DistanceType mx = 0;
+    for (IndexType a = distances_.begin(); a != distances_.end(); ++a)
+        for (IndexType b = boost::next(a); b != distances_.end(); ++b)
+            mx = std::max(mx, distances_(a,b));
+    return mx;
+}
+
+template<class Distances_, class Simplex_>
+typename WeightedRips<Distances_, Simplex_>::DistanceType
+WeightedRips<Distances_, Simplex_>::Evaluator::
+operator()(const Simplex& s) const
+{
+    return s.data();
+}