Implemented ImageZigzagPersistence
* Changed ZigzagPersistence to support a visitor, and implemented
ImageZigzagPersistence
* examples/rips/rips-zigzag now computes using ImageZigzagPersistence
* PersistenceDiagram no longer records zero persistence pairs
* Added utilities/memory.h with report_memory() function
--- a/examples/rips/rips-zigzag.cpp Mon Jan 12 15:33:04 2009 -0800
+++ b/examples/rips/rips-zigzag.cpp Tue Jan 20 10:53:35 2009 -0800
@@ -1,14 +1,17 @@
#include <topology/rips.h>
-#include <topology/zigzag-persistence.h>
+#include <topology/image-zigzag-persistence.h>
#include <utilities/types.h>
#include <utilities/containers.h>
#include <utilities/log.h>
+#include <utilities/memory.h>
#include <map>
#include <cmath>
#include <fstream>
+#include <stack>
+#include <boost/tuple/tuple.hpp>
#include <boost/program_options.hpp>
@@ -43,8 +46,9 @@
typedef RipsHelper::Evaluator SimplexEvaluator;
typedef std::pair<DistanceType, Dimension> BirthInfo;
-typedef ZigzagPersistence<BirthInfo> Zigzag;
+typedef ImageZigzagPersistence<BirthInfo> Zigzag;
typedef Zigzag::SimplexIndex Index;
+typedef Zigzag::Death Death;
typedef std::map<Smplx, Index,
Smplx::VertexDimensionComparison> Complex;
typedef Zigzag::ZColumn Boundary;
@@ -52,15 +56,36 @@
void make_boundary(const Smplx& s, Complex& c, const Zigzag& zz, Boundary& b)
{
- Dimension bdry_dim = s.dimension() - 1;
+ rDebug(" Boundary of <%s>", tostring(s).c_str());
for (Smplx::BoundaryIterator cur = s.boundary_begin(); cur != s.boundary_end(); ++cur)
+ {
b.append(c[*cur], zz.cmp);
+ rDebug(" %d (inL=%d)", c[*cur]->order, b.back()->subcomplex);
+ }
+}
- rDebug(" Boundary:");
- for (Boundary::const_iterator cur = b.begin(); cur != b.end(); ++cur)
- rDebug(" %d", (*cur)->order);
+bool face_leaving_subcomplex(Complex::reverse_iterator si, const SimplexEvaluator& size, DistanceType after, DistanceType before)
+{
+ const Smplx& s = si->first;
+ for (Smplx::VertexContainer::const_iterator v1 = s.vertices().begin(); v1 != s.vertices().end(); ++v1)
+ for (Smplx::VertexContainer::const_iterator v2 = boost::next(v1); v2 != s.vertices().end(); ++v2)
+ {
+ Smplx e; e.add(*v1); e.add(*v2);
+ if (size(e) > after && size(e) <= before)
+ return true;
+ }
+
+ return false;
}
+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;
+}
+
+
std::ostream& operator<<(std::ostream& out, const BirthInfo& bi)
{ return (out << bi.first); }
@@ -75,12 +100,12 @@
stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
#endif
- SetFrequency(cOperations, 500);
+ SetFrequency(cOperations, 25000);
SetTrigger(cOperations, cComplexSize);
unsigned ambient_dimension;
unsigned skeleton_dimension;
- float multiplier;
+ float from_multiplier, to_multiplier;
std::string infilename;
po::options_description hidden("Hidden options");
@@ -92,7 +117,8 @@
("help,h", "produce help message")
("ambient-dimsnion,a", po::value<unsigned>(&ambient_dimension)->default_value(3), "The ambient dimension of the point set")
("skeleton-dimsnion,s", po::value<unsigned>(&skeleton_dimension)->default_value(2), "Dimension of the Rips complex we want to compute")
- ("multiplier,m", po::value<float>(&multiplier)->default_value(4), "Multiplier for the epsilon (distance to next maxmin point) when computing the Rips complex");
+ ("from,f", po::value<float>(&from_multiplier)->default_value(4), "From multiplier for the epsilon (distance to next maxmin point) when computing the Rips complex")
+ ("to,t", po::value<float>(&to_multiplier)->default_value(16), "To multiplier for the epsilon (distance to next maxmin point) when computing the Rips complex");
#if LOGGING
std::vector<std::string> log_channels;
visible.add_options()
@@ -152,6 +178,7 @@
EpsilonVector dist(distances.size(), Infinity);
vertices.push_back(distances.begin());
+ //epsilons.push_back(Infinity);
while (vertices.size() < distances.size())
{
for (Vertex v = distances.begin(); v != distances.end(); ++v)
@@ -160,6 +187,7 @@
vertices.push_back(max - dist.begin());
epsilons.push_back(*max);
}
+ epsilons.push_back(0);
}
rInfo("Point and epsilon ordering:");
@@ -173,6 +201,10 @@
RipsHelper aux(distances);
SimplexEvaluator size(distances);
+ // TODO: it probably makes sense to do things in reverse.
+ // I.e., we should start from the smallest epsilon, and grow, rather than
+ // starting from the largest epsilon and shrinking since the interesting
+ // part of the computation is that with small epsilon.
rInfo("Commencing computation");
for (unsigned i = 0; i != vertices.size(); ++i)
{
@@ -181,7 +213,10 @@
// Add a point
Smplx sv; sv.add(vertices[i]);
rDebug("Added %s", tostring(sv).c_str());
- complex.insert(std::make_pair(sv, zz.add(Boundary(), std::make_pair(epsilons[i], 0)).first));
+ complex.insert(std::make_pair(sv,
+ zz.add(Boundary(),
+ true, // vertex is always in the subcomplex
+ std::make_pair(epsilons[i], 0)).first));
CountNum(cComplexSize, 0);
Count(cComplexSize);
Count(cOperations);
@@ -199,7 +234,7 @@
tostring(si->first).c_str(),
tostring(sv).c_str(),
aux.distance(si->first, sv));
- if (aux.distance(si->first, sv) <= multiplier*epsilons[i-1])
+ if (aux.distance(si->first, sv) <= to_multiplier*epsilons[i-1])
{
Boundary b;
Smplx s(si->first); s.join(sv);
@@ -208,51 +243,93 @@
rDebug("Adding %s", tostring(s).c_str());
make_boundary(s, complex, zz, b);
rDebug("Made boundary, %d", b.size());
- Zigzag::IndexDeathPair idp = zz.add(b, std::make_pair(epsilons[i], sv.dimension()));
+ 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()));
if (!zz.check_consistency())
{
//zz.show_all();
rError("Zigzag representation must be consistent after adding a simplex");
}
- complex.insert(std::make_pair(s, idp.first));
+ complex.insert(std::make_pair(s, idx));
CountNum(cComplexSize, s.dimension());
Count(cComplexSize);
Count(cOperations);
// Death
- if (idp.second) std::cout << (idp.second)->second << " " << (idp.second)->first << " " << epsilons[i] << std::endl;
+ if (d && ((d->first - epsilons[i-1]) != 0) && (d->second < skeleton_dimension))
+ std::cout << d->second << " " << d->first << " " << epsilons[i-1] << std::endl;
}
}
rDebug("Complex after addition:");
for (Complex::const_iterator si = complex.begin(); si != complex.end(); ++si)
- rDebug(" %s", tostring(si->first).c_str());
+ rDebug(" %s", tostring(si->first).c_str());
+ rInfo("Inserted point; complex size: %d", complex.size());
+ show_image_betti(zz, skeleton_dimension);
+ report_memory();
+
+ if (i == 0) continue; // want to skip the removal from the image check (involving epsilons[i-1]),
+ // and in any case, there is only one vertex at this point
rDebug("Removing simplices");
// Shrink epsilon
{
+ std::stack<Complex::reverse_iterator> leaving_subcomplex;
Complex::reverse_iterator si = complex.rbegin();
while(si != complex.rend())
{
rDebug(" Size of %s is %f", tostring(si->first).c_str(), size(si->first));
- if (size(si->first) > multiplier*epsilons[i])
+ if (size(si->first) > to_multiplier*epsilons[i])
{
//zz.show_all();
- rDebug(" Removing: %s", tostring(si->first).c_str());
- Zigzag::Death d = zz.remove(si->second,
- std::make_pair(epsilons[i], si->first.dimension() - 1));
+ 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));
AssertMsg(zz.check_consistency(), "Zigzag representation must be consistent after removing a simplex");
- if (d) std::cout << d->second << " " << d->first << " " << epsilons[i] << std::endl;
+ if (d && ((d->first - epsilons[i]) != 0) && (d->second < skeleton_dimension))
+ std::cout << d->second << " " << d->first << " " << epsilons[i] << std::endl;
CountNumBy(cComplexSize, si->first.dimension(), -1);
complex.erase(boost::prior(si.base()));
CountBy(cComplexSize, -1);
Count(cOperations);
+ } else if (face_leaving_subcomplex(si, size, from_multiplier*epsilons[i], from_multiplier*epsilons[i-1]))
+ {
+ // 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));
+ 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;
+ leaving_subcomplex.push(si++);
} else
++si;
}
+ while(!leaving_subcomplex.empty())
+ {
+ si = leaving_subcomplex.top(); // copying an iterator onto stack is probably Ok
+ Boundary b;
+ make_boundary(si->first, complex, zz, b);
+ 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()));
+ 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;
+ leaving_subcomplex.pop();
+ }
}
rDebug("Complex after removal:");
for (Complex::const_iterator si = complex.begin(); si != complex.end(); ++si)
rDebug(" %s", tostring(si->first).c_str());
+
+ rInfo("Shrunk epsilon; complex size: %d", complex.size());
+ show_image_betti(zz, skeleton_dimension);
+ report_memory();
}
}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/image-zigzag-persistence.h Tue Jan 20 10:53:35 2009 -0800
@@ -0,0 +1,106 @@
+#ifndef __IMAGE_ZIGZAG_PERSISTENCE_H__
+#define __IMAGE_ZIGZAG_PERSISTENCE_H__
+
+#include "zigzag-persistence.h"
+#include <limits>
+
+struct SimplexSubcomplexData
+{
+ SimplexSubcomplexData(bool sc = false):
+ subcomplex(sc) {}
+
+ bool subcomplex;
+};
+
+template<class BirthID_ = Empty<> >
+class ImageZigzagPersistence: public ZigzagPersistence<BirthID_, SimplexSubcomplexData>
+{
+ public:
+ typedef BirthID_ BirthID;
+ typedef ZigzagPersistence<BirthID, SimplexSubcomplexData> Parent;
+
+ typedef typename Parent::IndexDeathPair IndexDeathPair;
+ typedef typename Parent::Death Death;
+
+ typedef typename Parent::ZIndex ZIndex;
+ typedef typename Parent::BIndex BIndex;
+ typedef typename Parent::SimplexIndex SimplexIndex;
+ typedef typename Parent::ZColumn ZColumn;
+ typedef typename Parent::BColumn BColumn;
+ typedef typename Parent::BRow BRow;
+ typedef typename Parent::CRow CRow;
+ typedef typename Parent::ZNode ZNode;
+ typedef typename Parent::BNode BNode;
+ typedef typename Parent::SimplexNode SimplexNode;
+
+
+ ImageZigzagPersistence():
+ im_last(Parent::z_list.end()), cok_begin(Parent::z_list.end()),
+ im_order_begin(std::numeric_limits<int>::min()/2),
+ cok_order_begin(std::numeric_limits<int>::max()/2)
+ {}
+
+ IndexDeathPair add(ZColumn bdry,
+ bool subcomplex,
+ const BirthID& birth = BirthID()) { ImageZZVisitor zzv(subcomplex); return Parent::add(bdry, birth, zzv); }
+
+ Death remove(SimplexIndex s, const BirthID& birth = BirthID()) { ImageZZVisitor zzv(s->subcomplex); return Parent::remove(s, birth, zzv); }
+
+
+ ZIndex image_begin() { return Parent::z_list.begin(); }
+ ZIndex image_end() { return cok_begin; }
+ BIndex boundary_end() { return Parent::b_list.end(); }
+
+
+ // Class: ImageZZVisitor
+ // Contains all the tweaks to the normal zigzag algorithm to make it compute image zigzag
+ class ImageZZVisitor: public Parent::ZigzagVisitor
+ {
+ public:
+ ImageZZVisitor(bool sc = false):
+ subcomplex(sc), birth_in_image(false) {}
+
+ // Sets the subcomplex property of the new simplex
+ SimplexIndex new_simplex(Parent& zz);
+
+ // Decides where to put the new column (image or cokernel)
+ ZIndex new_z_in_add(Parent& zz, const ZColumn& z, const BRow& u);
+
+ // Checks if there is a boundary entirely in the subcomplex, and sets birth_in_image accordingly
+ BIndex select_j_in_remove(Parent& zz, const CRow& c_row);
+
+ ZIndex new_z_in_remove(Parent& zz);
+
+ // Updates im_last and cok_begin if necessary
+ void erasing_z(Parent& zz, ZIndex j);
+
+ // Determines if there is a death in the image
+ Death death(Parent& zz, ZIndex dying_z);
+
+ private:
+ ZIndex append_in_image(Parent& zz);
+ ZIndex append_in_cokernel(Parent& zz);
+ ZIndex prepend_in_image(Parent& zz);
+ ZIndex prepend_in_cokernel(Parent& zz);
+ bool in_subcomplex(const ZColumn& z);
+
+ bool subcomplex, birth_in_image;
+ };
+
+ const int im_order_begin;
+ const int cok_order_begin;
+
+ private:
+ using Parent::make_remover;
+ using Parent::make_appender;
+ using Parent::make_adder;
+
+ private:
+ ZIndex im_last; // index of the last image cycle
+ ZIndex cok_begin; // index of the first cokernel cycle
+};
+
+
+#include "image-zigzag-persistence.hpp"
+
+#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/image-zigzag-persistence.hpp Tue Jan 20 10:53:35 2009 -0800
@@ -0,0 +1,213 @@
+#ifdef LOGGING
+static rlog::RLogChannel* rlImageZigzag = DEF_CHANNEL("topology/persistence/zigzag/image", rlog::Log_Debug);
+#endif // LOGGING
+
+
+template<class BID>
+typename ImageZigzagPersistence<BID>::SimplexIndex
+ImageZigzagPersistence<BID>::ImageZZVisitor::
+new_simplex(Parent& zz)
+{
+ SimplexIndex si = Parent::ZigzagVisitor::new_simplex(zz);
+ si->subcomplex = subcomplex;
+ rLog(rlImageZigzag, "New simplex %d (inL=%d)", si->order, si->subcomplex);
+ return si;
+}
+
+template<class BID>
+typename ImageZigzagPersistence<BID>::ZIndex
+ImageZigzagPersistence<BID>::ImageZZVisitor::
+new_z_in_add(Parent& zz, const ZColumn& z, const BRow& u)
+{
+ ImageZigzagPersistence& izz = static_cast<ImageZigzagPersistence&>(zz);
+
+ // Check if z is entirely in the subcomplex
+ if (in_subcomplex(z))
+ return append_in_image(zz);
+ else
+ {
+ append_in_cokernel(zz);
+
+ // Simplex we are adding is in the subcomplex
+ if (subcomplex)
+ {
+ rLog(rlImageZigzag, "Modifying boundaries");
+
+ SimplexIndex s = boost::prior(izz.s_list.end());
+ AssertMsg(s->subcomplex, "The new simplex must be in the subcomplex");
+ rLog(rlImageZigzag, " s=%d", s->order);
+ rLog(rlImageZigzag, " u=[%s]", u.tostring(izz.out).c_str());
+
+ BIndex max = u.front();
+ for (typename BRow::const_iterator cur = u.begin(); cur != u.end(); ++cur)
+ if ((*cur)->b_column.back()->order > max->b_column.back()->order)
+ max = *cur;
+
+ // Replace B[max] with sum of b_columns from u
+ BColumn sum;
+ std::for_each(u.begin(), u.end(), izz.make_adder(&BNode::b_column, sum));
+#if ZIGZAG_CONSISTENCY
+ AssertMsg(in_subcomplex(s->boundary), "Boundary of s must be in the subcomplex");
+#endif
+ std::for_each(max->b_column.begin(), max->b_column.end(), izz.make_remover(&ZNode::b_row, max));
+ rLog(rlImageZigzag, "max->b_column=[%s]", max->b_column.tostring(izz.out).c_str());
+ rLog(rlImageZigzag, " sum=[%s]", sum.tostring(izz.out).c_str());
+ AssertMsg(sum.back() == max->b_column.back(), "Must be replacing by a column with the same low");
+ max->b_column.swap(sum);
+ std::for_each(max->b_column.begin(), max->b_column.end(), izz.make_appender(&ZNode::b_row, max));
+ // NB: low doesn't need to be adjusted (see AssertMsg above)
+
+ // Wipe out C[max], and replace it with s
+ std::for_each(max->c_column.begin(), max->c_column.end(), izz.make_remover(&SimplexNode::c_row, max));
+ max->c_column.clear();
+ max->c_column.append(s, izz.cmp);
+ AssertMsg(s->c_row.empty(), "s was just added, so it cannot appear in any bounding chain");
+ s->c_row.append(max, izz.cmp);
+ }
+
+ return boost::prior(izz.z_list.end());
+ }
+}
+
+
+template<class BID>
+typename ImageZigzagPersistence<BID>::BIndex
+ImageZigzagPersistence<BID>::ImageZZVisitor::
+select_j_in_remove(Parent& zz, const CRow& c_row)
+{
+ ImageZigzagPersistence& izz = static_cast<ImageZigzagPersistence&>(zz);
+ for (typename CRow::const_iterator cur = c_row.begin(); cur != c_row.end(); ++cur)
+ if ((*cur)->b_column.back()->order <= izz.im_last->order)
+ {
+ birth_in_image = true;
+ return *cur;
+ }
+
+ return Parent::ZigzagVisitor::select_j_in_remove(zz, c_row);
+}
+
+template<class BID>
+typename ImageZigzagPersistence<BID>::ZIndex
+ImageZigzagPersistence<BID>::ImageZZVisitor::
+new_z_in_remove(Parent& zz)
+{
+ if (birth_in_image)
+ return prepend_in_image(zz);
+ else
+ return prepend_in_cokernel(zz);
+}
+
+template<class BID>
+void
+ImageZigzagPersistence<BID>::ImageZZVisitor::
+erasing_z(Parent& zz, ZIndex j)
+{
+ ImageZigzagPersistence& izz = static_cast<ImageZigzagPersistence&>(zz);
+ if (j == izz.im_last) --(izz.im_last);
+ else if (j == izz.cok_begin) ++(izz.cok_begin);
+}
+
+template<class BID>
+typename ImageZigzagPersistence<BID>::Death
+ImageZigzagPersistence<BID>::ImageZZVisitor::
+death(Parent& zz, ZIndex dying_z)
+{
+ ImageZigzagPersistence& izz = static_cast<ImageZigzagPersistence&>(zz);
+ if (izz.im_last == izz.z_list.end() || dying_z->order > izz.im_last->order)
+ return Death();
+ else
+ return Death(dying_z->birth);
+}
+
+template<class BID>
+typename ImageZigzagPersistence<BID>::ZIndex
+ImageZigzagPersistence<BID>::ImageZZVisitor::
+append_in_image(Parent& zz)
+{
+ rLog(rlImageZigzag, "Appending in image");
+ ImageZigzagPersistence& izz = static_cast<ImageZigzagPersistence&>(zz);
+
+ // if no cycles in the image
+ if (izz.im_last == izz.z_list.end())
+ {
+ izz.z_list.push_front(ZNode(izz.im_order_begin, izz.b_list.end()));
+ return (izz.im_last = izz.z_list.begin());
+ } else
+ {
+ izz.z_list.insert(boost::next(izz.im_last), ZNode(izz.im_last->order + 1, izz.b_list.end()));
+ return ++(izz.im_last);
+ }
+}
+
+template<class BID>
+typename ImageZigzagPersistence<BID>::ZIndex
+ImageZigzagPersistence<BID>::ImageZZVisitor::
+append_in_cokernel(Parent& zz)
+{
+ rLog(rlImageZigzag, "Appending in cokernel");
+ ImageZigzagPersistence& izz = static_cast<ImageZigzagPersistence&>(zz);
+
+ // if no cycles in the cokernel
+ if (izz.cok_begin == izz.z_list.end())
+ {
+ izz.z_list.push_back(ZNode(izz.cok_order_begin, izz.b_list.end()));
+ izz.cok_begin = boost::prior(izz.z_list.end());
+ } else
+ {
+ izz.z_list.push_back(ZNode(boost::prior(izz.z_list.end())->order + 1, izz.b_list.end()));
+ }
+
+ return boost::prior(izz.z_list.end());
+}
+
+template<class BID>
+typename ImageZigzagPersistence<BID>::ZIndex
+ImageZigzagPersistence<BID>::ImageZZVisitor::
+prepend_in_image(Parent& zz)
+{
+ rLog(rlImageZigzag, "Prepending in image");
+ ImageZigzagPersistence& izz = static_cast<ImageZigzagPersistence&>(zz);
+
+ // if no cycles in the image
+ if (izz.im_last == izz.z_list.end())
+ {
+ izz.z_list.push_front(ZNode(izz.im_order_begin, izz.b_list.end()));
+ return (izz.im_last = izz.z_list.begin());
+ } else
+ {
+ izz.z_list.push_front(ZNode(izz.z_list.begin()->order - 1, izz.b_list.end()));
+ return izz.z_list.begin();
+ }
+}
+
+template<class BID>
+typename ImageZigzagPersistence<BID>::ZIndex
+ImageZigzagPersistence<BID>::ImageZZVisitor::
+prepend_in_cokernel(Parent& zz)
+{
+ rLog(rlImageZigzag, "Prepending in cokernel");
+ ImageZigzagPersistence& izz = static_cast<ImageZigzagPersistence&>(zz);
+
+ // if no cycles in the cokernel
+ if (izz.cok_begin == izz.z_list.end())
+ {
+ izz.z_list.push_back(ZNode(izz.cok_order_begin, izz.b_list.end()));
+ izz.cok_begin = boost::prior(izz.z_list.end());
+ } else
+ {
+ izz.z_list.insert(izz.cok_begin, ZNode(izz.cok_begin->order - 1, izz.b_list.end()));
+ --(izz.cok_begin);
+ }
+ return izz.cok_begin;
+}
+
+template<class BID>
+bool
+ImageZigzagPersistence<BID>::ImageZZVisitor::
+in_subcomplex(const ZColumn& z)
+{
+ for (typename ZColumn::const_iterator cur = z.begin(); cur != z.end(); ++cur)
+ if (!((*cur)->subcomplex))
+ return false;
+ return true;
+}
--- a/include/topology/persistence-diagram.h Mon Jan 12 15:33:04 2009 -0800
+++ b/include/topology/persistence-diagram.h Tue Jan 20 10:53:35 2009 -0800
@@ -8,6 +8,7 @@
#include <cmath>
#include <boost/compressed_pair.hpp>
+#include <boost/optional.hpp>
#include <boost/serialization/access.hpp>
/**
@@ -58,10 +59,13 @@
{ return (point.operator<<(out)); }
template<class Point, class Iterator, class Evaluator, class Visitor>
-Point make_point(Iterator i, const Evaluator& evaluator, const Visitor& visitor);
+boost::optional<Point>
+make_point(Iterator i, const Evaluator& evaluator, const Visitor& visitor);
template<class Point, class Iterator, class Evaluator>
-Point make_point(Iterator i, const Evaluator& evaluator) { return make_point<Point>(i, evaluator, Point::Visitor()); }
+boost::optional<Point>
+make_point(Iterator i, const Evaluator& evaluator)
+{ return make_point<Point>(i, evaluator, Point::Visitor()); }
/**
--- a/include/topology/persistence-diagram.hpp Mon Jan 12 15:33:04 2009 -0800
+++ b/include/topology/persistence-diagram.hpp Tue Jan 20 10:53:35 2009 -0800
@@ -48,11 +48,14 @@
{
for (Iterator cur = bg; cur != end; ++cur)
if (cur->sign())
- push_back(make_point(cur, evaluator, visitor));
+ {
+ boost::optional<Point> p = make_point(cur, evaluator, visitor);
+ if (p) push_back(*p);
+ }
}
template<class Point, class Iterator, class Evaluator, class Visitor>
-Point
+boost::optional<Point>
make_point(Iterator i, const Evaluator& evaluator, const Visitor& visitor)
{
RealType x = evaluator(i);
@@ -62,6 +65,8 @@
Point p(x,y);
visitor.point(i, p);
+
+ if (x == y) return boost::optional<Point>();
return p;
}
@@ -90,7 +95,11 @@
for (Iterator cur = bg; cur != end; ++cur)
if (cur->sign())
- diagrams[dimension(cur)].push_back(make_point<typename PDiagram::Point>(cur, evaluator, visitor));
+ {
+ boost::optional<typename PDiagram::Point> p = make_point<typename PDiagram::Point>(cur, evaluator, visitor);
+ if (p)
+ diagrams[dimension(cur)].push_back(*p);
+ }
}
template<class D>
--- a/include/topology/zigzag-persistence.h Mon Jan 12 15:33:04 2009 -0800
+++ b/include/topology/zigzag-persistence.h Tue Jan 20 10:53:35 2009 -0800
@@ -19,11 +19,12 @@
* Class: ZigzagPersistence
* TODO: this should probably be parametrized by Chain or Field
*/
-template<class BirthID_ = Empty<> >
+template<class BirthID_ = Empty<>, class SimplexData_ = Empty<> >
class ZigzagPersistence
{
public:
typedef BirthID_ BirthID;
+ typedef SimplexData_ SimplexData;
struct ZNode;
struct BNode;
@@ -49,8 +50,8 @@
struct ZNode
{
- ZNode(int o, const BirthID& b, BIndex l):
- order(o), birth(b), low(l) {}
+ ZNode(int o, BIndex l):
+ order(o), low(l) {}
int order;
ZColumn z_column;
@@ -70,7 +71,7 @@
CColumn c_column;
};
- struct SimplexNode
+ struct SimplexNode: public SimplexData
{
SimplexNode(unsigned o, ZIndex l):
order(o), low(l) {}
@@ -79,26 +80,60 @@
ZRow z_row;
CRow c_row;
ZIndex low; // which ZColumn has this SimplexNode as low
-#if !NDEBUG
+#if ZIGZAG_CONSISTENCY
ZColumn boundary; // NB: debug only
#endif
};
// Constructor: ZigzagPersistence()
ZigzagPersistence() {}
+
+ // Function: add(bdry, birth)
+ IndexDeathPair add(ZColumn bdry, const BirthID& birth = BirthID()) { ZigzagVisitor zzv; return add<ZigzagVisitor>(bdry, birth, zzv); }
+
+ // Function: remove(s, birth)
+ Death remove(SimplexIndex s, const BirthID& birth = BirthID()) { ZigzagVisitor zzv; return remove<ZigzagVisitor>(s, birth, zzv); }
+
+ protected:
// Function: add(s)
- IndexDeathPair add(ZColumn bdry, const BirthID& birth = BirthID());
+ template<class Visitor>
+ IndexDeathPair add(ZColumn bdry, const BirthID& birth, Visitor& visitor);
// Function: remove(s)
- Death remove(SimplexIndex s, const BirthID& birth = BirthID());
+ template<class Visitor>
+ Death remove(SimplexIndex s, const BirthID& birth, Visitor& visitor);
+
+ // Struct: ZigzagVisitor
+ // Various methods of an instance of this class are called at different stages of addition and removal algorithm.
+ // NB: currently the places where it's called are catered for image zigzags, in the future this could be expanded
+ // to provide simple support for other algorithms
+ // TODO: not obvious that the methods should be const (and therefore the reference passed to add() and remove())
+ // revisit when working on ImageZigzag
+ struct ZigzagVisitor
+ {
+ SimplexIndex new_simplex(ZigzagPersistence& zz);
+ // Function: new_z_in_add(zz, z, u)
+ // Called when a new cycle is born after adding a simplex. The method is expected to add an element to z_list, and return its ZIndex.
+ ZIndex new_z_in_add(ZigzagPersistence& zz, const ZColumn& z, const BRow& u);
+
+ BIndex select_j_in_remove(ZigzagPersistence& zz, const CRow& c_row);
+
+ ZIndex new_z_in_remove(ZigzagPersistence& zz);
+
+ void erasing_z(ZigzagPersistence& zz, ZIndex j) {}
+
+ Death death(ZigzagPersistence& zz, ZIndex dying_z);
+ };
+
+ public:
// Debug; non-const because Indices are iterators, and not const_iterators
void show_all();
bool check_consistency(SimplexIndex s_skip, ZIndex z_skip, BIndex b_skip);
bool check_consistency() { return check_consistency(s_list.end(), z_list.end(), b_list.end()); }
- private:
+ protected:
ZList z_list;
BList b_list;
SimplexList s_list;
--- a/include/topology/zigzag-persistence.hpp Mon Jan 12 15:33:04 2009 -0800
+++ b/include/topology/zigzag-persistence.hpp Tue Jan 20 10:53:35 2009 -0800
@@ -12,23 +12,20 @@
#endif // LOGGING
-template<class BID>
-typename ZigzagPersistence<BID>::IndexDeathPair
-ZigzagPersistence<BID>::
-add(ZColumn bdry, const BirthID& birth)
+template<class BID, class SD>
+template<class Visitor>
+typename ZigzagPersistence<BID,SD>::IndexDeathPair
+ZigzagPersistence<BID,SD>::
+add(ZColumn bdry, const BirthID& birth, Visitor& visitor)
{
rLog(rlZigzagAdd, "Entered ZigzagPersistence::add()");
rLog(rlZigzagAdd, " Boundary: %s", bdry.tostring(out).c_str());
rLog(rlZigzagAdd, " Boundary size: %d", bdry.size());
AssertMsg(check_consistency(), "Must be consistent before addition");
- { // scoping to not pollute with the name order
- unsigned order = s_list.empty() ? 0 : boost::prior(s_list.end())->order + 1;
- s_list.push_back(SimplexNode(order, z_list.end()));
- }
- SimplexIndex last_s = boost::prior(s_list.end());
+ SimplexIndex last_s = visitor.new_simplex(*this);
last_s->low = z_list.end();
-#if !NDEBUG
+#if ZIGZAG_CONSISTENCY
last_s->boundary = bdry; // NB: debug only
#endif
@@ -65,22 +62,24 @@
{
rLog(rlZigzagAdd, " Birth case in add");
- // Birth
- int order = z_list.empty() ? 0 : boost::prior(z_list.end())->order + 1;
- z_list.push_back(ZNode(order, birth, b_list.end()));
- ZIndex last_z = boost::prior(z_list.end());
-
- // Set z_column
- ZColumn& z = last_z->z_column;
+ // Figure out the new cycle z
+ ZColumn z;
std::for_each(u.begin(), u.end(), make_adder(&BNode::c_column, z));
z.append(last_s, cmp);
+ // Birth
+ ZIndex new_z = visitor.new_z_in_add(*this, z, u);
+ new_z->birth = birth;
+
// Set s_row
- std::for_each(z.begin(), z.end(), make_appender(&SimplexNode::z_row, last_z));
+ std::for_each(z.begin(), z.end(), make_appender(&SimplexNode::z_row, new_z));
+
+ // Set z_column
+ new_z->z_column.swap(z);
// Set low
- last_z->low = b_list.end();
- last_s->low = last_z;
+ new_z->low = b_list.end();
+ last_s->low = new_z;
return std::make_pair(last_s, Death());
} else
@@ -107,15 +106,16 @@
// Set c_row
std::for_each(c.begin(), c.end(), make_appender(&SimplexNode::c_row, last_b));
- return std::make_pair(last_s, Death(last_b->b_column.back()->birth));
+ return std::make_pair(last_s, visitor.death(*this, last_b->b_column.back()));
}
}
-template<class BID>
-typename ZigzagPersistence<BID>::Death
-ZigzagPersistence<BID>::
-remove(SimplexIndex s, const BirthID& birth)
+template<class BID, class SD>
+template<class Visitor>
+typename ZigzagPersistence<BID,SD>::Death
+ZigzagPersistence<BID,SD>::
+remove(SimplexIndex s, const BirthID& birth, Visitor& visitor)
{
rLog(rlZigzagRemove, "Entered ZigzagPersistence::remove(%d)", s->order);
AssertMsg(check_consistency(), "Must be consistent before removal");
@@ -128,18 +128,20 @@
//show_all();
rLog(rlZigzagRemove, "Birth case in remove");
- int order = z_list.empty() ? 0 : z_list.begin()->order - 1;
- z_list.push_front(ZNode(order, birth, b_list.end()));
- ZIndex first_z = z_list.begin();
- ZColumn& z = first_z->z_column;
- first_z->low = b_list.end();
-
// Prepend DC[j] = ZB[j] to Z
rLog(rlZigzagRemove, "Computing the column DC[j] = ZB[j] to prepend to Z");
- BIndex j = s->c_row.front();
+ BIndex j = visitor.select_j_in_remove(*this, s->c_row);
rLog(rlZigzagRemove, " j = %d", j->order);
+
+ ZColumn z;
std::for_each(j->b_column.begin(), j->b_column.end(), make_adder(&ZNode::z_column, z));
+
+ ZIndex first_z = visitor.new_z_in_remove(*this);
+ first_z->birth = birth;
std::for_each(z.begin(), z.end(), make_appender(&SimplexNode::z_row, first_z));
+ first_z->z_column.swap(z);
+ first_z->low = b_list.end();
+
rLog(rlZigzagRemove, " Prepended %d [%s]", first_z->order, z.tostring(out).c_str());
//AssertMsg(check_consistency(), "Must be consistent after prepending DC[j] = ZB[j] to Z");
@@ -150,7 +152,7 @@
std::for_each(first_z->b_row.begin(), first_z->b_row.end(), make_appender(&BNode::b_column, first_z));
//AssertMsg(check_consistency(), "Must be consistent after prepending row of s to B");
-#if !NDEBUG
+#if ZIGZAG_CONSISTENCY
{
ZColumn zz;
std::for_each(j->b_column.begin(), j->b_column.end(), make_adder(&ZNode::z_column, zz));
@@ -158,12 +160,18 @@
}
#endif
+ typedef std::not_equal_to<BIndex> NotEqualBIndex;
+
// Subtract C[j] from every column of C that contains s
AssertMsg(s->c_row == first_z->b_row, "s->c_row == first_z->b_row before subtracting C[j]");
rLog(rlZigzagRemove, "Subtracting C[j]=[%s] from every column of C that contains s=%d with row [%s]",
j->c_column.tostring(out).c_str(),
s->order, s->c_row.tostring(out).c_str());
- add_chains(first_z->b_row.rbegin(), first_z->b_row.rend(), j, &BNode::c_column, &SimplexNode::c_row);
+ add_chains(boost::make_filter_iterator(std::bind2nd(NotEqualBIndex(), j), first_z->b_row.begin(), first_z->b_row.end()),
+ boost::make_filter_iterator(std::bind2nd(NotEqualBIndex(), j), first_z->b_row.end(), first_z->b_row.end()),
+ j, &BNode::c_column, &SimplexNode::c_row);
+ add_chain(j, j, &BNode::c_column, &SimplexNode::c_row);
+ // TODO: remove add_chains(first_z->b_row.rbegin(), first_z->b_row.rend(), j, &BNode::c_column, &SimplexNode::c_row);
AssertMsg(check_consistency(s_list.end(), z_list.begin(), b_list.end()), "Must be consistent after subtracting C[j] in remove::birth");
// Subtract B[j] from every other column of B that has l
@@ -171,7 +179,6 @@
BRow l_row = l->b_row;
rLog(rlZigzagRemove, "Subtracting B[j], j is %d, l is %d, l_row: [%s]",
j->order, l->order, l_row.tostring(out).c_str());
- typedef std::not_equal_to<BIndex> NotEqualBIndex;
add_chains(boost::make_filter_iterator(std::bind2nd(NotEqualBIndex(), j), l_row.rbegin(), l_row.rend()),
boost::make_filter_iterator(std::bind2nd(NotEqualBIndex(), j), l_row.rend(), l_row.rend()),
j, &BNode::b_column, &ZNode::b_row);
@@ -192,7 +199,14 @@
AssertMsg(l->b_row.empty(), "b_row of l must be empty before erasing in remove::birth");
AssertMsg(s->z_row.empty(), "z_row of s must be empty before erasing in remove::birth");
+ rLog(rlZigzagRemove, "s->c_row: [%s]", s->c_row.tostring(out).c_str());
+ if (!s->c_row.empty())
+ {
+ rLog(rlZigzagRemove, "s->c_row[0]: [%s]", s->c_row.front()->c_column.tostring(out).c_str());
+ rLog(rlZigzagRemove, " b_column: [%s]", s->c_row.front()->b_column.tostring(out).c_str());
+ }
AssertMsg(s->c_row.empty(), "c_row of s must be empty before erasing in remove::birth");
+ visitor.erasing_z(*this, l);
b_list.erase(j);
z_list.erase(l);
s_list.erase(s);
@@ -224,6 +238,10 @@
ZIndex j = s->z_row.front();
CRow c_row = s->c_row;
+
+ rLog(rlZigzagRemove, "j=%d, j->b_row=[%s]", j->order, j->b_row.tostring(out).c_str());
+ rLog(rlZigzagRemove, "s=%d, s->c_row=[%s]", s->order, s->c_row.tostring(out).c_str());
+ rLog(rlZigzagRemove, "s=%d, s->z_row=[%s]", s->order, s->z_row.tostring(out).c_str());
// Subtract Z[j] from every chain in C that contains s
// (it's Ok to go in the forward order since we are subtracting a column in Z from C)
@@ -285,14 +303,21 @@
}
// Drop j and s
- BirthID birth = j->birth;
+ Death d = visitor.death(*this, j);
+
if (j->z_column.back()->low == j)
j->z_column.back()->low = z_list.end();
std::for_each(j->z_column.begin(), j->z_column.end(), make_remover(&SimplexNode::z_row, j));
rLog(rlZigzagRemove, "j->b_row: [%s]", j->b_row.tostring(out).c_str());
- AssertMsg(j->b_row.empty(), "b_row of j must be empty before erasing in remove()");
+ if (!j->b_row.empty())
+ {
+ rLog(rlZigzagRemove, "j->b_row[0]: [%s]", j->b_row.front()->b_column.tostring(out).c_str());
+ rLog(rlZigzagRemove, " c_column: [%s]", j->b_row.front()->c_column.tostring(out).c_str());
+ }
+ AssertMsg(j->b_row.empty(), "b_row of j must be empty before erasing in remove(). Most likely you are trying to remove a simplex whose coface is still in the complex.");
AssertMsg(s->z_row.empty(), "z_row of s must be empty before erasing in remove()");
AssertMsg(s->c_row.empty(), "c_row of s must be empty before erasing in remove()");
+ visitor.erasing_z(*this, j);
z_list.erase(j);
s_list.erase(s);
@@ -300,13 +325,13 @@
AssertMsg(check_consistency(), "Must be consistent when done in remove()");
- return Death(birth);
+ return d;
}
}
-template<class BID>
+template<class BID, class SD>
void
-ZigzagPersistence<BID>::
+ZigzagPersistence<BID,SD>::
show_all()
{
std::cout << "s_list:" << std::endl;
@@ -374,12 +399,12 @@
}
}
-template<class BID>
+template<class BID, class SD>
bool
-ZigzagPersistence<BID>::
+ZigzagPersistence<BID,SD>::
check_consistency(SimplexIndex s_skip, ZIndex z_skip, BIndex b_skip)
{
-#if !NDEBUG
+#ifdef ZIGZAG_CONSISTENCY
for (SimplexIndex cur = s_list.begin(); cur != s_list.end(); ++cur)
{
if (cur == s_skip) continue;
@@ -485,9 +510,9 @@
// Class: Appender
//
// Functor that appends given element to the given member of whatever parameter it is invoked with
-template<class BID>
+template<class BID, class SD>
template<class Member, class Element>
-struct ZigzagPersistence<BID>::Appender
+struct ZigzagPersistence<BID,SD>::Appender
{
Appender(Member mm, Element ee):
m(mm), e(ee) {}
@@ -503,9 +528,9 @@
// Class: Remover
//
// Functor that removes given element from the given member of whatever parameter it is invoked with
-template<class BID>
+template<class BID, class SD>
template<class Member, class Element>
-struct ZigzagPersistence<BID>::Remover
+struct ZigzagPersistence<BID,SD>::Remover
{
Remover(Member mm, Element ee):
m(mm), e(ee) {}
@@ -520,9 +545,9 @@
// Class: Adder
//
// Functor that adds the given member of whatever it is invoked with to the given chain
-template<class BID>
+template<class BID, class SD>
template<class Member, class Chain>
-struct ZigzagPersistence<BID>::Adder
+struct ZigzagPersistence<BID,SD>::Adder
{
Adder(Member mm, Chain& cc):
m(mm), c(cc) {}
@@ -540,10 +565,10 @@
//
// Special case of add_chains where all Indexes are the same, and
// therefore PrimaryMemberFrom and PrimaryMemberTo are the same
-template<class BID>
+template<class BID, class SD>
template<class Index, class IndexFrom, class PrimaryMember, class SecondaryMember>
void
-ZigzagPersistence<BID>::
+ZigzagPersistence<BID,SD>::
add_chains(Index bg, Index end, IndexFrom j, PrimaryMember pm, SecondaryMember sm)
{
add_chains(bg, end, j, pm, sm, pm);
@@ -555,10 +580,10 @@
// Fixes SecondaryMembers by adding and removing the corresponding elements.
// For example, if we add a column to a number of other columns, then PrimaryMember is that
// column member, and SecondaryMember is the corresponding row member.
-template<class BID>
+template<class BID, class SD>
template<class IndexTo, class IndexFrom, class PrimaryMemberTo, class SecondaryMemberTo, class PrimaryMemberFrom>
void
-ZigzagPersistence<BID>::
+ZigzagPersistence<BID,SD>::
add_chains(IndexTo bg, IndexTo end, IndexFrom j, PrimaryMemberTo pmt, SecondaryMemberTo smt, PrimaryMemberFrom pmf)
{
for (IndexTo cur = bg; cur != end; ++cur)
@@ -570,10 +595,10 @@
// Changes basis by adding PrimaryMember pm of j to pm of every element in range [bg, end).
// In parallel it performs the reverse (complementary) update on the dual members, i.e.
// column and row operations are performed in sync, so that the product of the two matrices doesn't change
-template<class BID>
+template<class BID, class SD>
template<class IndexTo, class IndexFrom, class PrimaryMember, class SecondaryMember, class DualPrimaryMember, class DualSecondaryMember>
void
-ZigzagPersistence<BID>::
+ZigzagPersistence<BID,SD>::
change_basis(IndexTo bg, IndexTo end, IndexFrom j, PrimaryMember pm, SecondaryMember sm, DualPrimaryMember dpm, DualSecondaryMember dsm)
{
for (IndexTo cur = bg; cur != end; ++cur)
@@ -583,10 +608,10 @@
}
}
-template<class BID>
+template<class BID, class SD>
template<class Index, class PrimaryMember, class SecondaryMember>
void
-ZigzagPersistence<BID>::
+ZigzagPersistence<BID,SD>::
add_chain(Index to, Index from, PrimaryMember pm, SecondaryMember sm)
{
add_chain(to, from, pm, sm, pm);
@@ -596,10 +621,10 @@
//
// Adds PrimaryMemberFrom pmf of `from` to PrimaryMemberTo pmt of `to`.
// Fixes SecondaryMemberTos. See add_chains().
-template<class BID>
+template<class BID, class SD>
template<class IndexTo, class IndexFrom, class PrimaryMemberTo, class SecondaryMemberTo, class PrimaryMemberFrom>
void
-ZigzagPersistence<BID>::
+ZigzagPersistence<BID,SD>::
add_chain(IndexTo to, IndexFrom from, PrimaryMemberTo pmt, SecondaryMemberTo smt, PrimaryMemberFrom pmf)
{
rLog(rlZigzagAddChain, "Adding %d [%s] to %d [%s]",
@@ -628,3 +653,51 @@
((*to).*pmt).add((*from).*pmf, cmp);
rLog(rlZigzagAddChain, "Got %s", ((*to).*pmt).tostring(out).c_str());
}
+
+
+/* ZigzagVisitor */
+template<class BID, class SD>
+typename ZigzagPersistence<BID,SD>::SimplexIndex
+ZigzagPersistence<BID,SD>::ZigzagVisitor::
+new_simplex(ZigzagPersistence& zz)
+{
+ unsigned order = zz.s_list.empty() ? 0 : boost::prior(zz.s_list.end())->order + 1;
+ zz.s_list.push_back(SimplexNode(order, zz.z_list.end()));
+ return boost::prior(zz.s_list.end());
+}
+
+template<class BID, class SD>
+typename ZigzagPersistence<BID,SD>::ZIndex
+ZigzagPersistence<BID,SD>::ZigzagVisitor::
+new_z_in_add(ZigzagPersistence& zz, const ZColumn& z, const BRow& u)
+{
+ int order = zz.z_list.empty() ? 0 : boost::prior(zz.z_list.end())->order + 1;
+ zz.z_list.push_back(ZNode(order, zz.b_list.end()));
+ return boost::prior(zz.z_list.end());
+}
+
+template<class BID, class SD>
+typename ZigzagPersistence<BID,SD>::BIndex
+ZigzagPersistence<BID,SD>::ZigzagVisitor::
+select_j_in_remove(ZigzagPersistence& zz, const CRow& c_row)
+{
+ return c_row.front();
+}
+
+template<class BID, class SD>
+typename ZigzagPersistence<BID,SD>::ZIndex
+ZigzagPersistence<BID,SD>::ZigzagVisitor::
+new_z_in_remove(ZigzagPersistence& zz)
+{
+ int order = zz.z_list.empty() ? 0 : zz.z_list.begin()->order - 1;
+ zz.z_list.push_front(ZNode(order, zz.b_list.end()));
+ return zz.z_list.begin();
+}
+
+template<class BID, class SD>
+typename ZigzagPersistence<BID,SD>::Death
+ZigzagPersistence<BID,SD>::ZigzagVisitor::
+death(ZigzagPersistence& zz, ZIndex dying_z)
+{
+ return Death(dying_z->birth);
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/utilities/memory.h Tue Jan 20 10:53:35 2009 -0800
@@ -0,0 +1,41 @@
+#ifndef __MEMORY_H__
+#define __MEMORY_H__
+
+#if LOGGING // TODO: add check for Linux (and preferably the right version of Linux)
+
+#include "log.h"
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <unistd.h>
+
+static rlog::RLogChannel* rlMemory = DEF_CHANNEL("memory", rlog::Log_Debug);
+
+void report_memory()
+{
+ pid_t pid = getpid();
+ std::stringstream smaps_name;
+ smaps_name << "/proc/" << pid << "/smaps";
+ std::ifstream in(smaps_name.str().c_str());
+
+ std::string str;
+ unsigned memory = 0, value;
+ while (in)
+ {
+ in >> str;
+ if (std::string(str, 0, 7) == "Private")
+ {
+ in >> value;
+ memory += value;
+ }
+ }
+ rLog(rlMemory, "Private memory usage: %d kB", memory);
+}
+
+#else
+
+void report_memory() {}
+
+#endif // LOGGING
+
+#endif // __MEMORY_H__