Instrumented code for counting:
* added counters to addition in cohomology and ChainWrapper
* rips-pairwise-cohomology counts the maximum elements stored in the cycles
* added alphashapes3d-cohomology
* moved progress_display from DynamicPersistence to StaticPersistence
--- a/examples/alphashapes/CMakeLists.txt Fri Nov 06 14:19:08 2009 -0800
+++ b/examples/alphashapes/CMakeLists.txt Sat Nov 28 16:45:42 2009 -0800
@@ -10,6 +10,7 @@
if (CGAL_FOUND)
set (targets alphashapes3d
alphashapes2d
+ alphashapes3d-cohomology
#alpharadius
)
add_definitions (${CGAL_CXX_FLAGS_INIT})
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/alphashapes/alphashapes3d-cohomology.cpp Sat Nov 28 16:45:42 2009 -0800
@@ -0,0 +1,130 @@
+#include "alphashapes3d.h"
+#include "../cohomology/wrappers.h"
+
+#include <topology/cohomology-persistence.h>
+
+#include <utilities/log.h>
+#include <utilities/timer.h>
+
+#include <iostream>
+#include <fstream>
+
+#include <boost/program_options.hpp>
+#include <boost/progress.hpp>
+#include <boost/foreach.hpp>
+
+
+typedef boost::tuple<Dimension, RealType> BirthInfo;
+typedef CohomologyPersistence<BirthInfo> Persistence;
+
+typedef Persistence::SimplexIndex Index;
+typedef Persistence::Death Death;
+
+namespace po = boost::program_options;
+
+int main(int argc, char** argv)
+{
+#ifdef LOGGING
+ rlog::RLogInit(argc, argv);
+
+ stdoutLog.subscribeTo( RLOG_CHANNEL("info") );
+ stdoutLog.subscribeTo( RLOG_CHANNEL("error") );
+ //stdoutLog.subscribeTo( RLOG_CHANNEL("topology/persistence") );
+ //stdoutLog.subscribeTo( RLOG_CHANNEL("topology/chain") );
+#endif
+
+ std::string infilename, outfilename;
+
+ po::options_description hidden("Hidden options");
+ hidden.add_options()
+ ("input-file", po::value<std::string>(&infilename), "Point set whose alpha shape filtration and persistence we want to compute")
+ ("output-file", po::value<std::string>(&outfilename), "Where to write the collection of persistence diagrams");
+
+ po::positional_options_description pos;
+ pos.add("input-file", 1);
+ pos.add("output-file", 2);
+
+ po::options_description all; all.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("input-file") || !vm.count("output-file"))
+ {
+ std::cout << "Usage: " << argv[0] << " input-file output-file" << std::endl;
+ // std::cout << hidden << std::endl;
+ return 1;
+ }
+
+ std::ofstream diagram_out(outfilename.c_str());
+
+ Timer total_timer; total_timer.start();
+
+ // Read in the point set and compute its Delaunay triangulation
+ std::ifstream in(infilename.c_str());
+ double x,y,z;
+ Delaunay3D Dt;
+ while(in)
+ {
+ in >> x >> y >> z;
+ Point p(x,y,z);
+ Dt.insert(p);
+ }
+ rInfo("Delaunay triangulation computed");
+
+ // Set up the alpha shape filtration
+ AlphaSimplex3DVector complex;
+ fill_complex(Dt, complex);
+ rInfo("Simplices: %d", complex.size());
+ std::sort(complex.begin(), complex.end(), AlphaSimplex3D::AlphaOrder());
+
+ Timer persistence_timer; persistence_timer.start();
+ std::map<AlphaSimplex3D, Index, AlphaSimplex3D::VertexComparison> complex_map;
+ Persistence p;
+ boost::progress_display show_progress(complex.size());
+
+ #ifdef COUNTERS
+ Counter::CounterType max_element_count = 0;
+ #endif
+
+ for(AlphaSimplex3DVector::const_iterator cur = complex.begin(); cur != complex.end(); ++cur)
+ {
+ const AlphaSimplex3D& s = *cur;
+ std::vector<Index> boundary;
+ for (AlphaSimplex3D::BoundaryIterator bcur = s.boundary_begin(); bcur != s.boundary_end(); ++bcur)
+ boundary.push_back(complex_map[*bcur]);
+
+ Index idx; Death d;
+ bool store = s.dimension() < 3;
+ boost::tie(idx, d) = p.add(boundary.begin(), boundary.end(), boost::make_tuple(s.dimension(), s.value()), store);
+
+ // c[*cur] = idx;
+ if (store)
+ complex_map[s] = idx;
+
+ if (d && (s.value() - d->get<1>()) > 0)
+ {
+ AssertMsg(d->get<0>() == s.dimension() - 1, "Dimensions must match");
+ diagram_out << (s.dimension() - 1) << " " << d->get<1>() << " " << s.value() << std::endl;
+ }
+ ++show_progress;
+
+ #ifdef COUNTERS
+ max_element_count = std::max(max_element_count, cCohomologyElementCount->count);
+ #endif
+ }
+ // output infinte persistence pairs
+ for (Persistence::CocycleIndex cur = p.begin(); cur != p.end(); ++cur)
+ diagram_out << cur->birth.get<0>() << " " << cur->birth.get<1>() << " inf" << std::endl;
+ persistence_timer.stop();
+
+ total_timer.stop();
+ persistence_timer.check("Persistence timer");
+ total_timer.check("Total timer");
+
+ #ifdef COUNTERS
+ std::cout << "Max element count: " << max_element_count << std::endl;
+ #endif
+}
--- a/examples/alphashapes/alphashapes3d.cpp Fri Nov 06 14:19:08 2009 -0800
+++ b/examples/alphashapes/alphashapes3d.cpp Sat Nov 28 16:45:42 2009 -0800
@@ -30,8 +30,8 @@
//stdoutLog.subscribeTo( RLOG_CHANNEL("topology/chain") );
#endif
- SetFrequency(GetCounter("persistence/pair"), 10000);
- SetTrigger(GetCounter("persistence/pair"), GetCounter(""));
+ // SetFrequency(GetCounter("persistence/pair"), 10000);
+ // SetTrigger(GetCounter("persistence/pair"), GetCounter(""));
std::string infilename, outfilename;
--- a/examples/cohomology/rips-pairwise-cohomology.cpp Fri Nov 06 14:19:08 2009 -0800
+++ b/examples/cohomology/rips-pairwise-cohomology.cpp Sat Nov 28 16:45:42 2009 -0800
@@ -8,6 +8,11 @@
#include <utilities/property-maps.h>
#include <utilities/timer.h>
#include <utilities/log.h>
+#include <utilities/counter.h>
+#include <utilities/memory.h>
+
+#include <sys/resource.h>
+#include <malloc.h>
#include <string>
@@ -90,6 +95,19 @@
ZpField zp(prime);
Persistence p(zp);
boost::progress_display show_progress(v.size());
+
+ #ifdef COUNTERS
+ Counter::CounterType max_element_count = 0;
+ unsigned max_memory = 0;
+ long max_rss = 0;
+ long max_ixrss = 0;
+ long max_idrss = 0;
+ long max_isrss = 0;
+
+ int max_uordblks = 0;
+ int max_fordblks = 0;
+ #endif
+
for (unsigned j = 0; j < index_in_v.size(); ++j)
{
SimplexVector::const_iterator cur = v.begin() + index_in_v[j];
@@ -111,6 +129,22 @@
diagram_out << (cur->dimension() - 1) << " " << d->get<1>() << " " << size(*cur) << std::endl;
}
++show_progress;
+
+ #ifdef COUNTERS
+ max_element_count = std::max(max_element_count, cCohomologyElementCount->count);
+ // max_memory = std::max(max_memory, report_memory());
+
+ // struct rusage usage;
+ // getrusage(RUSAGE_SELF, &usage);
+ // max_rss = std::max(max_rss, usage.ru_maxrss);
+ // max_ixrss = std::max(max_ixrss, usage.ru_ixrss);
+ // max_idrss = std::max(max_idrss, usage.ru_idrss);
+ // max_isrss = std::max(max_isrss, usage.ru_isrss);
+
+ // struct mallinfo info = mallinfo();
+ // max_uordblks = std::max(max_uordblks, info.uordblks);
+ // max_fordblks = std::max(max_fordblks, info.fordblks);
+ #endif
}
// output infinte persistence pairs
for (Persistence::CocycleIndex cur = p.begin(); cur != p.end(); ++cur)
@@ -120,18 +154,28 @@
// p.show_cocycles();
// Output alive cocycles of dimension 1
- unsigned i = 0;
- for (Persistence::Cocycles::const_iterator cur = p.begin(); cur != p.end(); ++cur)
+ if (!cocycle_prefix.empty())
{
- if (cur->birth.get<0>() != 1) continue;
- output_cocycle(cocycle_prefix, i, v, *cur, prime);
- // std::cout << "Cocycle of dimension: " << cur->birth.get<0>() << " born at " << cur->birth.get<1>() << std::endl;
- ++i;
+ unsigned i = 0;
+ for (Persistence::Cocycles::const_iterator cur = p.begin(); cur != p.end(); ++cur)
+ {
+ if (cur->birth.get<0>() != 1) continue;
+ output_cocycle(cocycle_prefix, i, v, *cur, prime);
+ // std::cout << "Cocycle of dimension: " << cur->birth.get<0>() << " born at " << cur->birth.get<1>() << std::endl;
+ ++i;
+ }
}
total_timer.stop();
rips_timer.check("Rips timer");
persistence_timer.check("Persistence timer");
total_timer.check("Total timer");
+
+ #ifdef COUNTERS
+ std::cout << "Max element count: " << max_element_count << std::endl;
+ // std::cout << "Max memory use: " << max_memory << " kB" << std::endl;
+ // std::cout << "Max RSS: " << max_rss << " " << max_ixrss << " " << max_idrss << " " << max_isrss << std::endl;
+ // std::cout << "Max Blks: " << max_uordblks << " " << max_fordblks << std::endl;
+ #endif
}
void program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, ZpField::Element& prime, std::string& boundary_name, std::string& cocycle_prefix, std::string& vertices_name, std::string& diagram_name)
--- a/examples/rips/rips-pairwise.cpp Fri Nov 06 14:19:08 2009 -0800
+++ b/examples/rips/rips-pairwise.cpp Sat Nov 28 16:45:42 2009 -0800
@@ -23,19 +23,21 @@
typedef Generator::Simplex Smplx;
typedef std::vector<Smplx> SimplexVector;
typedef Filtration<SimplexVector, unsigned> Fltr;
-//typedef StaticPersistence<> Persistence;
-typedef DynamicPersistenceChains<> Persistence;
+typedef StaticPersistence<> Persistence;
+// typedef DynamicPersistenceChains<> Persistence;
typedef PersistenceDiagram<> PDgm;
-void program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance);
+void program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, std::string& diagram_name);
int main(int argc, char* argv[])
{
Dimension skeleton;
DistanceType max_distance;
- std::string infilename;
+ std::string infilename, diagram_name;
- program_options(argc, argv, infilename, skeleton, max_distance);
+ program_options(argc, argv, infilename, skeleton, max_distance, diagram_name);
+ std::ofstream diagram_out(diagram_name.c_str());
+ std::cout << "Diagram: " << diagram_name << std::endl;
PointContainer points;
read_points(infilename, points);
@@ -48,7 +50,7 @@
// 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;
+ 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));
@@ -73,7 +75,7 @@
// if (b.dimension() != 1) continue;
// std::cout << "Pair: (" << size(b) << ", " << size(d) << ")" << std::endl;
if (b.dimension() >= skeleton) continue;
- std::cout << b.dimension() << " " << size(b) << " " << size(d) << std::endl;
+ diagram_out << b.dimension() << " " << size(b) << " " << size(d) << std::endl;
} else if (cur->pair == cur) // positive could be unpaired
{
const Smplx& b = f.simplex(f.begin() + (cur - p.begin()));
@@ -82,7 +84,7 @@
// std::cout << "Unpaired birth: " << size(b) << std::endl;
// cycle = cur->chain;
if (b.dimension() >= skeleton) continue;
- std::cout << b.dimension() << " " << size(b) << " inf" << std::endl;
+ diagram_out << b.dimension() << " " << size(b) << " inf" << std::endl;
}
// Iterate over the cycle
@@ -100,7 +102,7 @@
persistence_timer.check("# Persistence timer");
}
-void program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance)
+void program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance, std::string& diagram_name)
{
namespace po = boost::program_options;
@@ -112,7 +114,13 @@
visible.add_options()
("help,h", "produce help message")
("skeleton-dimsnion,s", po::value<Dimension>(&skeleton)->default_value(2), "Dimension of the Rips complex we want to compute")
- ("max-distance,m", po::value<DistanceType>(&max_distance)->default_value(Infinity), "Maximum value for the Rips complex construction");
+ ("max-distance,m", po::value<DistanceType>(&max_distance)->default_value(Infinity), "Maximum value for the Rips complex construction")
+ ("diagram,d", po::value<std::string>(&diagram_name), "Filename where to output the persistence diagram");
+#if LOGGING
+ std::vector<std::string> log_channels;
+ visible.add_options()
+ ("log,l", po::value< std::vector<std::string> >(&log_channels), "log channels to turn on (info, debug, etc)");
+#endif
po::positional_options_description pos;
pos.add("input-file", 1);
@@ -124,6 +132,11 @@
options(all).positional(pos).run(), vm);
po::notify(vm);
+#if LOGGING
+ for (std::vector<std::string>::const_iterator cur = log_channels.begin(); cur != log_channels.end(); ++cur)
+ stderrLog.subscribeTo( RLOG_CHANNEL(cur->c_str()) );
+#endif
+
if (vm.count("help") || !vm.count("input-file"))
{
std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl;
--- a/include/topology/chain.hpp Fri Nov 06 14:19:08 2009 -0800
+++ b/include/topology/chain.hpp Sat Nov 28 16:45:42 2009 -0800
@@ -141,6 +141,8 @@
CountingBackInserter<ChainRepresentation> bi(tmp);
std::set_symmetric_difference(begin(), end(), c.begin(), c.end(), bi, cmp);
+
+ CountBy(cChainAddBasic, size() + c.size() - (size() + c.size() - tmp.size())/2);
static_cast<ChainRepresentation*>(this)->swap(tmp);
static_cast<Size*>(this)->swap(bi);
--- a/include/topology/cohomology-persistence.hpp Fri Nov 06 14:19:08 2009 -0800
+++ b/include/topology/cohomology-persistence.hpp Sat Nov 28 16:45:42 2009 -0800
@@ -5,6 +5,7 @@
#include <utilities/log.h>
#include <utilities/indirect.h>
+#include <utilities/counter.h>
#include <boost/iterator/transform_iterator.hpp>
#include <boost/iterator/counting_iterator.hpp>
@@ -18,6 +19,12 @@
static rlog::RLogChannel* rlCohomology = DEF_CHANNEL("topology/cohomology", rlog::Log_Debug);
#endif
+#ifdef COUNTERS
+static Counter* cCohomologyAddBasic = GetCounter("cohomology/add/basic");
+static Counter* cCohomologyAddComparison = GetCounter("cohomology/add/comparison");
+static Counter* cCohomologyElementCount = GetCounter("cohomology/elements");
+#endif // COUNTERS
+
template<class BirthInfo, class SimplexData, class Field>
class CohomologyPersistence<BirthInfo, SimplexData, Field>::CompareSNode
{
@@ -67,7 +74,7 @@
FieldElement coefficient = field_.init(*coefficient_iter++);
SimplexIndex cursi = *cur;
- rLog(rlCohomology, " %d %d", cursi->order, sign);
+ rLog(rlCohomology, " %d %d", cursi->order, coefficient);
BOOST_FOREACH(const SNode& zcur, std::make_pair(cursi->row.begin(), cursi->row.end()))
candidates_bulk.push_back(std::make_pair(zcur.ci, field_.mul(coefficient, zcur.coefficient)));
}
@@ -140,6 +147,8 @@
cocycle.push_back(SNode(si, field_.id(), nw));
si->row.push_back(cocycle.front());
rLog(rlCohomology, " Cocyle: %d", si->order);
+
+ Count(cCohomologyElementCount);
return std::make_pair(si, Death());
}
@@ -166,11 +175,16 @@
// add z to everything else in candidates
for (typename Candidates::iterator cur = boost::next(candidates.begin());
cur != candidates.end(); ++cur)
+ {
+ CountBy(cCohomologyElementCount, -cur->first->zcolumn.size());
add_cocycle(*cur, z);
+ CountBy(cCohomologyElementCount, cur->first->zcolumn.size());
+ }
for (typename ZColumn::iterator cur = z.first->zcolumn.begin(); cur != z.first->zcolumn.end(); ++cur)
cur->unlink();
+ CountBy(cCohomologyElementCount, -z.first->zcolumn.size());
cocycles_.erase(candidates.front().first);
return std::make_pair(si, d);
@@ -208,6 +222,8 @@
while (tcur != to.first->zcolumn.end() && fcur != from.first->zcolumn.end())
{
rLog(rlCohomology, " %d %d", tcur->si->order, fcur->si->order);
+ Count(cCohomologyAddComparison);
+ Count(cCohomologyAddBasic);
if (cmp(*tcur, *fcur))
{
nw.push_back(*tcur);
@@ -230,11 +246,13 @@
for (; tcur != to.first->zcolumn.end(); ++tcur)
{
rLog(rlCohomology, " %d", tcur->si->order);
+ Count(cCohomologyAddBasic);
nw.push_back(SNode(*tcur));
}
for (; fcur != from.first->zcolumn.end(); ++fcur)
{
rLog(rlCohomology, " %d", fcur->si->order);
+ Count(cCohomologyAddBasic);
nw.push_back(SNode(fcur->si, field_.mul(multiplier, fcur->coefficient), ci));
}
--- a/include/topology/dynamic-persistence.h Fri Nov 06 14:19:08 2009 -0800
+++ b/include/topology/dynamic-persistence.h Sat Nov 28 16:45:42 2009 -0800
@@ -4,8 +4,6 @@
#include "static-persistence.h"
#include <utilities/types.h>
-#include <boost/progress.hpp>
-
#ifdef COUNTERS
static Counter* cTrailLength = GetCounter("persistence/pair/traillength"); // the size of matrix U in RU decomposition
static Counter* cChainLength = GetCounter("persistence/pair/chainlength"); // the size of matrix V in R=DV decomposition
@@ -274,17 +272,14 @@
{
// TODO: this is specialized for std::vector
PairingChainsVisitor(OrderIndex bg, ConsistencyComparison ccmp, unsigned size):
- bg_(bg), ccmp_(ccmp), show_progress(size) {}
+ Parent::PairVisitor(size), bg_(bg), ccmp_(ccmp) {}
void init(OrderIndex i) const { i->consistency = i - bg_; i->chain.append(i, ccmp_); }
void update(OrderIndex j, OrderIndex i) const { j->chain.add(i->pair->chain, ccmp_); }
- void finished(OrderIndex i) const { CountBy(cChainLength, i->chain.size()); ++show_progress; }
+ void finished(OrderIndex i) const { Parent::PairVisitor::finished(i); CountBy(cChainLength, i->chain.size()); }
OrderIndex bg_;
ConsistencyComparison ccmp_;
-
- mutable boost::progress_display
- show_progress;
};
ConsistencyComparison ccmp_;
--- a/include/topology/static-persistence.h Fri Nov 06 14:19:08 2009 -0800
+++ b/include/topology/static-persistence.h Sat Nov 28 16:45:42 2009 -0800
@@ -7,6 +7,8 @@
#include <utilities/types.h>
+#include <boost/progress.hpp>
+
template<class Data_, class ChainTraits_, class ContainerTraits_, class Element_ = use_default>
struct PairCycleData: public Data_
@@ -80,7 +82,7 @@
// Function: pair_simplices()
// Compute persistence of the filtration
- void pair_simplices() { pair_simplices<PairVisitor>(begin(), end()); }
+ void pair_simplices() { pair_simplices<PairVisitor>(begin(), end(), PairVisitor(size())); }
// Functions: Accessors
// begin() - returns OrderIndex of the first element
@@ -101,6 +103,7 @@
// Acts as an archetype and if necessary a base class for visitors passed to <pair_simplices(bg, end, visitor)>.
struct PairVisitor
{
+ PairVisitor(unsigned size): show_progress(size) {}
// Function: init(i)
// Called after OrderElement pointed to by `i` has been initialized
// (its cycle is set to be its boundary, and pair is set to self, i.e. `i`)
@@ -114,7 +117,10 @@
// Function: finished(j)
// Called after the processing of `j` is finished.
- void finished(OrderIndex j) const {}
+ void finished(OrderIndex j) const { ++show_progress; }
+
+ mutable boost::progress_display
+ show_progress;
};
const Order& order() const { return order_; }
--- a/include/utilities/counter.h Fri Nov 06 14:19:08 2009 -0800
+++ b/include/utilities/counter.h Sat Nov 28 16:45:42 2009 -0800
@@ -15,6 +15,7 @@
#define CountNumBy(x,y,z)
#define SetFrequency(x, freq)
#define SetTrigger(x, y)
+ #define Print(x)
#else // COUNTERS
#define GetCounter(path) get_counter(path)
#define Count(x) do { x->count++; if ((x->count % x->frequency == 0)) x->trigger->print(); } while (0)
@@ -23,6 +24,7 @@
#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)
+ #define Print(x) do { x->trigger->print(); } while(0)
#include <map>
--- a/include/utilities/memory.h Fri Nov 06 14:19:08 2009 -0800
+++ b/include/utilities/memory.h Sat Nov 28 16:45:42 2009 -0800
@@ -11,7 +11,7 @@
static rlog::RLogChannel* rlMemory = DEF_CHANNEL("memory", rlog::Log_Debug);
-void report_memory()
+unsigned report_memory()
{
pid_t pid = getpid();
std::stringstream smaps_name;
@@ -30,11 +30,13 @@
}
}
rLog(rlMemory, "Private memory usage: %d kB", memory);
+
+ return memory;
}
#else
-void report_memory() {}
+unsigned report_memory() { return 0; }
#endif // LOGGING