--- a/bindings/python/CMakeLists.txt Thu Jan 29 10:16:56 2009 -0800
+++ b/bindings/python/CMakeLists.txt Sat Jan 31 23:02:22 2009 -0800
@@ -12,8 +12,6 @@
static-persistence.cpp
simplex.cpp
zigzag-persistence.cpp
-
- test-optional.cpp
)
target_link_libraries (_dionysus ${libraries})
--- a/examples/fitness/avida-rips-distance.cpp Thu Jan 29 10:16:56 2009 -0800
+++ b/examples/fitness/avida-rips-distance.cpp Sat Jan 31 23:02:22 2009 -0800
@@ -9,9 +9,9 @@
typedef ExplicitDistances<AvidaPopulationDetail> ExplicitDist;
-typedef RipsGenerator<ExplicitDist> RipsGen;
+typedef Rips<ExplicitDist> RipsGen;
typedef RipsGen::Simplex Smplx;
-typedef RipsGen::SimplexVector Complex;
+typedef std::vector<Smplx> Complex;
typedef Filtration<Complex, unsigned> Fltr;
typedef StaticPersistence<> Persistence;
@@ -50,7 +50,7 @@
rInfo("Starting to generate rips complex");
Complex c;
- rips.generate(c, 1, rips.max_distance()/2);
+ rips.generate(1, rips.max_distance()/2, make_push_back_functor(c));
std::sort(c.begin(), c.end(), Smplx::VertexComparison());
rInfo("Generated Rips complex, filling filtration");
--- a/examples/rips/CMakeLists.txt Thu Jan 29 10:16:56 2009 -0800
+++ b/examples/rips/CMakeLists.txt Sat Jan 31 23:02:22 2009 -0800
@@ -4,5 +4,5 @@
foreach (t ${targets})
add_executable (${t} ${t}.cpp)
- target_link_libraries (${t} ${libraries} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+ target_link_libraries (${t} ${libraries} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${Boost_SERIALIZATION_LIBRARY})
endforeach (t ${targets})
--- a/examples/rips/rips-zigzag.cpp Thu Jan 29 10:16:56 2009 -0800
+++ b/examples/rips/rips-zigzag.cpp Sat Jan 31 23:02:22 2009 -0800
@@ -22,14 +22,13 @@
typedef std::vector<double> Point;
typedef std::vector<Point> PointContainer;
-struct L2Distance
+struct L2Distance:
+ public std::binary_function<const Point&, const Point&, double>
{
- typedef double value_type;
-
- value_type operator()(const Point& p1, const Point& p2) const
+ result_type operator()(const Point& p1, const Point& p2) const
{
AssertMsg(p1.size() == p2.size(), "Points must be in the same dimension (in L2Distance");
- value_type sum = 0;
+ result_type sum = 0;
for (size_t i = 0; i < p1.size(); ++i)
sum += (p1[i] - p2[i])*(p1[i] - p2[i]);
@@ -42,10 +41,16 @@
typedef PairDistances::IndexType Vertex;
typedef Simplex<Vertex> Smplx;
-typedef RipsBase<PairDistances, Smplx> RipsHelper;
-typedef RipsHelper::Evaluator SimplexEvaluator;
+typedef Rips<PairDistances, Smplx> RipsGenerator;
+typedef RipsGenerator::Evaluator SimplexEvaluator;
-typedef std::pair<DistanceType, Dimension> BirthInfo;
+struct BirthInfo
+{
+ BirthInfo(DistanceType dist = DistanceType(), Dimension dim = Dimension()):
+ distance(dist), dimension(dim) {}
+ DistanceType distance;
+ Dimension dimension;
+};
typedef ImageZigzagPersistence<BirthInfo> Zigzag;
typedef Zigzag::SimplexIndex Index;
typedef Zigzag::Death Death;
@@ -81,13 +86,13 @@
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.second < skeleton)
- std::cout << "Class in the image of dimension: " << cur->birth.second << std::endl;
+ 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.first); }
+{ return (out << bi.distance); }
namespace po = boost::program_options;
@@ -95,9 +100,9 @@
int main(int argc, char* argv[])
{
#ifdef LOGGING
- rlog::RLogInit(argc, argv);
+ rlog::RLogInit(argc, argv);
- stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
+ stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
#endif
SetFrequency(cOperations, 25000);
@@ -106,11 +111,12 @@
unsigned ambient_dimension;
unsigned skeleton_dimension;
float from_multiplier, to_multiplier;
- std::string infilename;
+ std::string infilename, outfilename;
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");
+ ("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()
@@ -127,6 +133,7 @@
po::positional_options_description pos;
pos.add("input-file", 1);
+ pos.add("output-file", 2);
po::options_description all; all.add(visible).add(hidden);
@@ -138,14 +145,14 @@
#if LOGGING
for (std::vector<std::string>::const_iterator cur = log_channels.begin(); cur != log_channels.end(); ++cur)
stdoutLog.subscribeTo( RLOG_CHANNEL(cur->c_str()) );
- /* Interesting channels
+ /* 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 << "Usage: " << argv[0] << " [options] input-file output-file" << std::endl;
std::cout << visible << std::endl;
return 1;
}
@@ -164,6 +171,9 @@
}
}
+ // Create output file
+ std::ofstream out(outfilename.c_str());
+
// Create pairwise distances
PairDistances distances(points);
@@ -198,7 +208,7 @@
// Construct zigzag
Complex complex;
Zigzag zz;
- RipsHelper aux(distances);
+ RipsGenerator aux(distances);
SimplexEvaluator size(distances);
// TODO: it probably makes sense to do things in reverse.
@@ -216,7 +226,7 @@
complex.insert(std::make_pair(sv,
zz.add(Boundary(),
true, // vertex is always in the subcomplex
- std::make_pair(epsilons[i], 0)).first));
+ BirthInfo(epsilons[i], 0)).first));
CountNum(cComplexSize, 0);
Count(cComplexSize);
Count(cOperations);
@@ -246,7 +256,7 @@
Index idx; Death d;
boost::tie(idx, d) = zz.add(b,
(size(s) <= from_multiplier*epsilons[i-1]),
- std::make_pair(epsilons[i-1], s.dimension()));
+ BirthInfo(epsilons[i-1], s.dimension()));
if (!zz.check_consistency())
{
//zz.show_all();
@@ -258,8 +268,8 @@
Count(cOperations);
// Death
- if (d && ((d->first - epsilons[i-1]) != 0) && (d->second < skeleton_dimension))
- std::cout << d->second << " " << d->first << " " << epsilons[i-1] << std::endl;
+ if (d && ((d->distance - epsilons[i-1]) != 0) && (d->dimension < skeleton_dimension))
+ out << d->dimension << " " << d->distance << " " << epsilons[i-1] << std::endl;
}
}
rDebug("Complex after addition:");
@@ -286,10 +296,10 @@
//zz.show_all();
rDebug(" Removing from complex: %s", tostring(si->first).c_str());
Death d = zz.remove(si->second,
- std::make_pair(epsilons[i], si->first.dimension() - 1));
+ BirthInfo(epsilons[i], si->first.dimension() - 1));
AssertMsg(zz.check_consistency(), "Zigzag representation must be consistent after removing a simplex");
- if (d && ((d->first - epsilons[i]) != 0) && (d->second < skeleton_dimension))
- std::cout << d->second << " " << d->first << " " << epsilons[i] << std::endl;
+ if (d && ((d->distance - epsilons[i]) != 0) && (d->dimension < skeleton_dimension))
+ out << d->dimension << " " << d->distance << " " << epsilons[i] << std::endl;
CountNumBy(cComplexSize, si->first.dimension(), -1);
complex.erase(boost::prior(si.base()));
CountBy(cComplexSize, -1);
@@ -299,11 +309,11 @@
// Remove from subcomplex (i.e., remove and reinsert as outside of the subcomplex)
rDebug(" Removing from subcomplex: %s", tostring(si->first).c_str());
Death d = zz.remove(si->second,
- std::make_pair(epsilons[i], si->first.dimension() - 1));
+ BirthInfo(epsilons[i], si->first.dimension() - 1));
Count(cOperations);
AssertMsg(zz.check_consistency(), "Zigzag representation must be consistent after removing a simplex");
- if (d && ((d->first - epsilons[i]) != 0) && (d->second < skeleton_dimension))
- std::cout << d->second << " " << d->first << " " << epsilons[i] << std::endl;
+ if (d && ((d->distance - epsilons[i]) != 0) && (d->dimension < skeleton_dimension))
+ out << d->dimension << " " << d->distance << " " << epsilons[i] << std::endl;
leaving_subcomplex.push(si++);
} else
++si;
@@ -316,11 +326,11 @@
Index idx; Death d;
boost::tie(idx, d) = zz.add(b,
false, // now it is outside of the subcomplex
- std::make_pair(epsilons[i], si->first.dimension()));
+ BirthInfo(epsilons[i], si->first.dimension()));
Count(cOperations);
si->second = idx;
- if (d && ((d->first - epsilons[i]) != 0) && (d->second < skeleton_dimension))
- std::cout << d->second << " " << d->first << " " << epsilons[i] << std::endl;
+ if (d && ((d->distance - epsilons[i]) != 0) && (d->dimension < skeleton_dimension))
+ out << d->dimension << " " << d->distance << " " << epsilons[i] << std::endl;
leaving_subcomplex.pop();
}
}
--- a/examples/rips/rips.cpp Thu Jan 29 10:16:56 2009 -0800
+++ b/examples/rips/rips.cpp Sat Jan 31 23:02:22 2009 -0800
@@ -1,4 +1,12 @@
#include <topology/rips.h>
+#include <topology/filtration.h>
+#include <topology/static-persistence.h>
+#include <topology/persistence-diagram.h>
+#include <utilities/containers.h> // for BackInsertFunctor
+
+#include <fstream>
+#include <boost/archive/binary_oarchive.hpp>
+#include <boost/serialization/map.hpp>
// Trivial example of size() points on a line with integer coordinates
struct Distances
@@ -13,33 +21,56 @@
IndexType end() const { return size(); }
};
+//typedef Rips<ExplicitDistances<Distances> > Generator;
+typedef Rips<Distances> Generator;
+typedef Generator::Simplex Smplx;
+typedef std::vector<Smplx> SimplexVector;
+typedef Filtration<SimplexVector, unsigned> Fltr;
+typedef StaticPersistence<> Persistence;
+typedef PersistenceDiagram<> PDgm;
+
+
int main(int argc, char* argv[])
{
#ifdef LOGGING
rlog::RLogInit(argc, argv);
stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
- stdoutLog.subscribeTo( RLOG_CHANNEL("rips/info") );
+ stdoutLog.subscribeTo( RLOG_CHANNEL("info") );
+ //stdoutLog.subscribeTo( RLOG_CHANNEL("rips/info") );
+ //stdoutLog.subscribeTo( RLOG_CHANNEL("rips/debug") );
#endif
Distances distances;
-#if 0
// Storing ExplicitDistances speeds up the computation (at the price of memory)
- typedef RipsGenerator<ExplicitDistances<Distances> > RipsGenerator;
- ExplicitDistances<Distances> explicit_distances(distances);
- RipsGenerator rips(explicit_distances);
-#else
- //typedef RipsGeneratorMemory<Distances> RipsGenerator;
- typedef RipsGenerator<Distances> RipsGenerator;
- RipsGenerator rips(distances);
-#endif
+ //ExplicitDistances<Distances> explicit_distances(distances);
+
+ Generator rips(distances);
+ Generator::Evaluator size(distances);
+ SimplexVector complex;
+
+ // Generate 2-skeleton of the Rips complex for epsilon = 50
+ rips.generate(2, 50, make_push_back_functor(complex));
+ std::sort(complex.begin(), complex.end(), Smplx::VertexComparison()); // unnecessary
+ rInfo("Generated complex of size: %d", complex.size());
- RipsGenerator::SimplexVector complex;
- //rips.generate(complex, 3, distances.size());
- rips.generate(complex, 3, 50);
-
- std::cout << "Size: " << complex.size() << std::endl;
-// for (RipsGenerator::SimplexVector::const_iterator cur = complex.begin(); cur != complex.end(); ++cur)
-// std::cout << *cur << std::endl;
+ // Generate filtration with respect to distance and compute its persistence
+ Fltr f(complex.begin(), complex.end(), Generator::Comparison(distances));
+ Persistence p(f);
+ p.pair_simplices();
+ rInfo("Simplices paired");
+
+ // Record the persistence intervals in the persistence diagrams
+ std::map<Dimension, PDgm> dgms;
+ init_diagrams(dgms, p.begin(), p.end(),
+ evaluate_through_map(make_offset_map(p.begin(), f.begin()),
+ evaluate_through_filtration(f, size)),
+ evaluate_through_map(make_offset_map(p.begin(), f.begin()),
+ evaluate_through_filtration(f, Smplx::DimensionExtractor())));
+
+ // Serialize the diagrams to a file
+ std::ofstream ofs("rips-diagrams");
+ boost::archive::binary_oarchive oa(ofs);
+ oa << dgms;
}
--- a/include/topology/rips.h Thu Jan 29 10:16:56 2009 -0800
+++ b/include/topology/rips.h Sat Jan 31 23:02:22 2009 -0800
@@ -4,11 +4,14 @@
#include <vector>
#include <string>
#include "simplex.h"
+#include <boost/iterator/counting_iterator.hpp>
+
/**
- * RipsBase class
+ * Rips class
*
- * Base class for the generator of Rips complexes.
+ * 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
@@ -16,7 +19,7 @@
* for iterating over IndexTypes as well as a method size().
*/
template<class Distances_, class Simplex_ = Simplex<typename Distances_::IndexType> >
-class RipsBase
+class Rips
{
public:
typedef Distances_ Distances;
@@ -24,72 +27,86 @@
typedef typename Distances::DistanceType DistanceType;
typedef Simplex_ Simplex;
- typedef std::vector<Simplex> SimplexVector;
+ typedef typename Simplex::Vertex Vertex; // should be the same as IndexType
+ typedef typename Simplex::VertexContainer VertexContainer;
class Evaluator;
class Comparison;
- struct ComparePair;
+ class ComparePair;
public:
- RipsBase(const Distances& distances):
+ Rips(const Distances& distances):
distances_(distances) {}
+ // Calls functor f on each simplex in the k-skeleton of the 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 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 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;
+
+
+ // No Iterator argument means IndexType and 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())); }
+
+
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) const;
+
private:
const Distances& distances_;
};
-template<class Distances_, class Simplex_ = Simplex<typename Distances_::IndexType> >
-class RipsGenerator: public RipsBase<Distances_, Simplex_>
+
+template<class Distances_, class Simplex_>
+class Rips<Distances_, Simplex_>::WithinDistance: public std::binary_function<Vertex, Vertex, bool>
{
public:
- typedef RipsBase<Distances_, Simplex_> Parent;
- typedef typename Parent::Distances Distances;
- typedef typename Parent::Simplex Simplex;
- typedef typename Parent::SimplexVector SimplexVector;
- typedef typename Parent::DistanceType DistanceType;
- typedef typename Parent::IndexType IndexType;
- typedef typename Parent::ComparePair ComparePair;
+ WithinDistance(const Distances_& distances,
+ DistanceType max):
+ distances_(distances), max_(max) {}
- RipsGenerator(const Distances& distances):
- Parent(distances) {}
+ bool operator()(Vertex u, Vertex v) const { return distances_(u, v) <= max_; }
- using Parent::distances;
-
- /// generate k-skeleton of the Rips complex
- void generate(SimplexVector& v, Dimension k, DistanceType max) const;
+ private:
+ const Distances& distances_;
+ DistanceType max_;
};
-// Much more memory efficient, but also much slower
-template<class Distances_, class Simplex_ = Simplex<typename Distances_::IndexType> >
-class RipsGeneratorMemory: public RipsBase<Distances_, Simplex_>
-{
- public:
- typedef RipsBase<Distances_, Simplex_> Parent;
- typedef typename Parent::Distances Distances;
- typedef typename Parent::Simplex Simplex;
- typedef typename Parent::SimplexVector SimplexVector;
- typedef typename Parent::DistanceType DistanceType;
- typedef typename Parent::IndexType IndexType;
- typedef typename Parent::ComparePair ComparePair;
-
- RipsGeneratorMemory(const Distances& distances):
- Parent(distances) {}
-
- using Parent::distances;
- using Parent::distance;
-
- /// generate k-skeleton of the Rips complex
- void generate(SimplexVector& v, Dimension k, DistanceType max) const;
-};
-
-
template<class Distances_, class Simplex_>
-class RipsBase<Distances_, Simplex_>::Evaluator
+class Rips<Distances_, Simplex_>::Evaluator: public std::unary_function<const Simplex&, DistanceType>
{
public:
typedef Simplex_ Simplex;
@@ -104,7 +121,7 @@
};
template<class Distances_, class Simplex_>
-class RipsBase<Distances_, Simplex_>::Comparison
+class Rips<Distances_, Simplex_>::Comparison: public std::binary_function<const Simplex&, const Simplex&, bool>
{
public:
typedef Simplex_ Simplex;
@@ -112,12 +129,37 @@
Comparison(const Distances& distances):
eval_(distances) {}
- bool operator()(const Simplex& s1, const Simplex& s2) const { return eval_(s1) < eval_(s2); }
+ bool operator()(const Simplex& s1, const Simplex& s2) const
+ {
+ DistanceType e1 = eval_(s1),
+ e2 = eval_(s2);
+ if (e1 == e2)
+ return s1.dimension() < s2.dimension();
+
+ return e1 < e2;
+ }
private:
Evaluator eval_;
};
+template<class Distances_, class Simplex_>
+struct Rips<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_;
+};
+
+
/**
* Class: ExplicitDistances
* Stores the pairwise distances of Distances_ instance passed at construction.
@@ -159,7 +201,7 @@
typedef Container_ Container;
typedef Distance_ Distance;
typedef Index_ IndexType;
- typedef typename Distance::value_type DistanceType;
+ typedef typename Distance::result_type DistanceType;
PairwiseDistances(const Container& container,
--- a/include/topology/rips.hpp Thu Jan 29 10:16:56 2009 -0800
+++ b/include/topology/rips.hpp Sat Jan 31 23:02:22 2009 -0800
@@ -4,6 +4,8 @@
#include <iostream>
#include <utilities/log.h>
#include <utilities/counter.h>
+#include <boost/iterator/counting_iterator.hpp>
+#include <functional>
#ifdef LOGGING
static rlog::RLogChannel* rlRips = DEF_CHANNEL("rips/info", rlog::Log_Debug);
@@ -14,129 +16,115 @@
static Counter* cClique = GetCounter("rips/clique");
#endif // COUNTERS
-template<class Distances_, class Simplex_>
-struct RipsBase<Distances_, Simplex_>::ComparePair
-{
- 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_;
-};
-
-template<class DistanceType_, class Simplex_>
+template<class D, class S>
+template<class Functor, class Iterator>
void
-RipsGenerator<DistanceType_, Simplex_>::
-generate(SimplexVector& simplices, Dimension k, DistanceType max) const
+Rips<D,S>::
+generate(Dimension k, DistanceType max, const Functor& f, Iterator bg, Iterator end) const
{
- // Order all the edges
- typedef std::vector< std::pair<IndexType, IndexType> > EdgeVector;
- EdgeVector edges;
- for (IndexType a = distances().begin(); a != distances().end(); ++a)
- {
- Simplex ssx; ssx.add(a);
- simplices.push_back(ssx);
- for (IndexType b = boost::next(a); b != distances().end(); ++b)
- {
- if (distances()(a,b) <= max)
- edges.push_back(std::make_pair(a,b));
- }
- }
- std::sort(edges.begin(), edges.end(), ComparePair(distances()));
+ rLog(rlRipsDebug, "Entered generate with %d indices", distances().size());
+
+ WithinDistance neighbor(distances(), max);
- // Generate simplices
- std::vector<std::vector<size_t> > vertex_star(distances().size());
- for(typename EdgeVector::const_iterator cur = edges.begin(); cur != edges.end(); ++cur)
- {
- rLog(rlRipsDebug, "Current edge: %d %d", cur->first, cur->second);
+ // current = empty
+ // candidates = everything
+ VertexContainer current;
+ VertexContainer candidates(bg, end);
+ bron_kerbosch(current, candidates, boost::prior(candidates.begin()), k, neighbor, f);
+}
- // Create the edge
- Simplex edge; edge.add(cur->first); edge.add(cur->second);
- simplices.push_back(edge);
- if (k <= 1) continue;
-
- vertex_star[cur->first].push_back(simplices.size() - 1);
- vertex_star[cur->second].push_back(simplices.size() - 1);
+template<class D, class S>
+template<class Functor, class Iterator>
+void
+Rips<D,S>::
+vertex_cofaces(IndexType v, Dimension k, DistanceType max, const Functor& f, Iterator bg, Iterator end) const
+{
+ WithinDistance neighbor(distances(), max);
- // Go through a star
- size_t sz = vertex_star[cur->first].size() - 1;
- for (size_t i = 0; i < sz; ++i)
- {
- const Simplex& ssx = simplices[vertex_star[cur->first][i]];
- // FIXME: eventually can uncomment, missing Empty::operator<<()
- // rLog(rlRipsDebug, " %s", tostring(ssx).c_str());
- bool accept = true;
- for (typename Simplex::VertexContainer::const_iterator v = ssx.vertices().begin(); v != ssx.vertices().end(); ++v)
- {
- if (*v == cur->first) continue;
-
- if ( distances()(*v, cur->second) > distances()(cur->first, cur->second) ||
- ((distances()(*v, cur->second) == distances()(cur->first, cur->second)) &&
- (*v > cur->first)))
- {
- accept = false;
- break;
- }
- }
- if (accept)
- {
- Simplex tsx(ssx); tsx.add(cur->second);
- simplices.push_back(tsx);
- // rLog(rlRipsDebug, " Accepting: %s", tostring(tsx).c_str());
-
- // Update stars
- if (tsx.dimension() < k - 1)
- for (typename Simplex::VertexContainer::const_iterator v = static_cast<const Simplex&>(tsx).vertices().begin();
- v != static_cast<const Simplex&>(tsx).vertices().end();
- ++v)
- vertex_star[*v].push_back(simplices.size() - 1);
- }
- }
- }
+ // 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 DistanceType_, class Simplex_>
+template<class D, class S>
+template<class Functor, class Iterator>
void
-RipsGeneratorMemory<DistanceType_, Simplex_>::
-generate(SimplexVector& simplices, Dimension k, DistanceType max) const
+Rips<D,S>::
+edge_cofaces(IndexType u, IndexType v, Dimension k, DistanceType max, const Functor& f, Iterator bg, Iterator end) const
{
- for (IndexType v = distances().begin(); v != distances().end(); ++v)
- {
- simplices.push_back(Simplex());
- simplices.back().add(v);
- }
- size_t last_vertex = simplices.size() - 1;
- size_t begin_previous_dimension = 0;
- size_t end_previous_dimension = simplices.size() - 1;
- typename Simplex::VertexComparison vcmp;
+ 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);
+ }
- for (Dimension d = 1; d < k; ++d)
+ bron_kerbosch(current, candidates, boost::prior(candidates.begin()), k, neighbor, f);
+}
+
+
+template<class D, class S>
+template<class Functor, class NeighborTest>
+void
+Rips<D,S>::
+bron_kerbosch(VertexContainer& current,
+ const VertexContainer& candidates,
+ typename VertexContainer::const_iterator excluded,
+ Dimension max_dim,
+ const NeighborTest& neighbor,
+ const Functor& functor) const
+{
+ rLog(rlRipsDebug, "Entered bron_kerbosch");
+
+ if (!current.empty())
{
- //rLog(rlRips, "Generating dimension %d", d);
- //rLog(rlRips, " Begin previous dimension: %d", begin_previous_dimension);
- //rLog(rlRips, " End previous dimension: %d", end_previous_dimension);
- for (size_t i = 0; i <= last_vertex; ++i)
- {
- for (size_t j = begin_previous_dimension; j <= end_previous_dimension; ++j)
- if (!simplices[j].contains(simplices[i]) &&
- vcmp(simplices[i], simplices[j]) &&
- distance(simplices[i], simplices[j]) <= max)
- {
- simplices.push_back(Simplex(simplices[j]));
- simplices.back().join(simplices[i]);
- }
- }
- begin_previous_dimension = end_previous_dimension + 1;
- end_previous_dimension = simplices.size() - 1;
+ 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 RipsBase<Distances_, Simplex_>::DistanceType
-RipsBase<Distances_, Simplex_>::
+typename Rips<Distances_, Simplex_>::DistanceType
+Rips<Distances_, Simplex_>::
distance(const Simplex& s1, const Simplex& s2) const
{
DistanceType mx = 0;
@@ -147,8 +135,8 @@
}
template<class Distances_, class Simplex_>
-typename RipsBase<Distances_, Simplex_>::DistanceType
-RipsBase<Distances_, Simplex_>::
+typename Rips<Distances_, Simplex_>::DistanceType
+Rips<Distances_, Simplex_>::
max_distance() const
{
DistanceType mx = 0;
@@ -159,8 +147,8 @@
}
template<class Distances_, class Simplex_>
-typename RipsBase<Distances_, Simplex_>::DistanceType
-RipsBase<Distances_, Simplex_>::Evaluator::
+typename Rips<Distances_, Simplex_>::DistanceType
+Rips<Distances_, Simplex_>::Evaluator::
operator()(const Simplex& s) const
{
DistanceType mx = 0;
--- a/include/topology/simplex.h Thu Jan 29 10:16:56 2009 -0800
+++ b/include/topology/simplex.h Sat Jan 31 23:02:22 2009 -0800
@@ -110,7 +110,7 @@
/* Classes: Comparisons
*
* VertexComparison - compare simplices based on the lexicographic ordering of their <vertices()>
- * VertexComparison - compare simplices based on the lexicographic ordering of their <vertices()> within each dimension
+ * VertexDimensionComparison - compare simplices based on the lexicographic ordering of their <vertices()> within each dimension
* DataComparison - compare simplices based on their <data()>
* DataDimensionComparison - compare simplices based on their <data()> within each <dimension()>
*/
--- a/include/utilities/containers.h Thu Jan 29 10:16:56 2009 -0800
+++ b/include/utilities/containers.h Sat Jan 31 23:02:22 2009 -0800
@@ -73,8 +73,8 @@
public SizeStorage<C>
{
typedef CountingBackInserter Self;
- typedef std::back_insert_iterator<C> ParentIterator;
- typedef SizeStorage<C> ParentSize;
+ typedef std::back_insert_iterator<C> ParentIterator;
+ typedef SizeStorage<C> ParentSize;
CountingBackInserter(C& c):
ParentIterator(c) {}
@@ -83,6 +83,56 @@
Self operator++(int i) { Self tmp = *this; ParentSize::operator++(i); ParentIterator::operator++(i); return tmp; }
};
+/**
+ * Class: PushBackFunctor<Container>
+ *
+ * Performs the same task as std::back_insert_iterator<Container>, but as a functor.
+ */
+template<class Container_>
+class PushBackFunctor
+{
+ public:
+ typedef Container_ Container;
+ typedef typename Container::value_type value_type;
+
+ PushBackFunctor(Container& container):
+ container_(container) {}
+
+ void operator()(const value_type& v) const { container_.push_back(v); }
+
+ private:
+ Container& container_;
+};
+
+template<class Container>
+PushBackFunctor<Container>
+make_push_back_functor(Container& container) { return PushBackFunctor<Container>(container); }
+
+/**
+ * Class: InsertFunctor<Container>
+ *
+ * Performs insertions of its arguments into the given container.
+ */
+template<class Container_>
+class InsertFunctor
+{
+ public:
+ typedef Container_ Container;
+ typedef typename Container::value_type value_type;
+
+ InsertFunctor(Container& container):
+ container_(container) {}
+
+ void operator()(const value_type& v) const { container_.insert(v); }
+
+ private:
+ Container& container_;
+};
+
+template<class Container>
+InsertFunctor<Container>
+make_insert_functor(Container& container) { return InsertFunctor<Container>(container); }
+
/* Specializations */
template<class T, class Comparison_>