Added consistency zigzag dev
authorDmitriy Morozov <dmitriy@mrzv.org>
Tue Mar 24 15:17:39 2009 -0700 (2009-03-24)
branchdev
changeset 118c4e25fb4082c
parent 117 06a4361bddaa
child 119 505b3795d239
Added consistency zigzag
examples/CMakeLists.txt
examples/consistency/CMakeLists.txt
examples/consistency/rips-consistency-zigzag.cpp
examples/rips/l2distance.h
examples/rips/rips-zigzag.cpp
include/topology/rips.h
include/topology/rips.hpp
     1.1 --- a/examples/CMakeLists.txt	Mon Mar 23 11:40:29 2009 -0700
     1.2 +++ b/examples/CMakeLists.txt	Tue Mar 24 15:17:39 2009 -0700
     1.3 @@ -1,6 +1,7 @@
     1.4  add_subdirectory			(alphashapes)
     1.5  #add_subdirectory			(ar-vineyard)
     1.6  add_subdirectory			(cech-complex)
     1.7 +add_subdirectory			(consistency)
     1.8  add_subdirectory			(cohomology)
     1.9  add_subdirectory			(fitness)
    1.10  #add_subdirectory			(grid)
     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/examples/consistency/CMakeLists.txt	Tue Mar 24 15:17:39 2009 -0700
     2.3 @@ -0,0 +1,7 @@
     2.4 +set                         (targets                        
     2.5 +                             rips-consistency-zigzag)
     2.6 +                             
     2.7 +foreach                     (t ${targets})
     2.8 +    add_executable          (${t} ${t}.cpp)
     2.9 +    target_link_libraries   (${t} ${libraries} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${Boost_SERIALIZATION_LIBRARY})
    2.10 +endforeach                  (t ${targets})
     3.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     3.2 +++ b/examples/consistency/rips-consistency-zigzag.cpp	Tue Mar 24 15:17:39 2009 -0700
     3.3 @@ -0,0 +1,338 @@
     3.4 +#include <topology/rips.h>
     3.5 +#include <topology/zigzag-persistence.h>
     3.6 +#include <utilities/types.h>
     3.7 +#include <utilities/containers.h>
     3.8 +
     3.9 +#include <utilities/log.h>
    3.10 +#include <utilities/memory.h>       // for report_memory()
    3.11 +#include <utilities/timer.h>
    3.12 +
    3.13 +#include <map>
    3.14 +#include <cmath>
    3.15 +#include <fstream>
    3.16 +#include <sstream>
    3.17 +#include <stack>
    3.18 +#include <cstdlib>
    3.19 +
    3.20 +#include <boost/tuple/tuple.hpp>
    3.21 +#include <boost/program_options.hpp>
    3.22 +#include <boost/progress.hpp>
    3.23 +
    3.24 +#include "../rips/l2distance.h"         // Point, PointContainer, L2DistanceType
    3.25 +
    3.26 +#ifdef COUNTERS
    3.27 +static Counter*  cComplexSize =                     GetCounter("rips/size");
    3.28 +static Counter*  cOperations =                      GetCounter("rips/operations");
    3.29 +#endif // COUNTERS
    3.30 +
    3.31 +typedef     PairwiseDistances<PointContainer, L2Distance>           PairDistances;
    3.32 +typedef     PairDistances::DistanceType                             DistanceType;
    3.33 +
    3.34 +typedef     PairDistances::IndexType                                Vertex;
    3.35 +typedef     std::vector<Vertex>                                     SubsampleVector;
    3.36 +
    3.37 +typedef     Simplex<Vertex>                                         Smplx;
    3.38 +typedef     std::vector<Smplx>                                      SimplexVector;
    3.39 +typedef     std::set<Smplx, Smplx::VertexDimensionComparison>       SimplexSet;
    3.40 +
    3.41 +typedef     std::vector<Vertex>                                     VertexVector;
    3.42 +typedef     std::vector<DistanceType>                               EpsilonVector;
    3.43 +typedef     std::vector<std::pair<Vertex, Vertex> >                 EdgeVector;
    3.44 +
    3.45 +typedef     Rips<PairDistances, Smplx>                              RipsGenerator;
    3.46 +typedef     RipsGenerator::Evaluator                                SimplexEvaluator;
    3.47 +
    3.48 +struct      BirthInfo;
    3.49 +typedef     ZigzagPersistence<BirthInfo>                            Zigzag;
    3.50 +typedef     Zigzag::SimplexIndex                                    Index;
    3.51 +typedef     Zigzag::Death                                           Death;
    3.52 +typedef     std::map<Smplx, Index, 
    3.53 +                            Smplx::VertexDimensionComparison>       Complex;
    3.54 +typedef     Zigzag::ZColumn                                         Boundary;
    3.55 +
    3.56 +// Information we need to know when a class dies
    3.57 +struct      BirthInfo
    3.58 +{
    3.59 +                    BirthInfo(Dimension d = 0, unsigned i = 0, bool u = false):
    3.60 +                        dimension(d), index(i), un(u)               {}
    3.61 +    
    3.62 +    bool            operator<(const BirthInfo& other) const         { if (index == other.index) return (!un && other.un); else return index < other.index; }
    3.63 +    bool            operator>(const BirthInfo& other) const         { return other.operator<(*this); }
    3.64 +
    3.65 +    Dimension       dimension;
    3.66 +    unsigned        index;
    3.67 +    bool            un;
    3.68 +};
    3.69 +std::ostream&       operator<<(std::ostream& out, const BirthInfo& bi)
    3.70 +{ out << "(d=" << bi.dimension << ") " << bi.index; if (bi.un) out << " U " << (bi.index + 1); return out; }
    3.71 +
    3.72 +// Forward declarations of auxilliary functions
    3.73 +void        report_death(std::ostream& out, const BirthInfo& birth, const BirthInfo& death);
    3.74 +void        make_boundary(const Smplx& s, Complex& c, const Zigzag& zz, Boundary& b);
    3.75 +void        process_command_line_options(int           argc,
    3.76 +                                         char*         argv[],
    3.77 +                                         unsigned&     ambient_dimension,
    3.78 +                                         unsigned&     skeleton_dimension,
    3.79 +                                         float&        epsilon,
    3.80 +                                         std::string&  points_filename,
    3.81 +                                         std::string&  subsample_filename,
    3.82 +                                         std::string&  outfilename);
    3.83 +void        add_simplex(const Smplx& s, const BirthInfo& birth, const BirthInfo& death, 
    3.84 +                        Complex& complex, Zigzag& zz, std::ostream& out, Timer& add);
    3.85 +void        remove_simplex(const Smplx& s, const BirthInfo& birth, const BirthInfo& death, 
    3.86 +                           Complex& complex, Zigzag& zz, std::ostream& out, Timer& remove);
    3.87 +
    3.88 +int main(int argc, char* argv[])
    3.89 +{
    3.90 +#ifdef LOGGING
    3.91 +    rlog::RLogInit(argc, argv);
    3.92 +
    3.93 +    stderrLog.subscribeTo( RLOG_CHANNEL("error") );
    3.94 +#endif
    3.95 +
    3.96 +    Timer total, remove, add, coface, generate;
    3.97 +    total.start();
    3.98 + 
    3.99 +#if 0
   3.100 +    SetFrequency(cOperations, 25000);
   3.101 +    SetTrigger(cOperations, cComplexSize);
   3.102 +#endif
   3.103 +
   3.104 +    unsigned        ambient_dimension;
   3.105 +    unsigned        skeleton_dimension;
   3.106 +    float           epsilon;
   3.107 +    std::string     points_filename, subsample_filename, outfilename;
   3.108 +    process_command_line_options(argc, argv, ambient_dimension, skeleton_dimension, epsilon, 
   3.109 +                                             points_filename, subsample_filename, outfilename);
   3.110 +
   3.111 +    // Read in points
   3.112 +    PointContainer      points;
   3.113 +    read_points(points_filename, points, ambient_dimension);
   3.114 +    
   3.115 +    // Read in subsamples
   3.116 +    std::ifstream subsample_in(subsample_filename.c_str());
   3.117 +    std::vector<SubsampleVector> subsamples;
   3.118 +    subsamples.push_back(SubsampleVector());    //  Empty subsample to start from
   3.119 +    while(subsample_in)
   3.120 +    {
   3.121 +        std::string line; std::getline(subsample_in, line);
   3.122 +        std::istringstream line_in(line);
   3.123 +        subsamples.push_back(SubsampleVector());
   3.124 +        unsigned sample;
   3.125 +        while (line_in >> sample)
   3.126 +            subsamples.back().push_back(sample);
   3.127 +    }
   3.128 +    
   3.129 +    std::cout << "Subsample size:" << std::endl;
   3.130 +    for (unsigned i = 0; i < subsamples.size(); ++i)
   3.131 +        std::cout << "  " << subsamples[i].size() << std::endl;
   3.132 +
   3.133 +    // Create output file
   3.134 +    std::ofstream out(outfilename.c_str());
   3.135 +
   3.136 +    // Create pairwise distances
   3.137 +    PairDistances distances(points);
   3.138 +    
   3.139 +
   3.140 +    // Construct zigzag
   3.141 +    Complex             complex;
   3.142 +    Zigzag              zz;
   3.143 +    RipsGenerator       rips(distances);
   3.144 +    SimplexEvaluator    size(distances);
   3.145 +    SimplexVector       subcomplex, across;
   3.146 +    
   3.147 +    rInfo("Commencing computation");
   3.148 +    boost::progress_display show_progress(subsamples.size());
   3.149 +    for (unsigned i = 0; i < subsamples.size() - 1; ++i)
   3.150 +    {
   3.151 +        // Take union of subsamples i and i+1
   3.152 +        SubsampleVector forward_difference, backward_difference;
   3.153 +        std::set_difference(subsamples[i+1].begin(), subsamples[i+1].end(),
   3.154 +                            subsamples[i].begin(),   subsamples[i].end(),
   3.155 +                            std::back_inserter(forward_difference));
   3.156 +        std::set_difference(subsamples[i].begin(),   subsamples[i].end(),
   3.157 +                            subsamples[i+1].begin(), subsamples[i+1].end(),
   3.158 +                            std::back_inserter(backward_difference));
   3.159 +        rInfo("Forward difference size:  %d", forward_difference.size());
   3.160 +        rInfo("Backward difference size: %d", backward_difference.size());
   3.161 +
   3.162 +#if LOGGING
   3.163 +        rDebug("Forward difference:");
   3.164 +        for (SubsampleVector::const_iterator cur = forward_difference.begin(); cur != forward_difference.end(); ++cur)
   3.165 +            rDebug("  %d", *cur);
   3.166 +        
   3.167 +        rDebug("Backward difference:");
   3.168 +        for (SubsampleVector::const_iterator cur = backward_difference.begin(); cur != backward_difference.end(); ++cur)
   3.169 +            rDebug("  %d", *cur);
   3.170 +#endif
   3.171 +
   3.172 +        // Add simplices:       take star of the vertices in subsamples[i+1] - subsamples[i]
   3.173 +        //                      by first computing the Rips complex on those vertices, and then 
   3.174 +        //                      taking the star of each individual simplex (perhaps not the most 
   3.175 +        //                      efficient procedure, but it will do for the first pass)
   3.176 +        rInfo("Adding");
   3.177 +        subcomplex.clear(); across.clear();
   3.178 +        generate.start();
   3.179 +        rips.generate(skeleton_dimension, epsilon, make_push_back_functor(subcomplex), forward_difference.begin(), forward_difference.end());
   3.180 +        std::sort(subcomplex.begin(), subcomplex.end(), Smplx::VertexDimensionComparison());
   3.181 +        generate.stop();
   3.182 +        rInfo("  Subcomplex size: %d", subcomplex.size());
   3.183 +        for (SimplexVector::const_iterator cur = subcomplex.begin(); cur != subcomplex.end(); ++cur)
   3.184 +        {
   3.185 +            add_simplex(*cur, BirthInfo(cur->dimension(), i, true), BirthInfo(cur->dimension() - 1, i), 
   3.186 +                        complex, zz, out, add);
   3.187 +
   3.188 +            // Record cofaces with all other vertices outside of the subcomplex
   3.189 +            coface.start();
   3.190 +            rips.cofaces(*cur, skeleton_dimension, epsilon, make_push_back_functor(across), subsamples[i].begin(), subsamples[i].end());
   3.191 +            coface.stop();
   3.192 +        }
   3.193 +        std::sort(across.begin(), across.end(), Smplx::VertexDimensionComparison());
   3.194 +        rInfo("  Subcomplex simplices added");
   3.195 +        rInfo("  Cross simplices size: %d", across.size());
   3.196 +
   3.197 +        // Add simplices that cut across VR(K_i) and VR(K_{i+1})
   3.198 +        for (SimplexVector::const_iterator cur = across.begin(); cur != across.end(); ++cur)
   3.199 +            add_simplex(*cur, BirthInfo(cur->dimension(), i, true), BirthInfo(cur->dimension() - 1, i), 
   3.200 +                        complex, zz, out, add);
   3.201 +        rInfo("  Cross simplices added");
   3.202 +
   3.203 +
   3.204 +        // Remove simplices:    take star of the vertices in subsamples[i] - subsamples[i+1]
   3.205 +        rInfo("Removing");
   3.206 +        subcomplex.clear(); across.clear();
   3.207 +        generate.start();
   3.208 +        rips.generate(skeleton_dimension, epsilon, make_push_back_functor(subcomplex), backward_difference.begin(), backward_difference.end());
   3.209 +        std::sort(subcomplex.begin(), subcomplex.end(), Smplx::VertexDimensionComparison());
   3.210 +        generate.stop();
   3.211 +        rInfo("  Subcomplex size: %d", subcomplex.size());
   3.212 +
   3.213 +        for (SimplexVector::const_iterator cur = subcomplex.begin(); cur != subcomplex.end(); ++cur)
   3.214 +            rips.cofaces(*cur, skeleton_dimension, epsilon, make_push_back_functor(across), subsamples[i+1].begin(), subsamples[i+1].end());
   3.215 +        std::sort(across.begin(), across.end(), Smplx::VertexDimensionComparison());
   3.216 +        rInfo("  Cross simplices size: %d", across.size());
   3.217 +
   3.218 +        for (SimplexVector::const_reverse_iterator cur = across.rbegin(); cur != across.rend(); ++cur)
   3.219 +            remove_simplex(*cur, BirthInfo(cur->dimension() - 1, i+1), BirthInfo(cur->dimension(), i, true),
   3.220 +                           complex, zz, out, remove);
   3.221 +        rInfo("  Cross simplices removed");
   3.222 +
   3.223 +        for (SimplexVector::const_reverse_iterator cur = subcomplex.rbegin(); cur != subcomplex.rend(); ++cur)
   3.224 +            remove_simplex(*cur, BirthInfo(cur->dimension() - 1, i+1), BirthInfo(cur->dimension(), i, true),
   3.225 +                           complex, zz, out, remove);
   3.226 +        rInfo("  Subcomplex simplices removed");
   3.227 +        
   3.228 +        rInfo("Complex size: %d", complex.size());
   3.229 +        
   3.230 +        ++show_progress;
   3.231 +    }
   3.232 +
   3.233 +    std::cout << std::endl;
   3.234 +    total.stop();
   3.235 +    remove.check("Remove timer          ");
   3.236 +    add.check   ("Add timer             ");
   3.237 +    coface.check("Coface timer          ");
   3.238 +    total.check ("Total timer           ");
   3.239 +}
   3.240 +
   3.241 +
   3.242 +void        make_boundary(const Smplx& s, Complex& c, const Zigzag& zz, Boundary& b)
   3.243 +{
   3.244 +    // rDebug("  Boundary of <%s>", tostring(s).c_str());
   3.245 +    for (Smplx::BoundaryIterator cur = s.boundary_begin(); cur != s.boundary_end(); ++cur)
   3.246 +    {
   3.247 +        b.append(c[*cur], zz.cmp);
   3.248 +        // rDebug("   %d", c[*cur]->order);
   3.249 +    }
   3.250 +}
   3.251 +
   3.252 +void        report_death(std::ostream& out, const BirthInfo& birth, const BirthInfo& death)
   3.253 +{
   3.254 +    if (death > birth)
   3.255 +        out << birth << " --- " << death << std::endl;
   3.256 +}
   3.257 +
   3.258 +void        add_simplex(const Smplx& s, const BirthInfo& birth, const BirthInfo& death, 
   3.259 +                        Complex& complex, Zigzag& zz, std::ostream& out, Timer& add)
   3.260 +{
   3.261 +    rDebug("Adding simplex: %s", tostring(s).c_str());
   3.262 +
   3.263 +    Index idx; Death d; Boundary b;
   3.264 +    make_boundary(s, complex, zz, b);
   3.265 +    add.start();
   3.266 +    boost::tie(idx, d)  = zz.add(b, birth);
   3.267 +    if (d) report_death(out, *d, death);
   3.268 +    add.stop();
   3.269 +    complex.insert(std::make_pair(s, idx));
   3.270 +}
   3.271 +
   3.272 +void        remove_simplex(const Smplx& s, const BirthInfo& birth, const BirthInfo& death, 
   3.273 +                           Complex& complex, Zigzag& zz, std::ostream& out, Timer& remove)
   3.274 +{
   3.275 +    rDebug("Removing simplex: %s", tostring(s).c_str());
   3.276 +
   3.277 +    Complex::iterator si = complex.find(s);
   3.278 +    remove.start();
   3.279 +    Death d = zz.remove(si->second, birth);
   3.280 +    remove.stop();
   3.281 +    if (d) report_death(out, *d, death);
   3.282 +    complex.erase(si);
   3.283 +}
   3.284 +
   3.285 +void        process_command_line_options(int           argc,
   3.286 +                                         char*         argv[],
   3.287 +                                         unsigned&     ambient_dimension,
   3.288 +                                         unsigned&     skeleton_dimension,
   3.289 +                                         float&        epsilon,
   3.290 +                                         std::string&  points_filename,
   3.291 +                                         std::string&  subsample_filename,
   3.292 +                                         std::string&  outfilename)
   3.293 +{
   3.294 +    namespace po = boost::program_options;
   3.295 +
   3.296 +    po::options_description hidden("Hidden options");
   3.297 +    hidden.add_options()
   3.298 +        ("points-file",         po::value<std::string>(&points_filename),   "Point set whose Rips consistency zigzag we want to compute")
   3.299 +        ("subsample-file",      po::value<std::string>(&subsample_filename),"Subsample list")
   3.300 +        ("output-file",         po::value<std::string>(&outfilename),       "Location to save persistence pairs");
   3.301 +    
   3.302 +    po::options_description visible("Allowed options", 100);
   3.303 +    visible.add_options()
   3.304 +        ("help,h",                                                                              "produce help message")
   3.305 +        ("ambient-dimsnion,a",  po::value<unsigned>(&ambient_dimension)->default_value(3),      "The ambient dimension of the point set")
   3.306 +        ("skeleton-dimsnion,s", po::value<unsigned>(&skeleton_dimension)->default_value(2),     "Dimension of the Rips complex we want to compute")
   3.307 +        ("epsilon,e",           po::value<float>(&epsilon)->default_value(0),                   "epsilon to use when computing the Rips complex");
   3.308 +#if LOGGING
   3.309 +    std::vector<std::string>    log_channels;
   3.310 +    visible.add_options()
   3.311 +        ("log,l",               po::value< std::vector<std::string> >(&log_channels),           "log channels to turn on (info, debug, etc)");
   3.312 +#endif
   3.313 +
   3.314 +    po::positional_options_description pos;
   3.315 +    pos.add("points-file", 1);
   3.316 +    pos.add("subsample-file", 1);
   3.317 +    pos.add("output-file", 1);
   3.318 +    
   3.319 +    po::options_description all; all.add(visible).add(hidden);
   3.320 +
   3.321 +    po::variables_map vm;
   3.322 +    po::store(po::command_line_parser(argc, argv).
   3.323 +                  options(all).positional(pos).run(), vm);
   3.324 +    po::notify(vm);
   3.325 +
   3.326 +#if LOGGING
   3.327 +    for (std::vector<std::string>::const_iterator cur = log_channels.begin(); cur != log_channels.end(); ++cur)
   3.328 +        stderrLog.subscribeTo( RLOG_CHANNEL(cur->c_str()) );
   3.329 +    /**
   3.330 +     * Interesting channels
   3.331 +     * "info", "debug", "topology/persistence"
   3.332 +     */
   3.333 +#endif
   3.334 +
   3.335 +    if (vm.count("help") || !vm.count("points-file") || !vm.count("subsample-file") || !vm.count("output-file"))
   3.336 +    { 
   3.337 +        std::cout << "Usage: " << argv[0] << " [options] points-file subsample-file output-file" << std::endl;
   3.338 +        std::cout << visible << std::endl; 
   3.339 +        std::abort();
   3.340 +    }
   3.341 +}
     4.1 --- a/examples/rips/l2distance.h	Mon Mar 23 11:40:29 2009 -0700
     4.2 +++ b/examples/rips/l2distance.h	Tue Mar 24 15:17:39 2009 -0700
     4.3 @@ -17,7 +17,7 @@
     4.4  {
     4.5      result_type     operator()(const Point& p1, const Point& p2) const
     4.6      {
     4.7 -        AssertMsg(p1.size() == p2.size(), "Points must be in the same dimension (in L2Distance)");
     4.8 +        AssertMsg(p1.size() == p2.size(), "Points must be in the same dimension (in L2Distance): dim1=%d, dim2=%d", p1.size(), p2.size());
     4.9          result_type sum = 0;
    4.10          for (size_t i = 0; i < p1.size(); ++i)
    4.11              sum += (p1[i] - p2[i])*(p1[i] - p2[i]);
    4.12 @@ -36,6 +36,7 @@
    4.13          {
    4.14              double      x;
    4.15              in >> x;
    4.16 +            if (!in) { points.pop_back(); break; }
    4.17              points.back().push_back(x);
    4.18          }
    4.19      }
     5.1 --- a/examples/rips/rips-zigzag.cpp	Mon Mar 23 11:40:29 2009 -0700
     5.2 +++ b/examples/rips/rips-zigzag.cpp	Tue Mar 24 15:17:39 2009 -0700
     5.3 @@ -297,21 +297,6 @@
     5.4      }
     5.5  }
     5.6  
     5.7 -bool        face_leaving_subcomplex(Complex::reverse_iterator si, const SimplexEvaluator& size, DistanceType after, DistanceType before)
     5.8 -{
     5.9 -    const Smplx& s = si->first;
    5.10 -    for (Smplx::VertexContainer::const_iterator v1 = s.vertices().begin(); v1 != s.vertices().end(); ++v1)
    5.11 -        for (Smplx::VertexContainer::const_iterator v2 = boost::next(v1); v2 != s.vertices().end(); ++v2)
    5.12 -        {
    5.13 -            Smplx e; e.add(*v1); e.add(*v2);
    5.14 -            if (size(e) > after && size(e) <= before)
    5.15 -                return true;
    5.16 -        }
    5.17 -
    5.18 -    return false;
    5.19 -}
    5.20 -
    5.21 -
    5.22  std::ostream&   operator<<(std::ostream& out, const BirthInfo& bi)
    5.23  { return (out << bi.distance); }
    5.24  
     6.1 --- a/include/topology/rips.h	Mon Mar 23 11:40:29 2009 -0700
     6.2 +++ b/include/topology/rips.h	Tue Mar 24 15:17:39 2009 -0700
     6.3 @@ -53,8 +53,14 @@
     6.4          void                edge_cofaces(IndexType u, IndexType v, Dimension k, DistanceType max, const Functor& f, 
     6.5                                           Iterator candidates_begin, Iterator candidates_end) const;
     6.6          
     6.7 +        // Calls functor f on all the simplices of the Rips complex that contain the given Simplex s
     6.8 +        // (unlike the previous methods it does not call the functor on the Simplex s itself)
     6.9 +        template<class Functor, class Iterator>
    6.10 +        void                cofaces(const Simplex& s, Dimension k, DistanceType max, const Functor& f,
    6.11 +                                    Iterator candidates_begin, Iterator candidates_end) const;
    6.12  
    6.13 -        // No Iterator argument means IndexType and distances().begin() - distances().end()
    6.14 +        
    6.15 +        /* No Iterator argument means Iterator = IndexType and the range is [distances().begin(), distances().end()) */
    6.16          template<class Functor>
    6.17          void                generate(Dimension k, DistanceType max, const Functor& f) const
    6.18          { generate(k, max, f, boost::make_counting_iterator(distances().begin()), boost::make_counting_iterator(distances().end())); }
    6.19 @@ -67,6 +73,10 @@
    6.20          void                edge_cofaces(IndexType u, IndexType v, Dimension k, DistanceType max, const Functor& f) const
    6.21          { edge_cofaces(u, v, k, max, f, boost::make_counting_iterator(distances().begin()), boost::make_counting_iterator(distances().end())); }
    6.22  
    6.23 +        template<class Functor>
    6.24 +        void                cofaces(const Simplex& s, Dimension k, DistanceType max, const Functor& f) const
    6.25 +        { cofaces(s, k, max, f, boost::make_counting_iterator(distances().begin()), boost::make_counting_iterator(distances().end())); }
    6.26 +
    6.27          
    6.28          const Distances&    distances() const                               { return distances_; }
    6.29          DistanceType        max_distance() const;
    6.30 @@ -83,7 +93,8 @@
    6.31                                            typename VertexContainer::const_iterator  excluded,
    6.32                                            Dimension                                 max_dim,
    6.33                                            const NeighborTest&                       neighbor,
    6.34 -                                          const Functor&                            functor) const;
    6.35 +                                          const Functor&                            functor,
    6.36 +                                          bool                                      check_initial = true) const;
    6.37          
    6.38      private:
    6.39          const Distances&    distances_;
     7.1 --- a/include/topology/rips.hpp	Mon Mar 23 11:40:29 2009 -0700
     7.2 +++ b/include/topology/rips.hpp	Tue Mar 24 15:17:39 2009 -0700
     7.3 @@ -4,6 +4,7 @@
     7.4  #include <iostream>
     7.5  #include <utilities/log.h>
     7.6  #include <utilities/counter.h>
     7.7 +#include <utilities/indirect.h>
     7.8  #include <boost/iterator/counting_iterator.hpp>
     7.9  #include <functional>
    7.10  
    7.11 @@ -77,6 +78,42 @@
    7.12      bron_kerbosch(current, candidates, boost::prior(candidates.begin()), k, neighbor, f);
    7.13  }
    7.14  
    7.15 +template<class D, class S>
    7.16 +template<class Functor, class Iterator>
    7.17 +void
    7.18 +Rips<D,S>::
    7.19 +cofaces(const Simplex& s, Dimension k, DistanceType max, const Functor& f, Iterator bg, Iterator end) const
    7.20 +{
    7.21 +    rLog(rlRipsDebug,   "In cofaces(%s)", tostring(s).c_str());
    7.22 +
    7.23 +    WithinDistance neighbor(distances(), max);
    7.24 +
    7.25 +    // current      = s.vertices()
    7.26 +    VertexContainer current(s.vertices().begin(), s.vertices().end());
    7.27 +    
    7.28 +    // candidates   = everything - s.vertices()     that is a neighbor() of every vertex in the simplex
    7.29 +    VertexContainer candidates;
    7.30 +    typedef difference_iterator<Iterator, 
    7.31 +                                typename VertexContainer::const_iterator, 
    7.32 +                                std::less<Vertex> >                     DifferenceIterator;
    7.33 +    for (DifferenceIterator cur =  DifferenceIterator(bg, end, s.vertices().begin(), s.vertices().end()); 
    7.34 +                            cur != DifferenceIterator(end, end, s.vertices().end(), s.vertices().end()); 
    7.35 +                            ++cur)
    7.36 +    {
    7.37 +        bool nghbr = true;
    7.38 +        for (typename VertexContainer::const_iterator v = s.vertices().begin(); v != s.vertices().end(); ++v)
    7.39 +            if (!neighbor(*v, *cur))        { nghbr = false; break; }
    7.40 +
    7.41 +        if (nghbr)  
    7.42 +        {
    7.43 +            candidates.push_back(*cur);
    7.44 +            rLog(rlRipsDebug,   "  added candidate: %d", *cur);
    7.45 +        }
    7.46 +    }
    7.47 +
    7.48 +    bron_kerbosch(current, candidates, boost::prior(candidates.begin()), k, neighbor, f, false);
    7.49 +}
    7.50 +
    7.51  
    7.52  template<class D, class S>
    7.53  template<class Functor, class NeighborTest>
    7.54 @@ -87,11 +124,12 @@
    7.55                typename VertexContainer::const_iterator  excluded,
    7.56                Dimension                                 max_dim,    
    7.57                const NeighborTest&                       neighbor,       
    7.58 -              const Functor&                            functor) const
    7.59 +              const Functor&                            functor,
    7.60 +              bool                                      check_initial) const
    7.61  {
    7.62      rLog(rlRipsDebug,       "Entered bron_kerbosch");
    7.63      
    7.64 -    if (!current.empty())
    7.65 +    if (check_initial && !current.empty())
    7.66      {
    7.67          Simplex s(current);
    7.68          rLog(rlRipsDebug,   "Reporting simplex: %s", tostring(s).c_str());