Modified cohomology code (unweighted and weighted) so that it doesn't use std::map for boundary computations, which should optimize memory usage.
--- a/examples/cohomology/CMakeLists.txt Tue Aug 04 13:23:16 2009 -0700
+++ b/examples/cohomology/CMakeLists.txt Tue Aug 04 16:58:44 2009 -0700
@@ -1,7 +1,7 @@
set (targets
rips-cohomology
rips-pairwise-cohomology
- rips-weighted-cohomology
+ rips-weighted-cohomology
triangle-cohomology
)
--- a/examples/cohomology/output.h Tue Aug 04 13:23:16 2009 -0700
+++ b/examples/cohomology/output.h Tue Aug 04 16:58:44 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 13:23:16 2009 -0700
+++ b/examples/cohomology/rips-pairwise-cohomology.cpp Tue Aug 04 16:58:44 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;
}
--- a/examples/cohomology/rips-weighted-cohomology.cpp Tue Aug 04 13:23:16 2009 -0700
+++ b/examples/cohomology/rips-weighted-cohomology.cpp Tue Aug 04 16:58:44 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>
@@ -14,21 +15,22 @@
#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 WeightedRips<PairDistances> Generator;
+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;
-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_*()
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);
@@ -64,41 +66,53 @@
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);
+
+ /* 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, 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);
boost::progress_display show_progress(v.size());
- 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)
{
AssertMsg(d->get<0>() == cur->dimension() - 1, "Dimensions must match");
diagram_out << (cur->dimension() - 1) << " " << d->get<1>() << " " << size(*cur) << std::endl;
}
- ++show_progress;
+ ++show_progress;
}
// output infinte persistence pairs
for (Persistence::CocycleIndex cur = p.begin(); cur != p.end(); ++cur)
@@ -112,7 +126,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/wrappers.h Tue Aug 04 16:58:44 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;
+
+};