Debugged ZigzagPersistence (having added heavier consistency checking)
* Added DEBUG_CONTAINERS option (uses std::__debug::* containers for chains
and in ZigzagPersistence)
* Added SizeStorage specialization for std::deque<T>
* ZigzagPersistence got a lot more consistency checking (in debug mode only,
which now crawls); as a result it's been debugged (running on non-trivial examples)
* examples/rips/rips-zigzag takes command-line options
* added ChainWrapper::clear()
* added Simplex::VertexDimensionComparison
* added PairwiseDistances class (for computing distances between points in a
container according to a distance functor)
--- a/examples/rips/CMakeLists.txt Fri Jan 02 14:54:15 2009 -0800
+++ b/examples/rips/CMakeLists.txt Mon Jan 12 15:33:04 2009 -0800
@@ -4,5 +4,5 @@
foreach (t ${targets})
add_executable (${t} ${t}.cpp)
- target_link_libraries (${t} ${libraries})
+ target_link_libraries (${t} ${libraries} ${Boost_PROGRAM_OPTIONS_LIBRARY})
endforeach (t ${targets})
--- a/examples/rips/rips-zigzag.cpp Fri Jan 02 14:54:15 2009 -0800
+++ b/examples/rips/rips-zigzag.cpp Mon Jan 12 15:33:04 2009 -0800
@@ -6,34 +6,53 @@
#include <utilities/log.h>
#include <map>
-
-// FIXME
-struct Distances
-{
- typedef int IndexType;
- typedef double DistanceType;
+#include <cmath>
+#include <fstream>
- 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(); }
-};
+#include <boost/program_options.hpp>
-typedef Distances::IndexType Vertex;
+#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
+{
+ typedef double value_type;
+
+ value_type operator()(const Point& p1, const Point& p2) const
+ {
+ AssertMsg(p1.size() == p2.size(), "Points must be in the same dimension (in L2Distance");
+ value_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;
+
+typedef PairDistances::IndexType Vertex;
typedef Simplex<Vertex> Smplx;
-typedef ZigzagPersistence<unsigned> Zigzag;
+
+typedef RipsBase<PairDistances, Smplx> RipsHelper;
+typedef RipsHelper::Evaluator SimplexEvaluator;
+
+typedef std::pair<DistanceType, Dimension> BirthInfo;
+typedef ZigzagPersistence<BirthInfo> Zigzag;
typedef Zigzag::SimplexIndex Index;
-typedef std::map<Smplx, Index, Smplx::VertexComparison> Complex;
+typedef std::map<Smplx, Index,
+ Smplx::VertexDimensionComparison> 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)
{
+ Dimension bdry_dim = s.dimension() - 1;
for (Smplx::BoundaryIterator cur = s.boundary_begin(); cur != s.boundary_end(); ++cur)
b.append(c[*cur], zz.cmp);
@@ -42,40 +61,110 @@
rDebug(" %d", (*cur)->order);
}
+std::ostream& operator<<(std::ostream& out, const BirthInfo& bi)
+{ return (out << bi.first); }
+
+namespace po = boost::program_options;
+
+
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") );
+ stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
+#endif
+
+ SetFrequency(cOperations, 500);
+ SetTrigger(cOperations, cComplexSize);
+
+ unsigned ambient_dimension;
+ unsigned skeleton_dimension;
+ float multiplier;
+ std::string infilename;
+
+ 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<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");
+#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
- Distances distances;
+ 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)
+ stdoutLog.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;
+ return 1;
+ }
+
+ // 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);
+ }
+ }
+
+ // Create pairwise distances
+ PairDistances distances(points);
// Order vertices and epsilons
typedef std::vector<Vertex> VertexVector;
- typedef std::vector<Distances::DistanceType> EpsilonVector;
+ typedef std::vector<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);
+ EpsilonVector dist(distances.size(), Infinity);
+
+ vertices.push_back(distances.begin());
+ while (vertices.size() < distances.size())
+ {
+ for (Vertex v = distances.begin(); v != distances.end(); ++v)
+ dist[v] = std::min(dist[v], distances(v, vertices.back()));
+ EpsilonVector::const_iterator max = std::max_element(dist.begin(), dist.end());
+ vertices.push_back(max - dist.begin());
+ epsilons.push_back(*max);
+ }
}
- std::cout << "Point and epsilon ordering:" << std::endl;
+ rInfo("Point and epsilon ordering:");
for (unsigned i = 0; i < vertices.size(); ++i)
- std::cout << " " << vertices[i] << " " << epsilons[i] << std::endl;
+ rInfo(" %d %f", vertices[i], epsilons[i]);
// Construct zigzag
@@ -92,7 +181,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(), i).first));
+ complex.insert(std::make_pair(sv, zz.add(Boundary(), std::make_pair(epsilons[i], 0)).first));
+ CountNum(cComplexSize, 0);
+ Count(cComplexSize);
+ Count(cOperations);
if (!zz.check_consistency())
{
//zz.show_all();
@@ -100,58 +192,67 @@
}
for (Complex::iterator si = complex.begin(); si != complex.end(); ++si)
{
- if (si->first.contains(sv)) continue;
- rInfo(" Distance between %s and %s: %f",
+ if (si->first.contains(sv)) continue;
+ if (si->first.dimension() + 1 > skeleton_dimension) continue;
+
+ rDebug(" 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])
+ if (aux.distance(si->first, sv) <= multiplier*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);
+ Zigzag::IndexDeathPair idp = zz.add(b, std::make_pair(epsilons[i], sv.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));
+ CountNum(cComplexSize, s.dimension());
+ Count(cComplexSize);
+ Count(cOperations);
// Death
- if (idp.second) std::cout << *(idp.second) << " " << i << std::endl;
+ if (idp.second) std::cout << (idp.second)->second << " " << (idp.second)->first << " " << epsilons[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());
+ for (Complex::const_iterator si = complex.begin(); si != complex.end(); ++si)
+ rDebug(" %s", tostring(si->first).c_str());
- rInfo("Removing simplices");
+ rDebug("Removing simplices");
// Shrink epsilon
{
- Complex::reverse_iterator si = complex.rbegin();
+ 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])
+ rDebug(" Size of %s is %f", tostring(si->first).c_str(), size(si->first));
+ if (size(si->first) > multiplier*epsilons[i])
{
//zz.show_all();
- rInfo(" Removing: %s", tostring(si->first).c_str());
- Zigzag::Death d = zz.remove(si->second, i);
+ rDebug(" Removing: %s", tostring(si->first).c_str());
+ Zigzag::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 << " " << i << std::endl;
+ if (d) 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
++si;
}
- rDebug("Complex after removal:");
- for (Complex::const_iterator cur = complex.begin(); cur != complex.end(); ++cur)
- rDebug(" %s", tostring(cur->first).c_str());
}
+ rDebug("Complex after removal:");
+ for (Complex::const_iterator si = complex.begin(); si != complex.end(); ++si)
+ rDebug(" %s", tostring(si->first).c_str());
}
}
--- a/examples/rips/rips.cpp Fri Jan 02 14:54:15 2009 -0800
+++ b/examples/rips/rips.cpp Mon Jan 12 15:33:04 2009 -0800
@@ -1,5 +1,6 @@
#include <topology/rips.h>
+// Trivial example of size() points on a line with integer coordinates
struct Distances
{
typedef int IndexType;
@@ -24,18 +25,19 @@
Distances distances;
#if 0
+ // Storing ExplicitDistances speeds up the computation (at the price of memory)
typedef RipsGenerator<ExplicitDistances<Distances> > RipsGenerator;
ExplicitDistances<Distances> explicit_distances(distances);
RipsGenerator rips(explicit_distances);
#else
- typedef RipsGeneratorMemory<Distances> RipsGenerator;
+ //typedef RipsGeneratorMemory<Distances> RipsGenerator;
+ typedef RipsGenerator<Distances> RipsGenerator;
RipsGenerator rips(distances);
#endif
RipsGenerator::SimplexVector complex;
//rips.generate(complex, 3, distances.size());
rips.generate(complex, 3, 50);
- //rips.print();
std::cout << "Size: " << complex.size() << std::endl;
// for (RipsGenerator::SimplexVector::const_iterator cur = complex.begin(); cur != complex.end(); ++cur)
--- a/include/topology/chain.h Fri Jan 02 14:54:15 2009 -0800
+++ b/include/topology/chain.h Mon Jan 12 15:33:04 2009 -0800
@@ -65,6 +65,7 @@
Self& add(const Self& c, const ConsistencyComparison& cmp);
void swap(ChainWrapper& c); ///< Swaps the contents of c and *this (like STL's swap destroys c)
+ void clear();
template<class ConsistencyComparison>
void sort(const ConsistencyComparison& cmp); ///< Sort elements to enforce consistency
@@ -72,7 +73,6 @@
size_t size() const { return Size::size(*this); }
using ChainRepresentation::empty;
- using ChainRepresentation::clear;
/// @}
/// \name Modifiers
--- a/include/topology/chain.hpp Fri Jan 02 14:54:15 2009 -0800
+++ b/include/topology/chain.hpp Mon Jan 12 15:33:04 2009 -0800
@@ -87,6 +87,15 @@
}
template<class C>
+void
+ChainWrapper<C>::
+clear()
+{
+ ChainRepresentation::clear();
+ Size::clear();
+}
+
+template<class C>
template<class ConsistencyComparison>
void
ChainWrapper<C>::
--- a/include/topology/cycles.h Fri Jan 02 14:54:15 2009 -0800
+++ b/include/topology/cycles.h Mon Jan 12 15:33:04 2009 -0800
@@ -2,15 +2,25 @@
#define __CYCLES_H__
#include "chain.h"
-#include <vector>
-#include <deque>
#include "utilities/circular_list.h"
+#if DEBUG_CONTAINERS
+ #include <debug/vector>
+ #include <debug/deque>
+ using std::__debug::vector;
+ using std::__debug::deque;
+#else
+ #include <vector>
+ #include <deque>
+ using std::vector;
+ using std::deque;
+#endif
+
template<class OrderIndex_ = int>
struct VectorChains
{
typedef OrderIndex_ OrderIndex;
- typedef ChainWrapper<std::vector<OrderIndex> > Chain;
+ typedef ChainWrapper<vector<OrderIndex> > Chain;
typedef Chain Cycle;
Cycle cycle;
@@ -28,7 +38,7 @@
struct DequeChains
{
typedef OrderIndex_ OrderIndex;
- typedef ChainWrapper<std::deque<OrderIndex> > Chain;
+ typedef ChainWrapper<deque<OrderIndex> > Chain;
typedef Chain Cycle;
Cycle cycle;
--- a/include/topology/rips.h Fri Jan 02 14:54:15 2009 -0800
+++ b/include/topology/rips.h Mon Jan 12 15:33:04 2009 -0800
@@ -119,7 +119,8 @@
};
/**
- * ExplicitDistances stores the pairwise distances of Distances_ instance passed at construction.
+ * Class: ExplicitDistances
+ * Stores the pairwise distances of Distances_ instance passed at construction.
* It's a protypical Distances template argument for the Rips complex.
*/
template<class Distances_>
@@ -143,6 +144,39 @@
size_t size_;
};
+
+/**
+ * Class: PairwiseDistances
+ * Given a Container_ of points and a Distance_, it computes distances between elements
+ * in the container (given as instances of Index_ defaulted to unsigned) using the Distance_ functor.
+ *
+ * Container_ is assumed to be an std::vector. That simplifies a number of things.
+ */
+template<class Container_, class Distance_, typename Index_ = unsigned>
+class PairwiseDistances
+{
+ public:
+ typedef Container_ Container;
+ typedef Distance_ Distance;
+ typedef Index_ IndexType;
+ typedef typename Distance::value_type DistanceType;
+
+
+ PairwiseDistances(const Container& container,
+ const Distance& distance = Distance()):
+ container_(container), distance_(distance) {}
+
+ DistanceType operator()(IndexType a, IndexType b) const { return distance_(container_[a], container_[b]); }
+
+ size_t size() const { return container_.size(); }
+ IndexType begin() const { return 0; }
+ IndexType end() const { return size(); }
+
+ private:
+ const Container& container_;
+ Distance distance_;
+};
+
#include "rips.hpp"
#endif // __RIPS_H__
--- a/include/topology/simplex.h Fri Jan 02 14:54:15 2009 -0800
+++ b/include/topology/simplex.h Mon Jan 12 15:33:04 2009 -0800
@@ -110,6 +110,7 @@
/* Classes: Comparisons
*
* VertexComparison - compare simplices based on the lexicographic ordering of their <vertices()>
+ * VertexComparison - compare simplices based on the lexicographic ordering of their <vertices()> within each dimension
* DataComparison - compare simplices based on their <data()>
* DataDimensionComparison - compare simplices based on their <data()> within each <dimension()>
*/
@@ -119,6 +120,7 @@
* DimensionExtractor - return dimesnion given a simplex
*/
struct VertexComparison;
+ struct VertexDimensionComparison;
struct DataComparison;
struct DataDimensionComparison;
@@ -150,6 +152,22 @@
};
template<class V, class T>
+struct Simplex<V,T>::VertexDimensionComparison
+{
+ typedef Self first_argument_type;
+ typedef Self second_argument_type;
+ typedef bool result_type;
+
+ bool operator()(const Self& a, const Self& b) const
+ {
+ if (a.dimension() == b.dimension())
+ return a.vertices() < b.vertices();
+ else
+ return a.dimension() < b.dimension();
+ }
+};
+
+template<class V, class T>
struct Simplex<V,T>::DataComparison
{
typedef Self first_argument_type;
--- a/include/topology/zigzag-persistence.h Fri Jan 02 14:54:15 2009 -0800
+++ b/include/topology/zigzag-persistence.h Mon Jan 12 15:33:04 2009 -0800
@@ -1,11 +1,20 @@
#ifndef __ZIGZAG_PERSISTENCE_H__
#define __ZIGZAG_PERSISTENCE_H__
-#include <list>
#include "cycles.h"
#include "utilities/types.h"
#include <sstream>
+#if DEBUG_CONTAINERS
+ #include <debug/list>
+ using std::__debug::list;
+ #warning "Using debug/list in ZigzagPersistence"
+#else
+ #include <list>
+ using std::list;
+#endif
+
+
/**
* Class: ZigzagPersistence
* TODO: this should probably be parametrized by Chain or Field
@@ -20,11 +29,11 @@
struct BNode;
struct SimplexNode;
- typedef std::list<ZNode> ZList;
+ typedef list<ZNode> ZList;
typedef typename ZList::iterator ZIndex;
- typedef std::list<BNode> BList;
+ typedef list<BNode> BList;
typedef typename BList::iterator BIndex;
- typedef std::list<SimplexNode> SimplexList;
+ typedef list<SimplexNode> SimplexList;
typedef typename SimplexList::iterator SimplexIndex;
// TODO: should all chains be DequeChains? probably not
@@ -70,6 +79,9 @@
ZRow z_row;
CRow c_row;
ZIndex low; // which ZColumn has this SimplexNode as low
+#if !NDEBUG
+ ZColumn boundary; // NB: debug only
+#endif
};
// Constructor: ZigzagPersistence()
@@ -83,7 +95,8 @@
// Debug; non-const because Indices are iterators, and not const_iterators
void show_all();
- bool check_consistency();
+ 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:
ZList z_list;
--- a/include/topology/zigzag-persistence.hpp Fri Jan 02 14:54:15 2009 -0800
+++ b/include/topology/zigzag-persistence.hpp Mon Jan 12 15:33:04 2009 -0800
@@ -2,10 +2,13 @@
#include <boost/utility.hpp>
#include <algorithm>
#include <utilities/indirect.h>
+#include <functional>
#ifdef LOGGING
-static rlog::RLogChannel* rlZigzagAdd = DEF_CHANNEL("topology/persistence/zigzag/add", rlog::Log_Debug);
-static rlog::RLogChannel* rlZigzagRemove = DEF_CHANNEL("topology/persistence/zigzag/remove", rlog::Log_Debug);
+static rlog::RLogChannel* rlZigzagAdd = DEF_CHANNEL("topology/persistence/zigzag/add", rlog::Log_Debug);
+static rlog::RLogChannel* rlZigzagRemove = DEF_CHANNEL("topology/persistence/zigzag/remove", rlog::Log_Debug);
+static rlog::RLogChannel* rlZigzagAddChain = DEF_CHANNEL("topology/persistence/zigzag/addchain", rlog::Log_Debug);
+static rlog::RLogChannel* rlZigzagCheckConsistency= DEF_CHANNEL("topology/persistence/zigzag/check", rlog::Log_Debug);
#endif // LOGGING
@@ -17,6 +20,7 @@
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;
@@ -24,9 +28,13 @@
}
SimplexIndex last_s = boost::prior(s_list.end());
last_s->low = z_list.end();
+#if !NDEBUG
+ last_s->boundary = bdry; // NB: debug only
+#endif
rLog(rlZigzagAdd, " Reducing among cycles");
// Reduce bdry among the cycles
+ rLog(rlZigzagAdd, " Boundary: %s", bdry.tostring(out).c_str());
BColumn v; // representation of the boundary in the cycle basis
while (!bdry.empty())
{
@@ -55,7 +63,7 @@
if (v.empty())
{
- rLog(rlZigzagAdd, " Birth");
+ rLog(rlZigzagAdd, " Birth case in add");
// Birth
int order = z_list.empty() ? 0 : boost::prior(z_list.end())->order + 1;
@@ -77,7 +85,7 @@
return std::make_pair(last_s, Death());
} else
{
- rLog(rlZigzagAdd, " Death");
+ rLog(rlZigzagAdd, " Death case in add");
// Death
unsigned order = b_list.empty() ? 0 : boost::prior(b_list.end())->order + 1;
@@ -117,6 +125,9 @@
AssertMsg(!(s->c_row.empty()), "Birth after removal, row in C must be non-empty");
// Birth
+ //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();
@@ -124,22 +135,50 @@
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();
+ rLog(rlZigzagRemove, " j = %d", j->order);
std::for_each(j->b_column.begin(), j->b_column.end(), make_adder(&ZNode::z_column, z));
std::for_each(z.begin(), z.end(), make_appender(&SimplexNode::z_row, first_z));
+ 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");
// Prepend row of s in C to B
+ rLog(rlZigzagRemove, "Prepending the row of s to B");
first_z->b_row = s->c_row; // copying instead of swapping is inefficient,
// but it simplifies logic when subtracting chains later
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
+ {
+ ZColumn zz;
+ std::for_each(j->b_column.begin(), j->b_column.end(), make_adder(&ZNode::z_column, zz));
+ AssertMsg(zz.empty(), "ZB[j] must be 0 after we prepended the row of s in C to B");
+ }
+#endif
// Subtract C[j] from every column of C that contains s
- add_chains(first_z->b_row.begin(), first_z->b_row.end(), j, &BNode::c_column, &SimplexNode::c_row);
+ 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);
+ 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
ZIndex l = j->b_column.back();
- BRow l_row = l->b_row; // again, copying is inefficient, but it simplifies the logic
- add_chains(l_row.begin(), l_row.end(), j, &BNode::b_column, &ZNode::b_row);
+ 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);
+ j->b_column.back()->low = b_list.end(); // redundant since l will be deleted (here for check_consistency only)
+ add_chain(j, j, &BNode::b_column, &ZNode::b_row);
+ AssertMsg(check_consistency(s_list.end(), first_z, b_list.end()), "Must be consistent after subtracting B[j] in remove::birth");
+
// Drop j, l, and s
//
@@ -147,15 +186,20 @@
// 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()");
+
+ //show_all();
+ rLog(rlZigzagRemove, "l=%d has z_column: [%s]", l->order, l->z_column.tostring(out).c_str());
+
+ 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");
+ AssertMsg(s->c_row.empty(), "c_row of s must be empty before erasing in remove::birth");
b_list.erase(j);
z_list.erase(l);
s_list.erase(s);
- AssertMsg(check_consistency(), "Must be consistent when done in add()");
+ AssertMsg(check_consistency(s_list.end(), first_z, b_list.end()), "Must be consistent before reducing Z in remove::birth");
// Reduce Z
+ rLog(rlZigzagRemove, "Reducing Z");
SimplexIndex ls = first_z->z_column.back();
while(ls->low != first_z)
{
@@ -164,59 +208,96 @@
// 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::b_row, &BNode::b_column);
add_chain(ls->low, first_z, &ZNode::z_column, &SimplexNode::z_row);
std::swap(ls->low, first_z);
ls = first_z->z_column.back();
}
+ AssertMsg(check_consistency(), "Must be consistent at the end of birth case in remove");
return Death();
} else
{
// Death
+ rLog(rlZigzagRemove, "Death case in remove");
+
ZIndex j = s->z_row.front();
CRow c_row = s->c_row;
// 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)
add_chains(c_row.begin(), c_row.end(), j, &BNode::c_column, &SimplexNode::c_row, &ZNode::z_column);
+ AssertMsg(check_consistency(), "Must be consistent after subtracting Z[j] from C");
- /**
- * Change basis to remove s from Z
- * All the processing is done from back to front, so that when elements are removed from s->z_row,
- * it doesn't invalidate the rest of the elements (iterators) in the row.
- */
+ // Change basis to remove s from Z
// Compute reducers --- columns that we will be adding to other columns
+ ZRow z_row = s->z_row;
typedef typename ZRow::reverse_iterator ZRowReverseIterator;
typedef std::list<ZRowReverseIterator> ReducersContainer;
ReducersContainer reducers; // list of ZColumns that we will be adding to other columns
- reducers.push_back(boost::prior(s->z_row.rend())); // j is the first reducer
+ reducers.push_back(boost::prior(z_row.rend())); // j is the first reducer
+ AssertMsg(*(reducers.back()) == j, "The first element of reducers should be j");
SimplexIndex low = j->z_column.back();
- for (typename ZRow::iterator cur = s->z_row.begin(); cur != s->z_row.end(); ++cur)
+ rLog(rlZigzagRemove, " Added reducer %d [%s] with low=%d",
+ j->order, j->z_column.tostring(out).c_str(),
+ low->order);
+ for (typename ZRow::iterator cur = z_row.begin(); cur != z_row.end(); ++cur)
if (cmp((*cur)->z_column.back(), low))
{
- reducers.push_back(ZRowReverseIterator(cur));
+ reducers.push_back(ZRowReverseIterator(boost::next(cur)));
low = (*(reducers.back()))->z_column.back();
+ rLog(rlZigzagRemove, " Added reducer %d [%s] with low=%d",
+ (*cur)->order, (*cur)->z_column.tostring(out).c_str(),
+ low->order);
+ rLog(rlZigzagRemove, " reducers.back(): %d [%s] with low=%d",
+ (*(reducers.back()))->order,
+ (*(reducers.back()))->z_column.tostring(out).c_str(),
+ (*(reducers.back()))->z_column.back()->order);
}
- reducers.push_back(s->z_row.rbegin());
- //rLog(rlZigzagRemove, " Reducers size: %d", reducers.size());
+ rLog(rlZigzagRemove, " Reducers size: %d, s is %d",
+ reducers.size(), s->order);
+
+ //show_all();
// 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)
+ typename ReducersContainer::reverse_iterator cur = reducers.rbegin();
+ ZRowReverseIterator zcur = z_row.rbegin();
+
+ while (cur != reducers.rend())
{
- change_basis(*boost::prior(cur), *cur, **cur,
+ rLog(rlZigzagRemove, " Cur reducer: %d [%s]", (**cur)->order,
+ (**cur)->z_column.tostring(out).c_str());
+ change_basis(zcur, *cur, **cur,
&ZNode::z_column, &SimplexNode::z_row,
&ZNode::b_row, &BNode::b_column);
- (**cur)->z_column.back()->low = **boost::prior(cur);
+ if (cur != reducers.rbegin())
+ {
+ AssertMsg((*zcur)->z_column.back() == (**cur)->z_column.back(),
+ "The back of the z_columns must be the same.");
+ (*zcur)->z_column.back()->low = *zcur;
+ }
+ else
+ (**cur)->z_column.back()->low = z_list.end();
+
+ zcur = *cur++;
+ // This makes it inconsistent until the next iteration of this update loop
}
// Drop j and s
BirthID birth = j->birth;
+ 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()");
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);
+
+ //show_all();
+
AssertMsg(check_consistency(), "Must be consistent when done in remove()");
return Death(birth);
@@ -296,39 +377,105 @@
template<class BID>
bool
ZigzagPersistence<BID>::
-check_consistency()
+check_consistency(SimplexIndex s_skip, ZIndex z_skip, BIndex b_skip)
{
+#if !NDEBUG
for (SimplexIndex cur = s_list.begin(); cur != s_list.end(); ++cur)
{
+ if (cur == s_skip) continue;
+ //rLog(rlZigzagCheckConsistency, "SimplexIndex cur: %d", cur->order);
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())
+ {
+ rError("In check_consistency(): SimplexNode %d not found in z_column of %d", cur->order, (*zcur)->order);
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())
+ {
+ rError("In check_consistency(): SimplexNode %d not found in c_column of %d", cur->order, (*ccur)->order);
return false;
- if (cur->low != z_list.end() && cur->low->z_column.back() != cur) return false;
+ }
+ if (cur->low != z_list.end())
+ AssertMsg(!(cur->low->z_column.empty()), "z_column must not be empty");
+ if (cur->low != z_list.end() && cur->low->z_column.back() != cur)
+ {
+ rError("low of SimplexNode %d is incorrect", cur->order);
+ return false;
+ }
}
for (ZIndex cur = z_list.begin(); cur != z_list.end(); ++cur)
{
+ if (cur == z_skip) continue;
+
+ //rLog(rlZigzagCheckConsistency, "ZIndex cur: %d", cur->order);
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())
+ {
+ rError("In check_consistency(): ZNode %d not found in z_row of %d", cur->order, (*scur)->order);
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())
+ {
+ rError("In check_consistency(): ZNode %d not found in b_column of %d", cur->order, (*bcur)->order);
return false;
- if (cur->low != b_list.end() && cur->low->b_column.back() != cur) return false;
+ }
+ if (cur->low != b_list.end() && cur->low->b_column.back() != cur)
+ {
+ rError("low of ZNode %d is incorrect", cur->order);
+ return false;
+ }
+ if (cur->z_column.back()->low != cur)
+ {
+ rError("The low of the back of the z_column must be set correctly");
+ rError(" %d [%s], its back %d with low=%d", cur->order,
+ cur->z_column.tostring(out).c_str(),
+ cur->z_column.back()->order,
+ (cur->z_column.back()->low == z_list.end()) ? 0 : cur->z_column.back()->low->order);
+ return false;
+ }
}
for (BIndex cur = b_list.begin(); cur != b_list.end(); ++cur)
{
+ if (cur == b_skip) continue;
+
+ //rLog(rlZigzagCheckConsistency, "BIndex cur: %d", cur->order);
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())
+ {
+ rError("In check_consistency(): BNode %d not found in b_row of %d", cur->order, (*zcur)->order);
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())
+ {
+ rError("In check_consistency(): BNode %d not found in c_row of %d", cur->order, (*scur)->order);
return false;
+ }
+ if (!(cur->b_column.empty() || cur->b_column.back()->low == cur))
+ {
+ rError("The low of the back of the b_column must be set correctly");
+ return false;
+ }
+
+ // ZB == DC
+ ZColumn zb, dc;
+ std::for_each(cur->b_column.begin(), cur->b_column.end(), make_adder(&ZNode::z_column, zb));
+ std::for_each(cur->c_column.begin(), cur->c_column.end(), make_adder(&SimplexNode::boundary, dc));
+ zb.add(dc, cmp);
+ if (!zb.empty())
+ {
+ rError(" b_column: [%s]", cur->b_column.tostring(out).c_str());
+ rError(" c_column: [%s]", cur->c_column.tostring(out).c_str());
+ rError(" zb - dc: [%s]", zb.tostring(out).c_str());
+ rError("ZB = DC");
+ return false;
+ }
}
+#endif
return true;
}
@@ -455,6 +602,12 @@
ZigzagPersistence<BID>::
add_chain(IndexTo to, IndexFrom from, PrimaryMemberTo pmt, SecondaryMemberTo smt, PrimaryMemberFrom pmf)
{
+ rLog(rlZigzagAddChain, "Adding %d [%s] to %d [%s]",
+ (*from).order,
+ ((*from).*pmf).tostring(out).c_str(),
+ (*to).order,
+ ((*to).*pmt).tostring(out).c_str());
+
// Fix secondaries
std::for_each(make_intersection_iterator(((*from).*pmf).begin(), ((*from).*pmf).end(),
((*to).*pmt).begin(), ((*to).*pmt).end(),
@@ -470,7 +623,8 @@
((*to).*pmt).end(), ((*to).*pmt).end(),
cmp),
make_appender(smt, to));
-
+
// Add primaries
((*to).*pmt).add((*from).*pmf, cmp);
+ rLog(rlZigzagAddChain, "Got %s", ((*to).*pmt).tostring(out).c_str());
}
--- a/include/utilities/containers.h Fri Jan 02 14:54:15 2009 -0800
+++ b/include/utilities/containers.h Mon Jan 12 15:33:04 2009 -0800
@@ -1,9 +1,20 @@
#ifndef __CONTAINERS_H__
#define __CONTAINERS_H__
-#include <vector>
#include "circular_list.h"
+#if DEBUG_CONTAINERS
+ #include <debug/vector>
+ #include <debug/deque>
+ using std::__debug::vector;
+ using std::__debug::deque;
+#else
+ #include <vector>
+ #include <deque>
+ using std::vector;
+ using std::deque;
+#endif
+
// TODO: write documentation
template<class Container_, class Comparison_ = std::less<typename Container_::value_type> >
@@ -45,6 +56,8 @@
size_t size(const Container& c) const { return size_; }
void swap(SizeStorage& other) { std::swap(size_, other.size_); }
+ void clear() { size_ = 0; }
+
private:
size_t size_;
};
@@ -73,10 +86,10 @@
/* Specializations */
template<class T, class Comparison_>
-struct ContainerTraits<std::vector<T>, Comparison_>
+struct ContainerTraits<vector<T>, Comparison_>
{
typedef T Item;
- typedef std::vector<T> Container;
+ typedef vector<T> Container;
typedef typename Container::const_reference const_reference;
typedef Comparison_ Comparison;
@@ -86,6 +99,19 @@
};
template<class T, class Comparison_>
+struct ContainerTraits<deque<T>, Comparison_>
+{
+ typedef T Item;
+ typedef deque<T> Container;
+ typedef typename Container::const_reference const_reference;
+ typedef Comparison_ Comparison;
+
+ static void reserve(Container& c, size_t sz) { }
+ static void sort(Container& c, const Comparison& cmp = Comparison()) { std::sort(c.begin(), c.end(), cmp); }
+ static void push_front(Container& c, const_reference x) { c.push_front(x); }
+};
+
+template<class T, class Comparison_>
struct ContainerTraits<List<T>, Comparison_>
{
typedef T Item;
@@ -96,7 +122,7 @@
static void reserve(Container& c, size_t sz) { }
static void sort(Container& c, const Comparison& cmp = Comparison())
{
- std::vector<Item> tmp(c.begin(), c.end());
+ vector<Item> tmp(c.begin(), c.end());
std::sort(tmp.begin(), tmp.end(), cmp);
std::copy(tmp.begin(), tmp.end(), c.begin());
}
@@ -106,12 +132,12 @@
// TODO: specialize for List (singly-linked list)
template<class T>
-class SizeStorage<std::vector<T> >
+class SizeStorage<vector<T> >
{
public:
- typedef std::vector<T> Container;
+ typedef vector<T> Container;
typedef SizeStorage<Container> Self;
-
+
SizeStorage(size_t size = 0) {}
Self& operator+=(size_t inc) { return *this; }
@@ -122,6 +148,29 @@
Self operator--(int) { return *this; }
size_t size(const Container& c) const { return c.size(); }
void swap(SizeStorage& other) {}
+
+ void clear() {}
+};
+
+template<class T>
+class SizeStorage<deque<T> >
+{
+ public:
+ typedef deque<T> Container;
+ typedef SizeStorage<Container> Self;
+
+ SizeStorage(size_t size = 0) {}
+
+ Self& operator+=(size_t inc) { return *this; }
+ Self& operator-=(size_t dec) { return *this; }
+ Self& operator++() { return *this; }
+ Self operator++(int) { return *this; }
+ Self& operator--() { return *this; }
+ Self operator--(int) { return *this; }
+ size_t size(const Container& c) const { return c.size(); }
+ void swap(SizeStorage& other) {}
+
+ void clear() {}
};
#endif // __CONTAINERS_H__
--- a/include/utilities/counter.h Fri Jan 02 14:54:15 2009 -0800
+++ b/include/utilities/counter.h Mon Jan 12 15:33:04 2009 -0800
@@ -12,6 +12,7 @@
#define Count(x)
#define CountNum(x,y)
#define CountBy(x,y)
+ #define CountNumBy(x,y,z)
#define SetFrequency(x, freq)
#define SetTrigger(x, y)
#else // COUNTERS
@@ -19,6 +20,7 @@
#define Count(x) do { x->count++; if ((x->count % x->frequency == 0)) x->trigger->print(); } while (0)
#define CountNum(x,y) do { x->subcount[y]++; } while (0)
#define CountBy(x,y) do { x->count += y; } while (0)
+ #define CountNumBy(x,y,z) do { x->subcount[y] += z; } while (0)
#define SetFrequency(x, freq) do { x->frequency = freq; } while (0)
#define SetTrigger(x, y) do { x->trigger = y; } while(0)