Added rips-zigzag; in the process caught a number of bugs in ZigzagPersistence (added check_consistency() to it)
--- a/examples/rips/CMakeLists.txt Wed Dec 31 11:54:34 2008 -0800
+++ b/examples/rips/CMakeLists.txt Fri Jan 02 14:54:15 2009 -0800
@@ -1,5 +1,6 @@
set (targets
- rips)
+ rips
+ rips-zigzag)
foreach (t ${targets})
add_executable (${t} ${t}.cpp)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/rips/rips-zigzag.cpp Fri Jan 02 14:54:15 2009 -0800
@@ -0,0 +1,157 @@
+#include <topology/rips.h>
+#include <topology/zigzag-persistence.h>
+#include <utilities/types.h>
+#include <utilities/containers.h>
+
+#include <utilities/log.h>
+
+#include <map>
+
+// FIXME
+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 10; }
+ IndexType begin() const { return 0; }
+ IndexType end() const { return size(); }
+};
+
+
+typedef Distances::IndexType Vertex;
+typedef Simplex<Vertex> Smplx;
+typedef ZigzagPersistence<unsigned> Zigzag;
+typedef Zigzag::SimplexIndex Index;
+typedef std::map<Smplx, Index, Smplx::VertexComparison> Complex;
+typedef Zigzag::ZColumn Boundary;
+
+typedef RipsBase<Distances, Smplx> RipsHelper;
+typedef RipsHelper::Evaluator SimplexEvaluator;
+
+
+void make_boundary(const Smplx& s, Complex& c, const Zigzag& zz, Boundary& b)
+{
+ for (Smplx::BoundaryIterator cur = s.boundary_begin(); cur != s.boundary_end(); ++cur)
+ b.append(c[*cur], zz.cmp);
+
+ rDebug(" Boundary:");
+ for (Boundary::const_iterator cur = b.begin(); cur != b.end(); ++cur)
+ rDebug(" %d", (*cur)->order);
+}
+
+int main(int argc, char* argv[])
+{
+#ifdef LOGGING
+ rlog::RLogInit(argc, argv);
+
+ stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+ stderrLog.subscribeTo( RLOG_CHANNEL("info") );
+ stderrLog.subscribeTo( RLOG_CHANNEL("debug") );
+ //stderrLog.subscribeTo( RLOG_CHANNEL("topology/persistence") );
+#endif
+
+ Distances distances;
+
+ // Order vertices and epsilons
+ typedef std::vector<Vertex> VertexVector;
+ typedef std::vector<Distances::DistanceType> EpsilonVector;
+
+ VertexVector vertices;
+ EpsilonVector epsilons;
+ EpsilonVector closest(distances.size(), Infinity);
+
+ vertices.push_back(0);
+ while (vertices.size() < distances.size())
+ {
+ for (Distances::IndexType v = distances.begin(); v != distances.end(); ++v)
+ closest[v] = std::min(closest[v], distances(v, vertices.back()));
+ EpsilonVector::const_iterator max = std::max_element(closest.begin(), closest.end());
+ vertices.push_back(max - closest.begin());
+ epsilons.push_back(*max);
+ }
+
+ std::cout << "Point and epsilon ordering:" << std::endl;
+ for (unsigned i = 0; i < vertices.size(); ++i)
+ std::cout << " " << vertices[i] << " " << epsilons[i] << std::endl;
+
+
+ // Construct zigzag
+ Complex complex;
+ Zigzag zz;
+ RipsHelper aux(distances);
+ SimplexEvaluator size(distances);
+
+ rInfo("Commencing computation");
+ for (unsigned i = 0; i != vertices.size(); ++i)
+ {
+ rInfo("Current stage %d: %d %f", i, vertices[i], epsilons[i]);
+
+ // 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(), i).first));
+ if (!zz.check_consistency())
+ {
+ //zz.show_all();
+ rError("Zigzag representation must be consistent after adding a vertex");
+ }
+ for (Complex::iterator si = complex.begin(); si != complex.end(); ++si)
+ {
+ if (si->first.contains(sv)) continue;
+ rInfo(" Distance between %s and %s: %f",
+ tostring(si->first).c_str(),
+ tostring(sv).c_str(),
+ aux.distance(si->first, sv));
+ if (aux.distance(si->first, sv) <= epsilons[i-1])
+ {
+ Boundary b;
+ Smplx s(si->first); s.join(sv);
+
+ //zz.show_all();
+ 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, i);
+ 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));
+
+ // Death
+ if (idp.second) std::cout << *(idp.second) << " " << i << std::endl;
+ }
+ }
+ rDebug("Complex after addition:");
+ for (Complex::const_iterator cur = complex.begin(); cur != complex.end(); ++cur)
+ rDebug(" %s", tostring(cur->first).c_str());
+
+ rInfo("Removing simplices");
+ // Shrink epsilon
+ {
+ Complex::reverse_iterator si = complex.rbegin();
+
+ while(si != complex.rend())
+ {
+ rInfo(" Size of %s is %f", tostring(si->first).c_str(), size(si->first));
+ if (size(si->first) > epsilons[i])
+ {
+ //zz.show_all();
+ rInfo(" Removing: %s", tostring(si->first).c_str());
+ Zigzag::Death d = zz.remove(si->second, i);
+ AssertMsg(zz.check_consistency(), "Zigzag representation must be consistent after removing a simplex");
+ if (d) std::cout << *d << " " << i << std::endl;
+ complex.erase(boost::prior(si.base()));
+ } else
+ ++si;
+ }
+ rDebug("Complex after removal:");
+ for (Complex::const_iterator cur = complex.begin(); cur != complex.end(); ++cur)
+ rDebug(" %s", tostring(cur->first).c_str());
+ }
+ }
+}
--- a/include/topology/chain.hpp Wed Dec 31 11:54:34 2008 -0800
+++ b/include/topology/chain.hpp Fri Jan 02 14:54:15 2009 -0800
@@ -125,7 +125,6 @@
if (cur != begin()) str += ", ";
str += outmap(*cur);
}
- // str += "(last: " + *last + ")"; // For debugging only
return str;
}
--- a/include/topology/zigzag-persistence.h Wed Dec 31 11:54:34 2008 -0800
+++ b/include/topology/zigzag-persistence.h Fri Jan 02 14:54:15 2009 -0800
@@ -38,19 +38,18 @@
typedef boost::optional<BirthID> Death;
typedef std::pair<SimplexIndex, Death> IndexDeathPair;
- // TODO: probably should store something to identify the birth to the outside world; this should probably
- // be a template parameter (perhaps a template parameter to ZigzagPersistence)
struct ZNode
{
- ZNode(int o, const BirthID& b):
- order(o), birth(b) {}
+ ZNode(int o, const BirthID& b, BIndex l):
+ order(o), birth(b), low(l) {}
int order;
ZColumn z_column;
BRow b_row;
BIndex low; // which BColumn has this ZIndex as low
- BirthID birth; // TODO: need to do empty-member optimization
+ BirthID birth; // TODO: do we need to do empty-member optimization?
+ // i.e., does it ever make sense for birth to be empty?
};
struct BNode
@@ -64,7 +63,8 @@
struct SimplexNode
{
- SimplexNode(unsigned o): order(o) {}
+ SimplexNode(unsigned o, ZIndex l):
+ order(o), low(l) {}
unsigned order;
ZRow z_row;
@@ -81,8 +81,9 @@
// Function: remove(s)
Death remove(SimplexIndex s, const BirthID& birth = BirthID());
- // Debug
+ // Debug; non-const because Indices are iterators, and not const_iterators
void show_all();
+ bool check_consistency();
private:
ZList z_list;
--- a/include/topology/zigzag-persistence.hpp Wed Dec 31 11:54:34 2008 -0800
+++ b/include/topology/zigzag-persistence.hpp Fri Jan 02 14:54:15 2009 -0800
@@ -20,11 +20,12 @@
{ // 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));
+ s_list.push_back(SimplexNode(order, z_list.end()));
}
SimplexIndex last_s = boost::prior(s_list.end());
last_s->low = z_list.end();
+ rLog(rlZigzagAdd, " Reducing among cycles");
// Reduce bdry among the cycles
BColumn v; // representation of the boundary in the cycle basis
while (!bdry.empty())
@@ -33,7 +34,9 @@
ZIndex k = l->low;
v.append(k, cmp);
bdry.add(k->z_column, cmp);
+ rLog(rlZigzagAdd, " Boundary: %s", bdry.tostring(out).c_str());
}
+ rLog(rlZigzagAdd, " Reduced among cycles");
// Reduce v among boundaries
BRow u;
@@ -48,12 +51,15 @@
u.append(k, cmp);
v.add(k->b_column, cmp);
}
+ rLog(rlZigzagAdd, " Reduced among boundaries");
if (v.empty())
{
+ rLog(rlZigzagAdd, " Birth");
+
// Birth
int order = z_list.empty() ? 0 : boost::prior(z_list.end())->order + 1;
- z_list.push_back(ZNode(order, birth));
+ z_list.push_back(ZNode(order, birth, b_list.end()));
ZIndex last_z = boost::prior(z_list.end());
// Set z_column
@@ -71,6 +77,8 @@
return std::make_pair(last_s, Death());
} else
{
+ rLog(rlZigzagAdd, " Death");
+
// Death
unsigned order = b_list.empty() ? 0 : boost::prior(b_list.end())->order + 1;
b_list.push_back(BNode(order));
@@ -102,7 +110,7 @@
remove(SimplexIndex s, const BirthID& birth)
{
rLog(rlZigzagRemove, "Entered ZigzagPersistence::remove(%d)", s->order);
- rLog(rlZigzagRemove, " s->z_row: %s", s->z_row.tostring(out).c_str());
+ AssertMsg(check_consistency(), "Must be consistent before removal");
if (s->z_row.empty())
{
@@ -110,9 +118,10 @@
// Birth
int order = z_list.empty() ? 0 : z_list.begin()->order - 1;
- z_list.push_front(ZNode(order, birth));
+ 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
BIndex j = s->c_row.front();
@@ -133,15 +142,18 @@
add_chains(l_row.begin(), l_row.end(), j, &BNode::b_column, &ZNode::b_row);
// Drop j, l, and s
- // TODO: add some assertions, like everything being empty
//
// l->z_column is the only non-empty thing, but we drop it,
// the basis is preserved because we added first_z
l->z_column.back()->low = z_list.end();
+ std::for_each(l->z_column.begin(), l->z_column.end(), make_remover(&SimplexNode::z_row, l));
+ AssertMsg(l->b_row.empty(), "b_row of l must be empty before erasing in add()");
+ AssertMsg(s->z_row.empty(), "z_row of s must be empty before erasing in add()");
+ AssertMsg(s->c_row.empty(), "c_row of s must be empty before erasing in add()");
b_list.erase(j);
- std::for_each(l->z_column.begin(), l->z_column.end(), make_remover(&SimplexNode::z_row, l));
z_list.erase(l);
s_list.erase(s);
+ AssertMsg(check_consistency(), "Must be consistent when done in add()");
// Reduce Z
SimplexIndex ls = first_z->z_column.back();
@@ -152,8 +164,10 @@
// if ls->low precedes first_z, swap them
if (cmp(ls->low, first_z)) std::swap(ls->low, first_z);
- add_chain(first_z, ls->low, &ZNode::z_column, &SimplexNode::z_row);
+ add_chain(ls->low, first_z, &ZNode::z_column, &SimplexNode::z_row);
std::swap(ls->low, first_z);
+
+ ls = first_z->z_column.back();
}
return Death();
@@ -188,16 +202,22 @@
// Add each reducer to the columns that follow them until the next reducer
for (typename ReducersContainer::reverse_iterator cur = boost::next(reducers.rbegin()); cur != reducers.rend(); ++cur)
+ {
change_basis(*boost::prior(cur), *cur, **cur,
&ZNode::z_column, &SimplexNode::z_row,
&ZNode::b_row, &BNode::b_column);
+ (**cur)->z_column.back()->low = **boost::prior(cur);
+ }
// Drop j and s
- // TODO: add some assertions
BirthID birth = j->birth;
std::for_each(j->z_column.begin(), j->z_column.end(), make_remover(&SimplexNode::z_row, j));
+ AssertMsg(j->b_row.empty(), "b_row of j must be empty before erasing in remove()");
+ 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()");
z_list.erase(j);
s_list.erase(s);
+ AssertMsg(check_consistency(), "Must be consistent when done in remove()");
return Death(birth);
}
@@ -222,6 +242,13 @@
for (typename CRow::const_iterator ccur = cur->c_row.begin(); ccur != cur->c_row.end(); ++ccur)
std::cout << (*ccur)->order << " ";
std::cout << std::endl;
+
+ std::cout << " low: ";
+ if (cur->low != z_list.end())
+ std::cout << cur->low->order;
+ else
+ std::cout << "none";
+ std::cout << std::endl;
}
std::cout << "z_list:" << std::endl;
@@ -241,6 +268,12 @@
std::cout << (*bcur)->order << " ";
std::cout << std::endl;
+ std::cout << " low: ";
+ if (cur->low != b_list.end())
+ std::cout << cur->low->order;
+ else
+ std::cout << "none";
+ std::cout << std::endl;
}
std::cout << "b_list:" << std::endl;
@@ -260,6 +293,46 @@
}
}
+template<class BID>
+bool
+ZigzagPersistence<BID>::
+check_consistency()
+{
+ for (SimplexIndex cur = s_list.begin(); cur != s_list.end(); ++cur)
+ {
+ for (typename ZRow::const_iterator zcur = cur->z_row.begin(); zcur != cur->z_row.end(); ++zcur)
+ if (std::find((*zcur)->z_column.begin(), (*zcur)->z_column.end(), cur) == (*zcur)->z_column.end())
+ return false;
+ for (typename CRow::const_iterator ccur = cur->c_row.begin(); ccur != cur->c_row.end(); ++ccur)
+ if (std::find((*ccur)->c_column.begin(), (*ccur)->c_column.end(), cur) == (*ccur)->c_column.end())
+ return false;
+ if (cur->low != z_list.end() && cur->low->z_column.back() != cur) return false;
+ }
+
+ for (ZIndex cur = z_list.begin(); cur != z_list.end(); ++cur)
+ {
+ for (typename ZColumn::const_iterator scur = cur->z_column.begin(); scur != cur->z_column.end(); ++scur)
+ if (std::find((*scur)->z_row.begin(), (*scur)->z_row.end(), cur) == (*scur)->z_row.end())
+ return false;
+ for (typename BRow::const_iterator bcur = cur->b_row.begin(); bcur != cur->b_row.end(); ++bcur)
+ if (std::find((*bcur)->b_column.begin(), (*bcur)->b_column.end(), cur) == (*bcur)->b_column.end())
+ return false;
+ if (cur->low != b_list.end() && cur->low->b_column.back() != cur) return false;
+ }
+
+ for (BIndex cur = b_list.begin(); cur != b_list.end(); ++cur)
+ {
+ for (typename BColumn::const_iterator zcur = cur->b_column.begin(); zcur != cur->b_column.end(); ++zcur)
+ if (std::find((*zcur)->b_row.begin(), (*zcur)->b_row.end(), cur) == (*zcur)->b_row.end())
+ return false;
+ for (typename CColumn::const_iterator scur = cur->c_column.begin(); scur != cur->c_column.end(); ++scur)
+ if (std::find((*scur)->c_row.begin(), (*scur)->c_row.end(), cur) == (*scur)->c_row.end())
+ return false;
+ }
+
+ return true;
+}
+
/* Private */
// Class: Appender