--- a/examples/CMakeLists.txt Sun Mar 22 11:18:07 2009 -0700
+++ b/examples/CMakeLists.txt Mon Mar 23 11:40:29 2009 -0700
@@ -1,6 +1,7 @@
add_subdirectory (alphashapes)
#add_subdirectory (ar-vineyard)
add_subdirectory (cech-complex)
+add_subdirectory (cohomology)
add_subdirectory (fitness)
#add_subdirectory (grid)
add_subdirectory (triangle)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/cohomology/CMakeLists.txt Mon Mar 23 11:40:29 2009 -0700
@@ -0,0 +1,10 @@
+set (targets
+ rips-cohomology
+ rips-pairwise-cohomology
+ triangle-cohomology
+ )
+
+foreach (t ${targets})
+ add_executable (${t} ${t}.cpp)
+ target_link_libraries (${t} ${libraries} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+endforeach (t ${targets})
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/cohomology/rips-cohomology.cpp Mon Mar 23 11:40:29 2009 -0700
@@ -0,0 +1,61 @@
+#include <topology/cohomology-persistence.h>
+#include <topology/rips.h>
+
+#include <utilities/containers.h> // for BackInsertFunctor
+
+#include <map>
+#include <iostream>
+
+#include <boost/tuple/tuple.hpp>
+
+// Trivial example of size() points on a line with integer coordinates
+struct Distances
+{
+ typedef int IndexType;
+ typedef double DistanceType;
+
+ DistanceType operator()(IndexType a, IndexType b) const { return std::abs(a - b); }
+
+ size_t size() const { return 2000; }
+ IndexType begin() const { return 0; }
+ IndexType end() const { return size(); }
+};
+
+typedef CohomologyPersistence<Distances::DistanceType> Persistence;
+typedef Persistence::SimplexIndex Index;
+typedef Persistence::Death Death;
+
+typedef Rips<Distances> Generator;
+typedef Generator::Simplex Smplx;
+
+typedef std::map<Smplx, Index,
+ Smplx::VertexComparison> Complex;
+typedef std::vector<Smplx> SimplexVector;
+
+int main()
+{
+ Distances distances;
+ Generator rips(distances);
+ Generator::Evaluator size(distances);
+ SimplexVector v;
+ Complex c;
+
+ rips.generate(2, 50, make_push_back_functor(v));
+ std::sort(v.begin(), v.end(), Generator::Comparison(distances));
+ std::cout << "Simplex vector generated, size: " << v.size() << std::endl;
+
+ Persistence p;
+ for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
+ {
+ std::vector<Index> boundary;
+ for (Smplx::BoundaryIterator bcur = cur->boundary_begin();
+ bcur != cur->boundary_end(); ++bcur)
+ boundary.push_back(c[*bcur]);
+
+ Index idx; Death d;
+ boost::tie(idx, d) = p.add(boundary.begin(), boundary.end(), size(*cur));
+ c[*cur] = idx;
+ if (d && (size(*cur) - *d) > 0)
+ std::cout << (cur->dimension() - 1) << " " << *d << " " << size(*cur) << std::endl;
+ }
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/cohomology/rips-pairwise-cohomology.cpp Mon Mar 23 11:40:29 2009 -0700
@@ -0,0 +1,103 @@
+#include <topology/cohomology-persistence.h>
+#include <topology/rips.h>
+
+#include <utilities/containers.h> // for BackInsertFunctor
+#include <utilities/timer.h>
+
+#include <string>
+
+#include <boost/tuple/tuple.hpp>
+#include <boost/program_options.hpp>
+
+#include "../rips/l2distance.h"
+
+typedef PairwiseDistances<PointContainer, L2Distance> PairDistances;
+typedef PairDistances::DistanceType DistanceType;
+typedef PairDistances::IndexType Vertex;
+
+typedef CohomologyPersistence<DistanceType> Persistence;
+typedef Persistence::SimplexIndex Index;
+typedef Persistence::Death Death;
+
+typedef Rips<PairDistances> Generator;
+typedef Generator::Simplex Smplx;
+
+typedef std::map<Smplx, Index,
+ Smplx::VertexComparison> Complex;
+typedef std::vector<Smplx> SimplexVector;
+
+void program_options(int argc, char* argv[], std::string& infilename, Dimension& ambient, Dimension& skeleton, DistanceType& max_distance);
+
+int main(int argc, char* argv[])
+{
+ Dimension ambient, skeleton;
+ DistanceType max_distance;
+ std::string infilename;
+
+ program_options(argc, argv, infilename, ambient, skeleton, max_distance);
+
+ PointContainer points;
+ read_points(infilename, points, ambient);
+
+ PairDistances distances(points);
+ Generator rips(distances);
+ Generator::Evaluator size(distances);
+ SimplexVector v;
+ Complex c;
+
+ rips.generate(skeleton, max_distance, make_push_back_functor(v));
+ std::sort(v.begin(), v.end(), Generator::Comparison(distances));
+ std::cout << "Simplex vector generated, size: " << v.size() << std::endl;
+
+ Timer persistence_timer; persistence_timer.start();
+ Persistence p;
+ for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
+ {
+ std::vector<Index> boundary;
+ for (Smplx::BoundaryIterator bcur = cur->boundary_begin();
+ bcur != cur->boundary_end(); ++bcur)
+ boundary.push_back(c[*bcur]);
+
+ Index idx; Death d;
+ boost::tie(idx, d) = p.add(boundary.begin(), boundary.end(), size(*cur));
+ c[*cur] = idx;
+ if (d && (size(*cur) - *d) > 0)
+ std::cout << (cur->dimension() - 1) << " " << *d << " " << size(*cur) << std::endl;
+ }
+ persistence_timer.stop();
+ persistence_timer.check("Persistence timer");
+}
+
+void program_options(int argc, char* argv[], std::string& infilename, Dimension& ambient, 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 Rips zigzag we want to compute");
+
+ po::options_description visible("Allowed options", 100);
+ visible.add_options()
+ ("help,h", "produce help message")
+ ("ambient-dimsnion,a", po::value<Dimension>(&ambient)->default_value(3), "The ambient dimension of the point set")
+ ("skeleton-dimsnion,s", po::value<Dimension>(&skeleton)->default_value(2), "Dimension of the Rips complex we want to compute")
+ ("max-distance,m", po::value<DistanceType>(&max_distance)->default_value(Infinity), "Maximum value for the Rips complex construction");
+
+ po::positional_options_description pos;
+ pos.add("input-file", 1);
+ pos.add("output-file", 2);
+
+ po::options_description all; all.add(visible).add(hidden);
+
+ po::variables_map vm;
+ po::store(po::command_line_parser(argc, argv).
+ options(all).positional(pos).run(), vm);
+ po::notify(vm);
+
+ if (vm.count("help") || !vm.count("input-file"))
+ {
+ std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl;
+ std::cout << visible << std::endl;
+ std::abort();
+ }
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/cohomology/triangle-cohomology.cpp Mon Mar 23 11:40:29 2009 -0700
@@ -0,0 +1,90 @@
+#include <topology/simplex.h>
+#include <topology/cohomology-persistence.h>
+
+#include <utilities/log.h>
+#include <utilities/indirect.h>
+
+#include <vector>
+#include <map>
+#include <iostream>
+
+#include <boost/tuple/tuple.hpp>
+
+#if 1
+#include <fstream>
+#include <boost/archive/text_oarchive.hpp>
+#include <boost/archive/text_iarchive.hpp>
+#include <boost/serialization/vector.hpp>
+#endif
+
+typedef CohomologyPersistence<unsigned> Persistence;
+typedef Persistence::SimplexIndex Index;
+typedef Persistence::Death Death;
+
+typedef unsigned Vertex;
+typedef Simplex<Vertex, double> Smplx;
+
+typedef std::map<Smplx, Index,
+ Smplx::VertexComparison> Complex;
+typedef std::vector<Smplx> SimplexVector;
+
+
+void fillTriangleSimplices(SimplexVector& c)
+{
+ typedef std::vector<Vertex> VertexVector;
+ VertexVector vertices(4);
+ vertices[0] = 0; vertices[1] = 1; vertices[2] = 2;
+ vertices[3] = 0;
+
+ VertexVector::const_iterator bg = vertices.begin();
+ VertexVector::const_iterator end = vertices.end();
+ c.push_back(Smplx(bg, bg + 1, 0)); // 0 = A
+ c.push_back(Smplx(bg + 1, bg + 2, 1)); // 1 = B
+ c.push_back(Smplx(bg + 2, bg + 3, 2)); // 2 = C
+ c.push_back(Smplx(bg, bg + 2, 2.5)); // AB
+ c.push_back(Smplx(bg + 1, bg + 3, 2.9)); // BC
+ c.push_back(Smplx(bg + 2, end, 3.5)); // CA
+ c.push_back(Smplx(bg, bg + 3, 5)); // ABC
+}
+
+int main(int argc, char** argv)
+{
+#ifdef LOGGING
+ rlog::RLogInit(argc, argv);
+
+ stderrLog.subscribeTo(RLOG_CHANNEL("info"));
+ stderrLog.subscribeTo(RLOG_CHANNEL("error"));
+ //stderrLog.subscribeTo(RLOG_CHANNEL("topology"));
+#endif
+
+ SimplexVector v;
+ fillTriangleSimplices(v);
+ std::sort(v.begin(), v.end(), Smplx::DataComparison());
+ std::cout << "Simplices filled" << std::endl;
+ for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
+ std::cout << " " << *cur << std::endl;
+
+ // Compute persistence
+ Complex c;
+ Persistence p;
+ unsigned i = 0;
+ for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
+ {
+ std::vector<Index> boundary;
+ for (Smplx::BoundaryIterator bcur = cur->boundary_begin();
+ bcur != cur->boundary_end(); ++bcur)
+ boundary.push_back(c[*bcur]);
+
+ Index idx; Death d;
+ boost::tie(idx, d) = p.add(boundary.begin(), boundary.end(), i++);
+ c[*cur] = idx;
+ if (d)
+ std::cout << (cur->dimension() - 1) << " " << *d << " " << (i-1) << std::endl;
+ // the dimension above is adjusted for what it would look like in homology
+ // (i.e. when a 1 class kills 0, it's really that in cohomology forward 0 kills 1,
+ // in cohomology backward 1 kills 0, and in homology 1 kills 0)
+
+ //p.show_cycles();
+ }
+}
+
--- a/examples/rips/CMakeLists.txt Sun Mar 22 11:18:07 2009 -0700
+++ b/examples/rips/CMakeLists.txt Mon Mar 23 11:40:29 2009 -0700
@@ -1,5 +1,6 @@
set (targets
rips
+ rips-pairwise
rips-image-zigzag
rips-zigzag)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/rips/l2distance.h Mon Mar 23 11:40:29 2009 -0700
@@ -0,0 +1,44 @@
+#ifndef __L2_DISTANCE_H__
+#define __L2_DISTANCE_H__
+
+#include <utilities/types.h>
+
+#include <vector>
+#include <fstream>
+#include <functional>
+#include <cmath>
+
+
+typedef std::vector<double> Point;
+typedef std::vector<Point> PointContainer;
+
+struct L2Distance:
+ 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)");
+ result_type sum = 0;
+ for (size_t i = 0; i < p1.size(); ++i)
+ sum += (p1[i] - p2[i])*(p1[i] - p2[i]);
+
+ return sqrt(sum);
+ }
+};
+
+void read_points(const std::string& infilename, PointContainer& points, Dimension ambient)
+{
+ std::ifstream in(infilename.c_str());
+ while(in)
+ {
+ points.push_back(Point());
+ for (unsigned i = 0; i < ambient; ++i)
+ {
+ double x;
+ in >> x;
+ points.back().push_back(x);
+ }
+ }
+}
+
+#endif // __L2_DISTANCE_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/rips/rips-pairwise.cpp Mon Mar 23 11:40:29 2009 -0700
@@ -0,0 +1,129 @@
+#include <topology/rips.h>
+#include <topology/filtration.h>
+#include <topology/static-persistence.h>
+#include <topology/dynamic-persistence.h>
+#include <topology/persistence-diagram.h>
+
+#include <utilities/containers.h> // for BackInsertFunctor
+#include <utilities/timer.h>
+
+#include <vector>
+
+#include <boost/program_options.hpp>
+
+#include "l2distance.h"
+
+typedef PairwiseDistances<PointContainer, L2Distance> PairDistances;
+typedef PairDistances::DistanceType DistanceType;
+typedef PairDistances::IndexType Vertex;
+
+typedef Rips<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& ambient, Dimension& skeleton, DistanceType& max_distance);
+
+int main(int argc, char* argv[])
+{
+ Dimension ambient, skeleton;
+ DistanceType max_distance;
+ std::string infilename;
+
+ program_options(argc, argv, infilename, ambient, skeleton, max_distance);
+
+ PointContainer points;
+ read_points(infilename, points, ambient);
+
+ PairDistances distances(points);
+ Generator rips(distances);
+ Generator::Evaluator size(distances);
+ SimplexVector complex;
+
+ // Generate 2-skeleton of the 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(), Generator::Comparison(distances));
+
+ 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::OrderDescriptor::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() != 1) continue;
+ std::cout << "Pair: (" << 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() != 1) continue;
+
+ std::cout << "Unpaired birth: " << size(b) << std::endl;
+ cycle = cur->chain;
+ }
+
+ // Iterate over the cycle
+ for (Persistence::OrderDescriptor::Cycle::const_iterator si = cycle.begin();
+ si != cycle.end(); ++si)
+ {
+ const Smplx& s = f.simplex(f.begin() + (*si - p.begin()));
+ //std::cout << s.dimension() << std::endl;
+ const Smplx::VertexContainer& vertices = s.vertices(); // std::vector<Vertex> where Vertex = Distances::IndexType
+ AssertMsg(vertices.size() == s.dimension() + 1, "dimension of a simplex is one less than the number of its vertices");
+ std::cout << vertices[0] << " " << vertices[1] << std::endl;
+ }
+ }
+
+ persistence_timer.check("Persistence timer");
+}
+
+void program_options(int argc, char* argv[], std::string& infilename, Dimension& ambient, 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 Rips zigzag we want to compute");
+
+ po::options_description visible("Allowed options", 100);
+ visible.add_options()
+ ("help,h", "produce help message")
+ ("ambient-dimsnion,a", po::value<Dimension>(&ambient)->default_value(3), "The ambient dimension of the point set")
+ ("skeleton-dimsnion,s", po::value<Dimension>(&skeleton)->default_value(2), "Dimension of the Rips complex we want to compute")
+ ("max-distance,m", po::value<DistanceType>(&max_distance)->default_value(Infinity), "Maximum value for the Rips complex construction");
+
+ po::positional_options_description pos;
+ pos.add("input-file", 1);
+ pos.add("output-file", 2);
+
+ po::options_description all; all.add(visible).add(hidden);
+
+ po::variables_map vm;
+ po::store(po::command_line_parser(argc, argv).
+ options(all).positional(pos).run(), vm);
+ po::notify(vm);
+
+ if (vm.count("help") || !vm.count("input-file"))
+ {
+ std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl;
+ std::cout << visible << std::endl;
+ std::abort();
+ }
+}
--- a/examples/rips/rips-zigzag.cpp Sun Mar 22 11:18:07 2009 -0700
+++ b/examples/rips/rips-zigzag.cpp Mon Mar 23 11:40:29 2009 -0700
@@ -17,27 +17,13 @@
#include <boost/program_options.hpp>
#include <boost/progress.hpp>
+#include "l2distance.h" // Point, PointContainer, L2DistanceType
#ifdef COUNTERS
static Counter* cComplexSize = GetCounter("rips/size");
static Counter* cOperations = GetCounter("rips/operations");
#endif // COUNTERS
-typedef std::vector<double> Point;
-typedef std::vector<Point> PointContainer;
-struct L2Distance:
- 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)");
- result_type sum = 0;
- for (size_t i = 0; i < p1.size(); ++i)
- sum += (p1[i] - p2[i])*(p1[i] - p2[i]);
-
- return sqrt(sum);
- }
-};
typedef PairwiseDistances<PointContainer, L2Distance> PairDistances;
typedef PairDistances::DistanceType DistanceType;
@@ -106,18 +92,8 @@
process_command_line_options(argc, argv, ambient_dimension, skeleton_dimension, multiplier, infilename, outfilename);
// Read in points
- std::ifstream in(infilename.c_str());
PointContainer points;
- while(in)
- {
- points.push_back(Point());
- for (unsigned i = 0; i < ambient_dimension; ++i)
- {
- DistanceType x;
- in >> x;
- points.back().push_back(x);
- }
- }
+ read_points(infilename, points, ambient_dimension);
// Create output file
std::ofstream out(outfilename.c_str());
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/cohomology-persistence.h Mon Mar 23 11:40:29 2009 -0700
@@ -0,0 +1,109 @@
+#ifndef __COHOMOLOGY_PERSISTENCE_H__
+#define __COHOMOLOGY_PERSISTENCE_H__
+
+#if DEBUG_CONTAINERS
+ #include <debug/list>
+ #include <debug/vector>
+ namespace s = std::__debug;
+ #warning "Using debug/list and debug/vector in CohomologyPersistence"
+#else
+ #include <list>
+ #include <vector>
+ namespace s = std;
+#endif
+
+#include <vector>
+#include <list>
+#include <utility>
+
+#include "utilities/types.h"
+
+#include <boost/optional.hpp>
+#include <boost/intrusive/list.hpp>
+namespace bi = boost::intrusive;
+
+
+template<class BirthInfo_, class SimplexData_ = Empty<> >
+class CohomologyPersistence
+{
+ public:
+ typedef BirthInfo_ BirthInfo;
+ typedef SimplexData_ SimplexData;
+
+
+ struct SNode;
+ typedef bi::list<SNode, bi::constant_time_size<false> > ZRow;
+
+ // Simplex representation
+ struct SHead: public SimplexData
+ {
+ SHead(const SHead& other):
+ SimplexData(other), order(other.order) {} // don't copy row since we can't
+ SHead(const SimplexData& sd, unsigned o):
+ SimplexData(sd), order(o) {}
+
+ // intrusive list corresponding to row of s in Z^*, not ordered in any particular order
+ ZRow row;
+ unsigned order;
+ };
+
+ typedef s::list<SHead> Simplices;
+ typedef typename Simplices::iterator SimplexIndex;
+
+ struct Cocycle;
+ typedef s::list<Cocycle> Cocycles;
+ typedef typename Cocycles::iterator CocycleIndex;
+
+ // An entry in a cocycle column; it's also an element in an intrusive list, hence the list_base_hook<>
+ typedef bi::list_base_hook<bi::link_mode<bi::auto_unlink> > auto_unlink_hook;
+ struct SNode: public auto_unlink_hook
+ {
+ SNode() {}
+ SNode(SimplexIndex sidx): si(sidx) {}
+
+ // eventually store a field element
+
+ SimplexIndex si;
+ CocycleIndex cocycle; // TODO: is there no way to get rid of this overhead?
+
+ void unlink() { auto_unlink_hook::unlink(); }
+ };
+ class CompareSNode;
+
+ typedef s::vector<SNode> ZColumn;
+ struct Cocycle
+ {
+ Cocycle(const BirthInfo& b, unsigned o):
+ birth(b), order(o) {}
+
+ ZColumn cocycle;
+ BirthInfo birth;
+ unsigned order;
+
+ bool operator<(const Cocycle& other) const { return order > other.order; }
+ bool operator==(const Cocycle& other) const { return order == other.order; }
+ };
+
+ typedef boost::optional<BirthInfo> Death;
+ typedef std::pair<SimplexIndex,
+ Death> IndexDeathPair;
+
+ // return either a SimplexIndex or a Death
+ // BI = BoundaryIterator; it should dereference to a SimplexIndex
+ template<class BI>
+ IndexDeathPair add(BI begin, BI end, BirthInfo b, const SimplexData& sd = SimplexData());
+
+ void show_cycles() const;
+
+
+ private:
+ void add_cocycle(Cocycle& z1, Cocycle& z2);
+
+ private:
+ Simplices simplices_;
+ Cocycles cocycles_;
+};
+
+#include "cohomology-persistence.hpp"
+
+#endif // __COHOMOLOGY_PERSISTENCE_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/cohomology-persistence.hpp Mon Mar 23 11:40:29 2009 -0700
@@ -0,0 +1,146 @@
+#include <boost/utility.hpp>
+#include <queue>
+#include <vector>
+
+#include <utilities/log.h>
+#include <utilities/indirect.h>
+
+#ifdef LOGGING
+static rlog::RLogChannel* rlCohomology = DEF_CHANNEL("topology/cohomology", rlog::Log_Debug);
+#endif
+
+template<class BirthInfo, class SimplexData>
+class CohomologyPersistence<BirthInfo, SimplexData>::CompareSNode
+{
+ public:
+ bool operator()(const SNode& s1, const SNode& s2) const { return s1.si->order < s2.si->order; }
+};
+
+template<class BirthInfo, class SimplexData>
+template<class BI>
+typename CohomologyPersistence<BirthInfo, SimplexData>::IndexDeathPair
+CohomologyPersistence<BirthInfo, SimplexData>::
+add(BI begin, BI end, BirthInfo birth, const SimplexData& sd)
+{
+ // Create simplex representation
+ simplices_.push_back(SHead(sd, simplices_.empty() ? 0 : (simplices_.back().order + 1)));
+ SimplexIndex si = boost::prior(simplices_.end());
+
+ // Find out if there are cocycles that evaluate to non-zero on the new simplex
+ typedef std::list<CocycleIndex> Candidates;
+ Candidates candidates, candidates_bulk;
+ rLog(rlCohomology, "Boundary");
+ for (BI cur = begin; cur != end; ++cur)
+ {
+ rLog(rlCohomology, " %d", (*cur)->order);
+ for (typename ZRow::const_iterator zcur = (*cur)->row.begin(); zcur != (*cur)->row.end(); ++zcur)
+ candidates_bulk.push_back(zcur->cocycle);
+ }
+
+ candidates_bulk.sort(make_indirect_comparison(std::less<Cocycle>()));
+
+ rLog(rlCohomology, " Candidates bulk");
+ for (typename Candidates::iterator cur = candidates_bulk.begin();
+ cur != candidates_bulk.end(); ++cur)
+ rLog(rlCohomology, " %d", (*cur)->order);
+
+ // Remove duplicates --- this is really Z_2, we need a more sophisticated
+ {
+ typename Candidates::const_iterator cur = candidates_bulk.begin();
+ while (cur != candidates_bulk.end())
+ {
+ typename Candidates::const_iterator next = cur;
+ unsigned count = 0;
+ while (next != candidates_bulk.end() && *next == *cur) { ++next; ++count; }
+
+ if (count % 2)
+ candidates.push_back(*cur);
+
+ cur = next;
+ }
+ }
+
+ // Birth
+ if (candidates.empty())
+ {
+ rLog(rlCohomology, "Birth");
+
+ unsigned order = cocycles_.empty() ? 0 : cocycles_.front().order + 1;
+ cocycles_.push_front(Cocycle(birth, order));
+
+ // set up the cocycle
+ ZColumn& cocycle = cocycles_.front().cocycle;
+ cocycle.push_back(si);
+ cocycle.front().cocycle = cocycles_.begin();
+ si->row.push_back(cocycles_.front().cocycle.front());
+
+ return std::make_pair(si, Death());
+ }
+
+ // Death
+ rLog(rlCohomology, "Death");
+
+#if 0
+ // Debug only, output candidates
+ rLog(rlCohomology, " Candidates");
+ for (typename Candidates::iterator cur = candidates.begin();
+ cur != candidates.end(); ++cur)
+ rLog(rlCohomology, " %d", (*cur)->order);
+#endif
+
+ Cocycle& z = *candidates.front();
+ Death d = z.birth;
+
+ // add z to everything else in candidates
+ for (typename Candidates::iterator cur = boost::next(candidates.begin());
+ cur != candidates.end(); ++cur)
+ add_cocycle(**cur, z);
+
+ for (typename ZColumn::iterator cur = z.cocycle.begin(); cur != z.cocycle.end(); ++cur)
+ cur->unlink();
+
+ cocycles_.erase(candidates.front());
+
+ return std::make_pair(si, d);
+}
+
+template<class BirthInfo, class SimplexData>
+void
+CohomologyPersistence<BirthInfo, SimplexData>::
+show_cycles() const
+{
+ std::cout << "Cocycles" << std::endl;
+ for (typename Cocycles::const_iterator cur = cocycles_.begin(); cur != cocycles_.end(); ++cur)
+ {
+ std::cout << cur->order << ": ";
+ for (typename ZColumn::const_iterator zcur = cur->cocycle.begin(); zcur != cur->cocycle.end(); ++zcur)
+ std::cout << zcur->si->order << ", ";
+ std::cout << std::endl;
+ }
+}
+
+template<class BirthInfo, class SimplexData>
+void
+CohomologyPersistence<BirthInfo, SimplexData>::
+add_cocycle(Cocycle& to, Cocycle& from)
+{
+ rLog(rlCohomology, "Adding cocycle %d to %d", from.order, to.order);
+
+ ZColumn nw;
+
+ std::set_symmetric_difference(to.cocycle.begin(), to.cocycle.end(), // this is catered to Z_2
+ from.cocycle.begin(), from.cocycle.end(),
+ std::back_insert_iterator<ZColumn>(nw),
+ CompareSNode());
+
+ for (typename ZColumn::iterator cur = to.cocycle.begin(); cur != to.cocycle.end(); ++cur)
+ cur->unlink();
+
+ to.cocycle.swap(nw);
+
+ for (typename ZColumn::iterator cur = to.cocycle.begin(); cur != to.cocycle.end(); ++cur)
+ {
+ cur->si->row.push_back(*cur);
+ cur->cocycle = nw[0].cocycle;
+ }
+}