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