--- a/bindings/python/cohomology-persistence.cpp Sun Dec 20 07:45:29 2009 -0800
+++ b/bindings/python/cohomology-persistence.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -2,6 +2,7 @@
#include <boost/python.hpp>
#include <boost/python/stl_iterator.hpp>
+#include <boost/python/overloads.hpp>
#include <boost/shared_ptr.hpp>
namespace bp = boost::python;
@@ -19,16 +20,34 @@
return chp;
}
-bp::tuple chp_add(dp::CohomPersistence& chp, bp::object bdry, dp::BirthID birth)
+
+bp::tuple chp_add(dp::CohomPersistence& chp,
+ bp::object bdry,
+ dp::BirthID birth,
+ bool store,
+ bool image,
+ bp::object coefficients)
{
dp::CohomPersistence::SimplexIndex i;
dp::CohomPersistence::Death d;
- boost::tie(i,d) = chp.add(bp::stl_input_iterator<dp::CohomPersistence::SimplexIndex>(bdry),
- bp::stl_input_iterator<dp::CohomPersistence::SimplexIndex>(),
- birth);
+
+ if (coefficients)
+ {
+ boost::tie(i,d) = chp.add(bp::stl_input_iterator<int>(coefficients),
+ bp::stl_input_iterator<dp::CohomPersistence::SimplexIndex>(bdry),
+ bp::stl_input_iterator<dp::CohomPersistence::SimplexIndex>(),
+ birth, store, dp::CohomPersistence::SimplexData(), image);
+ } else
+ {
+ boost::tie(i,d) = chp.add(bp::stl_input_iterator<dp::CohomPersistence::SimplexIndex>(bdry),
+ bp::stl_input_iterator<dp::CohomPersistence::SimplexIndex>(),
+ birth, store, dp::CohomPersistence::SimplexData(), image);
+ }
+
return bp::make_tuple(i,d);
}
+
dp::CohomPersistence::ZColumn::const_iterator
cocycle_zcolumn_begin(dp::CohomPersistence::Cocycle& ccl)
{ return ccl.zcolumn.begin(); }
@@ -46,6 +65,8 @@
}
+BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(add_overloads, add, 2, 4)
+
void export_cohomology_persistence()
{
bp::class_<dp::CohomPersistence::SimplexIndex>("CHSimplexIndex")
@@ -59,7 +80,7 @@
bp::class_<dp::CohomPersistence>("CohomologyPersistence")
.def("__init__", bp::make_constructor(&init_from_prime))
- .def("add", &chp_add)
+ .def("add", &chp_add, (bp::arg("bdry"), bp::arg("birth"), bp::arg("store")=true, bp::arg("image")=true, bp::arg("coefficients")=false))
.def("__iter__", bp::range(&dp::CohomPersistence::begin, &dp::CohomPersistence::end))
;
--- a/doc/get-build-install.rst Sun Dec 20 07:45:29 2009 -0800
+++ b/doc/get-build-install.rst Sun Dec 20 08:27:01 2009 -0800
@@ -1,3 +1,5 @@
+.. _download:
+
Get, Build, Install
===================
@@ -35,14 +37,25 @@
:Boost_: C++ utilities (version :math:`\geq` 1.36; including Boost.Python used to create
Python bindings)
-There also seems to be a dependence on the version of GCC, although I don't
-entirely understand it. GCC :math:`\geq` 4.3 definitely works, but some versions
-below that seem not to. There were some major changes in 4.3, so it's not
-entirely surprising, but I don't understand all the subtleties.
+.. warning::
+
+ There also seems to be a dependence on the version of GCC, although I don't
+ entirely understand it. GCC 4.3 and above definitely work, but some versions
+ below that seem not to.
+
+ One particular catch is that the default compiler on many current Mac OS X
+ is GCC 4.0 that has a well-known bug making Dionysus unusable. Fortunately
+ the problem is easy to solve by using GCC 4.2 that is often available on a
+ Mac under the name ``gcc-4.2``.
+
+ One can check the compiler version with ``g++ --version`` command.
Optional dependencies:
:CGAL_: for alpha shapes (version :math:`\geq` 3.4)
+ :CVXOPT_: for :ref:`circle-valued parametrization <cohomology-parametrization>` using LSQR
+ :PyX_: :sfile:`tools/draw-diagram/draw.py` uses `PyX`_ to
+ produce a PDF of the diagram
:rlog_: used for logging only (not needed by default)
.. :dsrpdb_: for reading PDB files
@@ -52,6 +65,8 @@
.. _CMake: http://www.cmake.org
.. _Boost: http://www.boost.org
.. _CGAL: http://www.cgal.org
+.. _CVXOPT: http://abel.ee.ucla.edu/cvxopt/
+.. _PyX: http://pyx.sourceforge.net/
.. _rlog: http://www.arg0.net/rlog
.. _dsrpdb: http://www.salilab.org/~drussel/pdb/
.. _SYNAPS: http://www-sop.inria.fr/galaad/synaps/
@@ -68,6 +83,11 @@
cmake ..
make
+.. tip::
+
+ To use GCC 4.2 on a Mac one can try ``CXX=g++-4.2 cmake ..`` instead of
+ ``cmake ..``.
+
Instead of ``cmake``, one can run ``ccmake`` for a curses interface. The
following configuration options are available. One can set them either through
the curses interface or by passing a flag of the form ``-Doptimize:bool=on`` to
--- a/doc/index.rst Sun Dec 20 07:45:29 2009 -0800
+++ b/doc/index.rst Sun Dec 20 08:27:01 2009 -0800
@@ -17,14 +17,14 @@
* Persistent homology computation [ELZ02]_ [ZC05]_
* Vineyards [CEM06]_ |cpp-only|
-* Persistent cohomology computation (described in [dSVJ09]_) |cpp-only|
+* Persistent cohomology computation (described in [dSVJ09]_)
* Zigzag persistent homology [CdSM09]_
* :ref:`examples` provide useful functionality in and of themselves:
* :ref:`Alpha shape construction <alpha-shape-example>` in 2D and 3D
* :ref:`Rips complex construction <rips-example>`
* Cech complex construction |cpp-only|
- * :ref:`Circle-valued parametrization using persistent cohomology <cohomology-parametrization>`
+ * :ref:`Circle-valued parametrization <cohomology-parametrization>`
.. todo::
Document more examples.
@@ -39,6 +39,8 @@
simplicity, elegance, and interactivity of Python. Since they mimick the C++
functionality, their documentation may be a helpful resource for the latter.
+:ref:`Download <download>` the software, or read a :ref:`tutorial` to
+get acquainted with its structure.
.. include:: substitutions.aux
--- a/doc/python/cohomology-persistence.rst Sun Dec 20 07:45:29 2009 -0800
+++ b/doc/python/cohomology-persistence.rst Sun Dec 20 08:27:01 2009 -0800
@@ -11,12 +11,28 @@
this point on all the computation will be performed with coefficients
coming from :math:`\mathbb{Z}/prime \mathbb{Z}`.
- .. method:: add(boundary, birth)
+ .. method:: add(boundary, birth, [store = True], [image = True], [coefficients = []])
Adds a simplex with the given `boundary` to the complex, i.e.
:math:`K_{i+1} = K_i \cup \sigma` and `boundary` = :math:`\partial \sigma`.
If a new class is born as a result of the addition, `birth` is stored with
- it for future reference.
+ it for future reference.
+
+ If `store` is ``False`` and a class is born, it will not be stored in
+ :class:`CohomologyPersistence`. This is useful to not waste space on the
+ classes of the dimension equal to the maximum-dimensional simplices of
+ the complex since such classes will never die.
+
+ The `image` parameter allows one to work with a case of a space
+ :math:`L \subseteq K` where the filtration of :math:`K` induces a
+ filtration of :math:`L`. In this case, one may want to compute **image
+ persistence** (i.e. the persistence of the sequences of the images given
+ by the inclusion of :math:`L` in :math:`K`). `image` indicates whether
+ the simplex added belongs to :math:`L` or not.
+
+ If given, `coefficients` is a list parallel to `boundary` that provides
+ coefficients for the corresponding boundary elements. If empty, it is
+ assumed to be :math:`(-1)^i`.
:returns: a pair (`i`, `d`). The first element is the index `i`.
It is the internal representation of the newly added simplex,
--- a/doc/tutorial.rst Sun Dec 20 07:45:29 2009 -0800
+++ b/doc/tutorial.rst Sun Dec 20 08:27:01 2009 -0800
@@ -32,8 +32,9 @@
Each one has its attribute ``data`` set to a pair: the
smallest value of the squared distance function on the dual Voronoi cell and
whether the simplex is critical or not (i.e. whether its dual cell does or
-does not contain a critical point of the distance function). See :ref:`alphashapes`
-for more details.
+does not contain a critical point of the distance function).
+See :ref:`alphashapes` for more details, and :ref:`alpha-shape-example` for a
+full example.
As a result, if one wanted only those simplices whose alpha shape value did not
exceed 10, one could obtain them as follows::
@@ -70,7 +71,7 @@
the entire filtration at once, compute its persistence in one operation, and
then examine the pairings. The second way is to feed simplices one by one in an
*online* manner and manually keep track of the pairs that are created.
-:class:`ZigzagPersistence` and `CohomologyPersistence` |cpp-only| accept their
+:class:`ZigzagPersistence` and :class:`CohomologyPersistence` accept their
input this way,
@@ -121,7 +122,10 @@
^^^^^^
Class :class:`ZigzagPersistence` accepts additions and removals of the simplices
-one by one, and returns an internal representation of the simplices. When one
+one by one, and returns an internal representation of the simplices.
+(:class:`CohomologyPersistence` works in the same way, but accepts only
+additions of the simplices.)
+When one
adds a simplex via :meth:`ZigzagPersistence.add`, one must provide its boundary
in this internal representation together with a *birth value* which
:class:`ZigzagPersistence` will store in case a class is born as a result of the
@@ -136,7 +140,7 @@
complex = {}
zz = ZigzagPersistence()
for s in simplices:
- i,d = zz.add([complex[sb] for sb in s.boundary, (s.dimension(), s.data))
+ i,d = zz.add([complex[sb] for sb in s.boundary], (s.dimension(), s.data))
complex[s] = i
if d is not None: # we have a death
dimension, birth = d # we previously stored the (dimension, data) pair
@@ -159,7 +163,9 @@
:func:`add_simplices` and :func:`remove_simplices`.
See the :ref:`bindings reference <python-bindings>` for more details, and
-:ref:`triangle-zigzag-example` for an example.
+:ref:`triangle-zigzag-example` for an example of :class:`ZigzagPersistence`.
+See :ref:`rips-pairwise-cohomology` for an example of
+:class:`CohomologyPersistence`.
.. _speed-up-suggestions:
--- a/examples/alphashapes/CMakeLists.txt Sun Dec 20 07:45:29 2009 -0800
+++ b/examples/alphashapes/CMakeLists.txt Sun Dec 20 08:27:01 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 Sun Dec 20 08:27:01 2009 -0800
@@ -0,0 +1,131 @@
+#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
+ typedef std::vector<AlphaSimplex3D> AlphaSimplex3DVector;
+ 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 Sun Dec 20 07:45:29 2009 -0800
+++ b/examples/alphashapes/alphashapes3d.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -1,4 +1,5 @@
#include <utilities/log.h>
+#include <utilities/timer.h>
#include "alphashapes3d.h"
#include <topology/filtration.h>
@@ -30,8 +31,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;
@@ -82,8 +83,11 @@
Persistence p(af);
rInfo("Persistence initializaed");
+ Timer persistence_timer; persistence_timer.start();
p.pair_simplices();
+ persistence_timer.stop();
rInfo("Simplices paired");
+ persistence_timer.check("Persistence timer");
Persistence::SimplexMap<AlphaFiltration> m = p.make_simplex_map(af);
std::map<Dimension, PDgm> dgms;
--- a/examples/cohomology/CMakeLists.txt Sun Dec 20 07:45:29 2009 -0800
+++ b/examples/cohomology/CMakeLists.txt Sun Dec 20 08:27:01 2009 -0800
@@ -1,6 +1,7 @@
set (targets
rips-cohomology
rips-pairwise-cohomology
+ rips-weighted-cohomology
triangle-cohomology
)
--- a/examples/cohomology/cocycle.py Sun Dec 20 07:45:29 2009 -0800
+++ b/examples/cohomology/cocycle.py Sun Dec 20 08:27:01 2009 -0800
@@ -6,7 +6,7 @@
from sys import argv, exit
import os.path
-def smooth(boundary_list, cocycle_list, vertexmap_list):
+def smooth(boundary_list, cocycle_list, vertices):
dimension = max((max(d[1], d[2]) for d in boundary_list))
dimension += 1
@@ -46,8 +46,8 @@
print "Expected a harmonic cocycle:", sum((D*v)**2), sum((D.T*v)**2)
values = [None]*len(vertices)
- for i in xrange(len(vertices)):
- values[vertices[i]] = solution[0][i]
+ for i,v in vertices:
+ values[v] = solution[0][i]
return values
--- a/examples/cohomology/lsqr.py Sun Dec 20 07:45:29 2009 -0800
+++ b/examples/cohomology/lsqr.py Sun Dec 20 08:27:01 2009 -0800
@@ -24,7 +24,7 @@
c=0.
s=b/ab
r=ab
- elif ab>aa:
+ elif ab>=aa:
sb=1
if b<0: sb=-1
tau=a/b
--- a/examples/cohomology/output.h Sun Dec 20 07:45:29 2009 -0800
+++ b/examples/cohomology/output.h Sun Dec 20 08:27:01 2009 -0800
@@ -13,14 +13,16 @@
return cmp(s1, s2) || cmp(s2, s1);
}
-unsigned index(const SimplexVector& v, const Smplx& s, const Generator::Comparison& cmp)
+template<class Comparison>
+unsigned index(const SimplexVector& v, const Smplx& s, const Comparison& cmp)
{
SimplexVector::const_iterator it = std::lower_bound(v.begin(), v.end(), s, cmp);
while (neq(*it, s)) ++it;
return it - v.begin();
}
-void output_boundary_matrix(std::ostream& out, const SimplexVector& v, const Generator::Comparison& cmp)
+template<class Comparison>
+void output_boundary_matrix(std::ostream& out, const SimplexVector& v, const Comparison& cmp)
{
unsigned i = 0;
for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
@@ -49,7 +51,7 @@
}
}
-void output_cocycle(std::string cocycle_prefix, unsigned i, const SimplexVector& v, const Persistence::Cocycle& c, ZpField::Element prime, const Generator::Comparison& cmp)
+void output_cocycle(std::string cocycle_prefix, unsigned i, const SimplexVector& v, const Persistence::Cocycle& c, ZpField::Element prime)
{
std::ostringstream istr; istr << '-' << i;
std::string filename = cocycle_prefix + istr.str() + ".ccl";
@@ -57,8 +59,8 @@
out << "# Cocycle born at " << c.birth.get<1>() << std::endl;
for (Persistence::ZColumn::const_iterator zcur = c.zcolumn.begin(); zcur != c.zcolumn.end(); ++zcur)
{
- const Smplx& s = **(zcur->si);
+ //const Smplx& s = **(zcur->si);
out << (zcur->coefficient <= prime/2 ? zcur->coefficient : zcur->coefficient - prime) << " ";
- out << index(v, s, cmp) << "\n";
+ out << zcur->si->getValue() << "\n";
}
}
--- a/examples/cohomology/rips-pairwise-cohomology.cpp Sun Dec 20 07:45:29 2009 -0800
+++ b/examples/cohomology/rips-pairwise-cohomology.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -5,27 +5,37 @@
#include <geometry/distances.h>
#include <utilities/containers.h> // for BackInsertFunctor
+#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>
#include <boost/tuple/tuple.hpp>
#include <boost/program_options.hpp>
+#include <boost/progress.hpp>
+
+#include "wrappers.h"
typedef PairwiseDistances<PointContainer, L2Distance> PairDistances;
typedef PairDistances::DistanceType DistanceType;
typedef PairDistances::IndexType Vertex;
-typedef Rips<PairDistances> Generator;
+typedef boost::tuple<Dimension, DistanceType> BirthInfo;
+typedef CohomologyPersistence<BirthInfo, Wrapper<unsigned> > Persistence;
+typedef Persistence::SimplexIndex Index;
+typedef Persistence::Death Death;
+
+typedef Rips<PairDistances, Simplex<Vertex, Index> > Generator;
typedef Generator::Simplex Smplx;
typedef std::vector<Smplx> SimplexVector;
typedef SimplexVector::const_iterator SV_const_iterator;
-typedef boost::tuple<Dimension, DistanceType> BirthInfo;
-typedef CohomologyPersistence<BirthInfo, SV_const_iterator> Persistence;
-typedef Persistence::SimplexIndex Index;
-typedef Persistence::Death Death;
typedef std::map<Smplx, Index, Smplx::VertexComparison> Complex;
#include "output.h" // for output_*()
@@ -63,39 +73,78 @@
Generator::Evaluator size(distances);
Generator::Comparison cmp(distances);
SimplexVector v;
- Complex c;
Timer rips_timer; rips_timer.start();
rips.generate(skeleton, max_distance, make_push_back_functor(v));
- std::sort(v.begin(), v.end(), cmp);
+ std::sort(v.begin(), v.end(), Smplx::VertexComparison());
+
+ std::vector<unsigned> index_in_v(v.size());
+ for (unsigned idx = 0; idx < v.size(); ++idx)
+ index_in_v[idx] = idx;
+ std::sort(index_in_v.begin(), index_in_v.end(), IndirectIndexComparison<SimplexVector, Generator::Comparison>(v, cmp));
+
+ BinarySearchMap<Smplx, SimplexVector::iterator, Smplx::VertexComparison> map_of_v(v.begin(), v.end());
+
rips_timer.stop();
std::cout << "Simplex vector generated, size: " << v.size() << std::endl;
- output_boundary_matrix(bdry_out, v, cmp);
+ output_boundary_matrix(bdry_out, v, Smplx::VertexComparison());
output_vertex_indices(vertices_out, v);
Timer persistence_timer; persistence_timer.start();
ZpField zp(prime);
Persistence p(zp);
- for (SimplexVector::const_iterator cur = v.begin(); cur != v.end(); ++cur)
+ 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];
std::vector<Index> boundary;
for (Smplx::BoundaryIterator bcur = cur->boundary_begin(); bcur != cur->boundary_end(); ++bcur)
- boundary.push_back(c[*bcur]);
+ boundary.push_back(map_of_v[*bcur]->data());
Index idx; Death d;
bool store = cur->dimension() < skeleton;
- boost::tie(idx, d) = p.add(boundary.begin(), boundary.end(), boost::make_tuple(cur->dimension(), size(*cur)), store, cur);
+ boost::tie(idx, d) = p.add(boundary.begin(), boundary.end(), boost::make_tuple(cur->dimension(), size(*cur)), store, index_in_v[j]);
// c[*cur] = idx;
if (store)
- c.insert(std::make_pair(*cur, idx));
+ map_of_v[*cur]->data() = idx;
if (d && (size(*cur) - d->get<1>()) > 0)
{
AssertMsg(d->get<0>() == cur->dimension() - 1, "Dimensions must match");
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)
@@ -105,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, cmp);
- // 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/cohomology/rips-pairwise-cohomology.py Sun Dec 20 07:45:29 2009 -0800
+++ b/examples/cohomology/rips-pairwise-cohomology.py Sun Dec 20 08:27:01 2009 -0800
@@ -19,7 +19,7 @@
# While this step is unnecessary (Filtration below can be passed rips.cmp),
# it greatly speeds up the running times
for s in simplices: s.data = rips.eval(s)
- print time.asctime(), simplices[0], '...', simplices[-1]
+ print '#', time.asctime(), simplices[0], '...', simplices[-1]
simplices.sort(data_dim_cmp)
print '#', time.asctime(), "Simplices sorted"
@@ -28,7 +28,7 @@
complex = {}
for s in simplices:
- i,d = ch.add([complex[sb] for sb in s.boundary], (s.dimension(), s.data))
+ i,d = ch.add([complex[sb] for sb in s.boundary], (s.dimension(), s.data), store = (s.dimension() < skeleton))
complex[s] = i
if d:
dimension, birth = d
@@ -41,7 +41,7 @@
print dimension, birth, 'inf' # dimension, simplex data = birth
print "# Cocycle at (dim=%d, birth=%f)" % ccl.birth
for e in ccl:
- print "# ", e.si.order(), normalized(e.coefficient, prime)
+ print "# ", e.si.order, normalized(e.coefficient, prime)
def normalized(coefficient, prime):
if coefficient > prime/2:
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/cohomology/rips-weighted-cohomology.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -0,0 +1,184 @@
+#include <topology/cohomology-persistence.h>
+#include <topology/weighted-rips.h>
+
+#include <geometry/weighted-l2distance.h>
+#include <geometry/distances.h>
+
+#include <utilities/containers.h> // for BackInsertFunctor
+#include <utilities/property-maps.h>
+#include <utilities/timer.h>
+#include <utilities/log.h>
+
+#include <string>
+
+#include <boost/tuple/tuple.hpp>
+#include <boost/program_options.hpp>
+#include <boost/progress.hpp>
+
+#include "wrappers.h"
+
+typedef PairwiseDistances<PointContainer, WeightedL2Distance> PairDistances;
+typedef PairDistances::DistanceType DistanceType;
+typedef PairDistances::IndexType Vertex;
+
+typedef boost::tuple<Dimension, DistanceType> BirthInfo;
+typedef CohomologyPersistence<BirthInfo, Wrapper<unsigned> > Persistence;
+typedef Persistence::SimplexIndex Index;
+typedef Persistence::Death Death;
+
+typedef WeightedRips<PairDistances, Simplex<Vertex, Index> > Generator;
+typedef Generator::Simplex Smplx;
+typedef std::vector<Smplx> SimplexVector;
+typedef SimplexVector::const_iterator SV_const_iterator;
+
+#include "output.h" // for output_*()
+
+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);
+
+int main(int argc, char* argv[])
+{
+#ifdef LOGGING
+ rlog::RLogInit(argc, argv);
+
+ stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+#endif
+
+ Dimension skeleton;
+ DistanceType max_distance;
+ ZpField::Element prime;
+ std::string infilename, boundary_name, cocycle_prefix, vertices_name, diagram_name;
+
+ program_options(argc, argv, infilename, skeleton, max_distance, prime, boundary_name, cocycle_prefix, vertices_name, diagram_name);
+ std::ofstream bdry_out(boundary_name.c_str());
+ std::ofstream vertices_out(vertices_name.c_str());
+ std::ofstream diagram_out(diagram_name.c_str());
+ std::cout << "Boundary matrix: " << boundary_name << std::endl;
+ std::cout << "Cocycles: " << cocycle_prefix << "*.ccl" << std::endl;
+ std::cout << "Vertices: " << vertices_name << std::endl;
+ std::cout << "Diagram: " << diagram_name << std::endl;
+
+ Timer total_timer; total_timer.start();
+ PointContainer points;
+ read_weighted_points(infilename, points);
+
+ PairDistances distances(points);
+ Generator rips(distances);
+ Generator::Evaluator size(distances);
+ Generator::Comparison cmp(distances);
+ SimplexVector v;
+
+ Timer rips_timer; rips_timer.start();
+ rips.generate(skeleton, max_distance, make_push_back_functor(v));
+
+ /* Keep simplices sorted lexicographically (so that we can binary search through them) */
+ std::sort(v.begin(), v.end(), Smplx::VertexComparison());
+
+ /* We also need the simplices sorted by value though for the filtration:
+ index_in_v[j] refers to the simplex v[index_in_v[j]] */
+ std::vector<unsigned> index_in_v(v.size());
+ for (unsigned idx = 0; idx < v.size(); ++idx)
+ index_in_v[idx] = idx;
+ std::sort(index_in_v.begin(), index_in_v.end(), IndirectIndexComparison<SimplexVector, Generator::Comparison>(v, cmp));
+
+ /* Set up map access to the lexicographically sorted simplices */
+ BinarySearchMap<Smplx, SimplexVector::iterator, Smplx::VertexComparison> map_of_v(v.begin(), v.end());
+
+ rips_timer.stop();
+ std::cout << "Simplex vector generated, size: " << v.size() << std::endl;
+
+ /*output_boundary_matrix(bdry_out, v, Smplx::VertexComparison());
+ output_vertex_indices(vertices_out, v);*/
+
+ Timer persistence_timer; persistence_timer.start();
+ ZpField zp(prime);
+ Persistence p(zp);
+ boost::progress_display show_progress(v.size());
+ for (unsigned j = 0; j < index_in_v.size(); ++j)
+ {
+ SimplexVector::const_iterator cur = v.begin() + index_in_v[j];
+ std::vector<Index> boundary;
+ for (Smplx::BoundaryIterator bcur = cur->boundary_begin(); bcur != cur->boundary_end(); ++bcur)
+ boundary.push_back(map_of_v[*bcur]->data());
+
+ Index idx; Death d;
+ bool store = cur->dimension() < skeleton;
+ boost::tie(idx, d) = p.add(boundary.begin(), boundary.end(), boost::make_tuple(cur->dimension(), size(*cur)), store, index_in_v[j]);
+
+ if (store)
+ map_of_v[*cur]->data() = idx;
+
+ if (d && (size(*cur) - d->get<1>()) > 0)
+ {
+ AssertMsg(d->get<0>() == cur->dimension() - 1, "Dimensions must match");
+ diagram_out << (cur->dimension() - 1) << " " << d->get<1>() << " " << size(*cur) << std::endl;
+ }
+ ++show_progress;
+ }
+ // 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();
+
+
+ // 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 (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");
+}
+
+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)
+{
+ namespace po = boost::program_options;
+
+ po::options_description hidden("Hidden options");
+ hidden.add_options()
+ ("input-file", po::value<std::string>(&infilename), "Point set whose Rips zigzag we want to compute");
+
+ po::options_description visible("Allowed options", 100);
+ visible.add_options()
+ ("help,h", "produce help message")
+ ("skeleton-dimsnion,s", po::value<Dimension>(&skeleton)->default_value(2), "Dimension of the Rips complex we want to compute")
+ ("prime,p", po::value<ZpField::Element>(&prime)->default_value(11), "Prime p for the field F_p")
+ ("max-distance,m", po::value<DistanceType>(&max_distance)->default_value(Infinity), "Maximum value for the Rips complex construction")
+ ("boundary,b", po::value<std::string>(&boundary_name), "Filename where to output the boundary matrix")
+ ("cocycle,c", po::value<std::string>(&cocycle_prefix), "Prefix of the filename where to output the 1-dimensional cocycles")
+ ("vertices,v", po::value<std::string>(&vertices_name), "Filename where to output the simplex-vertex mapping")
+ ("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);
+
+ 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)
+ 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;
+ std::cout << visible << std::endl;
+ std::abort();
+ }
+}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/cohomology/wrappers.h Sun Dec 20 08:27:01 2009 -0800
@@ -0,0 +1,59 @@
+/**
+ * Wrapper class
+ *
+ * At points we need to wrap primitive data types in classes, so that we can
+ * pass them as template arguments to classes that like to inherit from their
+ * template arguments--this is done in CohomologyPersistence, for example.
+ */
+
+template<typename Primitive>
+class Wrapper
+{
+
+ public:
+
+ Wrapper () {}
+ Wrapper (Primitive v) { value = v; }
+
+ void setValue (const Primitive &v) { value = v; }
+ Primitive &getValue () { return value; }
+
+ /* provide seemless integration */
+ Wrapper &operator =(const Primitive &v) { setValue(v); return *this; }
+ operator Primitive() { return getValue; }
+
+ protected:
+
+ Primitive value;
+
+};
+
+/**
+ * IndirectIndexComparison class
+ *
+ * This class serves as a comparison function for arrays that are being sorted
+ * even though they only contain *indices* to the real data. Therefore, a reference
+ * to the original data as well as the data comparison function needs to be passed
+ * to this class for it to be functional.
+ */
+
+template<class DataContainer, class DataComparison>
+class IndirectIndexComparison
+{
+
+ public:
+
+ IndirectIndexComparison(const DataContainer &dstor, const DataComparison &dcmp) :
+ container(dstor), comparison(dcmp) { }
+
+ bool operator()(const unsigned &idx_1, const unsigned &idx_2) const
+ {
+ return comparison(container[idx_1], container[idx_2]);
+ }
+
+ private:
+
+ const DataContainer &container;
+ const DataComparison &comparison;
+
+};
--- a/examples/consistency/make-zigzag-subsamples.py Sun Dec 20 07:45:29 2009 -0800
+++ b/examples/consistency/make-zigzag-subsamples.py Sun Dec 20 08:27:01 2009 -0800
@@ -27,10 +27,10 @@
for c in count:
counts.append(' '.join(map(str, xrange(cur, cur+c))) + '\n')
cur += c
- counts.append(' '.join(map(str, xrange(cur-c, cur+c))) + '\n')
+ # counts.append(' '.join(map(str, xrange(cur-c, cur+c))) + '\n')
with open(subsamples_fn, 'w') as f:
- f.writelines(counts[:-1])
+ f.writelines(counts)
if __name__ == '__main__':
@@ -38,7 +38,7 @@
print "Usage: %s POINTS SUBSAMPLES POINTS1 [POINTS2 [POINTS3 [...]]]" % argv[0]
print
print "Creates a file POINTS with the union of POINTS* and SUBSAMPLES which lists"
- print "the indices of the points and their pairwise unions, one per line"
+ print "the indices of the points one per line"
exit()
points_fn = argv[1]
--- a/examples/consistency/rips-consistency-zigzag.cpp Sun Dec 20 07:45:29 2009 -0800
+++ b/examples/consistency/rips-consistency-zigzag.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -59,6 +59,8 @@
bool operator<(const BirthInfo& other) const { if (index == other.index) return (!un && other.un); else return index < other.index; }
bool operator>(const BirthInfo& other) const { return other.operator<(*this); }
+ bool operator>=(const BirthInfo& other) const { return !operator<(other); }
+ bool operator<=(const BirthInfo& other) const { return !operator>(other); }
Dimension dimension;
unsigned index;
@@ -121,6 +123,8 @@
while (line_in >> sample)
subsamples.back().push_back(sample);
}
+ AssertMsg(subsamples.front().size() == 0, "The first subsample should be empty");
+ AssertMsg(subsamples.back().size() == 0, "The last subsample should be empty"); // it's a convenient artifact of the stream processing above
std::cout << "Subsample size:" << std::endl;
for (unsigned i = 0; i < subsamples.size(); ++i)
@@ -141,7 +145,7 @@
SimplexVector subcomplex, across;
rInfo("Commencing computation");
- boost::progress_display show_progress(subsamples.size());
+ boost::progress_display show_progress(subsamples.size() - 1);
for (unsigned i = 0; i < subsamples.size() - 1; ++i)
{
// Take union of subsamples i and i+1
@@ -251,7 +255,7 @@
void report_death(std::ostream& out, const BirthInfo& birth, const BirthInfo& death, unsigned skeleton_dimension)
{
- if (birth.dimension < skeleton_dimension && death > birth)
+ if (birth.dimension < skeleton_dimension && death >= birth)
out << birth << " --- " << death << std::endl;
}
--- a/examples/rips/CMakeLists.txt Sun Dec 20 07:45:29 2009 -0800
+++ b/examples/rips/CMakeLists.txt Sun Dec 20 08:27:01 2009 -0800
@@ -1,6 +1,7 @@
set (targets
rips
rips-pairwise
+ rips-weighted
rips-image-zigzag
rips-zigzag)
--- a/examples/rips/rips-pairwise.cpp Sun Dec 20 07:45:29 2009 -0800
+++ b/examples/rips/rips-pairwise.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -26,15 +26,17 @@
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);
@@ -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->unpaired()) // positive could be unpaired
{
const Smplx& b = m[cur];
@@ -82,7 +84,7 @@
// std::cout << "Unpaired birth: " << size(b) << std::endl;
// cycle = cur->chain; // TODO
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
@@ -101,7 +103,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;
@@ -113,7 +115,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);
@@ -125,6 +133,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;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/rips/rips-weighted.cpp Sun Dec 20 08:27:01 2009 -0800
@@ -0,0 +1,144 @@
+#include <topology/weighted-rips.h>
+#include <topology/filtration.h>
+#include <topology/static-persistence.h>
+#include <topology/dynamic-persistence.h>
+#include <topology/persistence-diagram.h>
+
+//#define RIPS_CLOSURE_CECH_SKELETON
+#ifndef RIPS_CLOSURE_CECH_SKELETON
+#include <geometry/weighted-l2distance.h>
+#else
+#include <geometry/weighted-cechdistance.h>
+#endif
+
+#include <geometry/distances.h>
+
+#include <utilities/containers.h> // for BackInsertFunctor
+#include <utilities/timer.h>
+
+#include <vector>
+
+#include <boost/program_options.hpp>
+
+
+typedef PairwiseDistances<PointContainer, WeightedL2Distance> PairDistances;
+typedef PairDistances::DistanceType DistanceType;
+typedef PairDistances::IndexType Vertex;
+
+typedef WeightedRips<PairDistances> Generator;
+typedef Generator::Simplex Smplx;
+typedef Filtration<Smplx> Fltr;
+typedef StaticPersistence<> Persistence;
+//typedef DynamicPersistenceChains<> Persistence;
+typedef PersistenceDiagram<> PDgm;
+
+void program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance);
+
+int main(int argc, char* argv[])
+{
+#ifdef LOGGING
+ rlog::RLogInit(argc, argv);
+
+ stderrLog.subscribeTo( RLOG_CHANNEL("error") );
+#endif
+
+ Dimension skeleton;
+ DistanceType max_distance;
+ std::string infilename;
+
+ program_options(argc, argv, infilename, skeleton, max_distance);
+
+ PointContainer points;
+ read_weighted_points(infilename, points);
+
+ PairDistances distances(points);
+ Generator rips(distances);
+ Generator::Evaluator size(distances);
+ Generator::Comparison cmp (distances);
+ Fltr complex;
+
+ // Generate skeleton of the weighted Rips complex for epsilon = 50
+ rips.generate(skeleton, max_distance, make_push_back_functor(complex));
+ std::cout << "# Generated complex of size: " << complex.size() << std::endl;
+
+ // Generate filtration with respect to distance and compute its persistence
+ complex.sort(cmp);
+
+ Timer persistence_timer; persistence_timer.start();
+ Persistence p(complex);
+ p.pair_simplices();
+ persistence_timer.stop();
+
+ // Output cycles
+ Persistence::SimplexMap<Fltr> m = p.make_simplex_map(complex);
+ for (Persistence::iterator cur = p.begin(); cur != p.end(); ++cur)
+ {
+ const Persistence::Cycle& cycle = cur->cycle;
+
+ if (!cur->sign()) // only negative simplices have non-empty cycles
+ {
+ Persistence::OrderIndex birth = cur->pair; // the cycle that cur killed was born when we added birth (another simplex)
+
+ const Smplx& b = m[birth];
+ const Smplx& d = m[cur];
+
+ if (b.dimension() >= skeleton) continue;
+ std::cout << b.dimension() << " " << size(b) << " " << size(d) << std::endl;
+ } else if (cur->unpaired()) // positive could be unpaired
+ {
+ const Smplx& b = m[cur];
+ if (b.dimension() >= skeleton) continue;
+
+ std::cout << b.dimension() << " " << size(b) << " inf" << std::endl;
+ //cycle = cur->chain;
+ }
+ }
+
+ persistence_timer.check("# Persistence timer");
+}
+
+void program_options(int argc, char* argv[], std::string& infilename, Dimension& skeleton, DistanceType& max_distance)
+{
+ namespace po = boost::program_options;
+
+ po::options_description hidden("Hidden options");
+ hidden.add_options()
+ ("input-file", po::value<std::string>(&infilename), "Point set whose weighed Rips zigzag we want to compute");
+
+ po::options_description visible("Allowed options", 100);
+ visible.add_options()
+ ("help,h", "produce help message")
+ ("skeleton-dimsnion,s", po::value<Dimension>(&skeleton)->default_value(2), "Dimension of the weighted Rips complex we want to compute")
+ ("max-distance,m", po::value<DistanceType>(&max_distance)->default_value(Infinity), "Maximum value for the weighted Rips complex construction");
+#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);
+
+ 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)
+ stderrLog.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;
+ std::abort();
+ }
+}
--- a/examples/triangle/triangle-zigzag.py Sun Dec 20 07:45:29 2009 -0800
+++ b/examples/triangle/triangle-zigzag.py Sun Dec 20 08:27:01 2009 -0800
@@ -27,9 +27,9 @@
b += 1
# Remove all the simplices
-for s in sorted(complex.keys(), reverse = True):
+for s in sorted(complex.keys(), data_cmp, reverse = True):
print "%d: Removing %s" % (b, s)
d = zz.remove(complex[s], b)
- complex[s] = None
+ del complex[s]
if d: print "Interval (%d, %d)" % (d, b-1)
b += 1
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/weighted-cechdistance.h Sun Dec 20 08:27:01 2009 -0800
@@ -0,0 +1,55 @@
+#ifndef __L2_DISTANCE_H__
+#define __L2_DISTANCE_H__
+
+#include <utilities/types.h>
+
+#include <vector>
+#include <fstream>
+#include <functional>
+#include <cmath>
+#include <string>
+#include <sstream>
+
+
+typedef std::vector<double> Point;
+typedef std::vector<Point> PointContainer;
+
+struct WeightedL2Distance:
+ public std::binary_function<const Point&, const Point&, double>
+{
+ result_type operator()(const Point& p1, const Point& p2) const
+ {
+ AssertMsg(p1.size() == p2.size(), "Points must be in the same dimension (in L2Distance): dim1=%d, dim2=%d", p1.size()-1, p2.size()-1);
+
+ /* the distance of a point to itself is the radius at which the "power distance" ball around it appears:
+ d(p,p) := sqrt(-w_p) */
+ if (p1 == p2)
+ return sqrt(-p1[p1.size()-1]);
+
+ /* otherwise, the distance is the radius at which the power balls around p, q intersect:
+ d(p,q) := sqrt( ( (w_p-w_q)^2 / 4||p-q||^2 ) - (w_p+w_q) / 2 + ||p-q||^2 / 4 ) */
+ result_type sq_l2dist = 0;
+ result_type wp = p1[p1.size()-1];
+ result_type wq = p2[p2.size()-1];
+ for (size_t i = 0; i < p1.size()-1; ++i)
+ sq_l2dist += (p1[i] - p2[i])*(p1[i] - p2[i]);
+ return sqrt( (wp-wq)*(wp-wq) / (4 * sq_l2dist) - (wp+wq) / 2 + sq_l2dist / 4 );
+ }
+};
+
+void read_weighted_points(const std::string& infilename, PointContainer& points)
+{
+ std::ifstream in(infilename.c_str());
+ std::string line;
+ while(std::getline(in, line))
+ {
+ if (line[0] == '#') continue; // comment line in the file
+ std::stringstream linestream(line);
+ double x;
+ points.push_back(Point());
+ while (linestream >> x)
+ points.back().push_back(x);
+ }
+}
+
+#endif // __L2_DISTANCE_H__
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/geometry/weighted-l2distance.h Sun Dec 20 08:27:01 2009 -0800
@@ -0,0 +1,53 @@
+#ifndef __L2_DISTANCE_H__
+#define __L2_DISTANCE_H__
+
+#include <utilities/types.h>
+
+#include <vector>
+#include <fstream>
+#include <functional>
+#include <cmath>
+#include <string>
+#include <sstream>
+
+
+typedef std::vector<double> Point;
+typedef std::vector<Point> PointContainer;
+
+struct WeightedL2Distance:
+ public std::binary_function<const Point&, const Point&, double>
+{
+ result_type operator()(const Point& p1, const Point& p2) const
+ {
+ AssertMsg(p1.size() == p2.size(), "Points must be in the same dimension (in L2Distance): dim1=%d, dim2=%d", p1.size()-1, p2.size()-1);
+
+ /* the distance of a point to itself is the radius at which the "power distance" ball around it appears:
+ d(p,p) := sqrt(-w_p) */
+ if (p1 == p2)
+ return sqrt(-p1[p1.size()-1]);
+
+ /* otherwise, the distance is the radius at which the power balls around p, q intersect:
+ d(p,q) := sqrt( ||p-q||^2 - w_p - w_q ) */
+ result_type sq_l2dist = 0;
+ for (size_t i = 0; i < p1.size()-1; ++i)
+ sq_l2dist += (p1[i] - p2[i])*(p1[i] - p2[i]);
+ return sqrt(sq_l2dist - p1[p1.size()-1] - p2[p2.size()-1]);
+ }
+};
+
+void read_weighted_points(const std::string& infilename, PointContainer& points)
+{
+ std::ifstream in(infilename.c_str());
+ std::string line;
+ while(std::getline(in, line))
+ {
+ if (line[0] == '#') continue; // comment line in the file
+ std::stringstream linestream(line);
+ double x;
+ points.push_back(Point());
+ while (linestream >> x)
+ points.back().push_back(x);
+ }
+}
+
+#endif // __L2_DISTANCE_H__
--- a/include/topology/chain.hpp Sun Dec 20 07:45:29 2009 -0800
+++ b/include/topology/chain.hpp Sun Dec 20 08:27:01 2009 -0800
@@ -155,6 +155,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.h Sun Dec 20 07:45:29 2009 -0800
+++ b/include/topology/cohomology-persistence.h Sun Dec 20 08:27:01 2009 -0800
@@ -36,7 +36,7 @@
CohomologyPersistence(const Field& field = Field()):
- field_(field) {}
+ field_(field), image_begin_(cocycles_.end()) {}
// An entry in a cocycle column
@@ -60,10 +60,14 @@
// return either a SimplexIndex or a Death
// BI = BoundaryIterator; it should dereference to a SimplexIndex
template<class BI>
- IndexDeathPair add(BI begin, BI end, BirthInfo b, bool store = true, const SimplexData& sd = SimplexData());
+ IndexDeathPair add(BI begin, BI end, BirthInfo b, bool store = true, const SimplexData& sd = SimplexData(), bool image = true);
+
+ // if sign needs to be specified explicitly, provide (parallel) coefficient_iter
+ template<class BI, class CI>
+ IndexDeathPair add(CI coefficient_iter, BI begin, BI end, BirthInfo b, bool store = true, const SimplexData& sd = SimplexData(), bool image = true);
void show_cocycles() const;
- CocycleIndex begin() { return cocycles_.begin(); }
+ CocycleIndex begin() { return image_begin_; }
CocycleIndex end() { return cocycles_.end(); }
private:
@@ -72,6 +76,7 @@
private:
Simplices simplices_;
Cocycles cocycles_;
+ CocycleIndex image_begin_;
Field field_;
};
@@ -117,7 +122,7 @@
ZColumn zcolumn;
BirthInfo birth;
- unsigned order;
+ signed order;
bool operator<(const Cocycle& other) const { return order > other.order; }
bool operator==(const Cocycle& other) const { return order == other.order; }
--- a/include/topology/cohomology-persistence.hpp Sun Dec 20 07:45:29 2009 -0800
+++ b/include/topology/cohomology-persistence.hpp Sun Dec 20 08:27:01 2009 -0800
@@ -1,26 +1,64 @@
#include <boost/utility.hpp>
#include <queue>
#include <vector>
+#include <limits>
#include <utilities/log.h>
#include <utilities/indirect.h>
+#include <utilities/counter.h>
+
+#include <boost/iterator/transform_iterator.hpp>
+#include <boost/iterator/counting_iterator.hpp>
+#include <boost/foreach.hpp>
+
+#include <boost/lambda/lambda.hpp>
+#include <boost/lambda/bind.hpp>
+namespace bl = boost::lambda;
#ifdef LOGGING
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
{
public:
bool operator()(const SNode& s1, const SNode& s2) const { return s1.si->order < s2.si->order; }
};
+
+
+struct Alternator
+{
+ typedef int result_type;
+ int operator()(int x) const { return (x % 2)*2 - 1; }
+};
+
template<class BirthInfo, class SimplexData, class Field>
template<class BI>
typename CohomologyPersistence<BirthInfo, SimplexData, Field>::IndexDeathPair
CohomologyPersistence<BirthInfo, SimplexData, Field>::
-add(BI begin, BI end, BirthInfo birth, bool store, const SimplexData& sd)
+add(BI begin, BI end, BirthInfo birth, bool store, const SimplexData& sd, bool image)
+{
+ // Set coefficient to be an iterator over (-1)^i
+
+ return add(boost::make_transform_iterator(boost::make_counting_iterator(1), Alternator()),
+ begin, end,
+ birth, store, sd, image);
+}
+
+
+template<class BirthInfo, class SimplexData, class Field>
+template<class BI, class CI>
+typename CohomologyPersistence<BirthInfo, SimplexData, Field>::IndexDeathPair
+CohomologyPersistence<BirthInfo, SimplexData, Field>::
+add(CI coefficient_iter, BI begin, BI end, BirthInfo birth, bool store, const SimplexData& sd, bool image)
{
// Create simplex representation
simplices_.push_back(SHead(sd, simplices_.empty() ? 0 : (simplices_.back().order + 1)));
@@ -30,14 +68,15 @@
typedef std::list<CocycleCoefficientPair> Candidates;
Candidates candidates, candidates_bulk;
rLog(rlCohomology, "Boundary");
- bool sign = true; // TODO: this is very crude, fix later
+
for (BI cur = begin; cur != end; ++cur)
{
- rLog(rlCohomology, " %d %d", (*cur)->order, sign);
- for (typename ZRow::const_iterator zcur = (*cur)->row.begin(); zcur != (*cur)->row.end(); ++zcur)
- candidates_bulk.push_back(std::make_pair(zcur->ci,
- (sign ? zcur->coefficient : field_.neg(zcur->coefficient))));
- sign = !sign;
+ FieldElement coefficient = field_.init(*coefficient_iter++);
+ SimplexIndex cursi = *cur;
+
+ 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)));
}
candidates_bulk.sort(make_first_comparison(make_indirect_comparison(std::less<Cocycle>())));
@@ -73,22 +112,43 @@
// Birth
if (candidates.empty())
{
+ // rLog(rlCohomology, " Birth occurred");
if (!store)
{
simplices_.pop_back();
- return std::make_pair(simplices_.begin(), Death()); // TODO: shouldn't return front
+ return std::make_pair(simplices_.end(), Death());
}
- unsigned order = cocycles_.empty() ? 0 : cocycles_.front().order + 1;
- cocycles_.push_front(Cocycle(birth, order));
+ signed order = 0;
+ if (image)
+ if (image_begin_ == cocycles_.end())
+ order = std::numeric_limits<signed>::min();
+ else
+ order = image_begin_->order + 1;
+ else
+ if (!cocycles_.empty() && cocycles_.front().order >= 0) // we have something outside the image
+ order = cocycles_.front().order + 1;
+
+ CocycleIndex nw;
+ if (image)
+ {
+ image_begin_ = cocycles_.insert(image_begin_, Cocycle(birth, order));
+ nw = image_begin_;
+ } else
+ {
+ cocycles_.push_front(Cocycle(birth, order));
+ nw = cocycles_.begin();
+ }
- rLog(rlCohomology, "Birth: %d", cocycles_.front().order);
+ rLog(rlCohomology, "Birth: %d", nw->order);
// set up the cocycle
- ZColumn& cocycle = cocycles_.front().zcolumn;
- cocycle.push_back(SNode(si, field_.id(), cocycles_.begin()));
- si->row.push_back(cocycles_.front().zcolumn.front());
+ ZColumn& cocycle = nw->zcolumn;
+ 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());
}
@@ -105,15 +165,26 @@
CocycleCoefficientPair& z = candidates.front();
Death d = z.first->birth;
+ rLog(rlCohomology, " Order: %d", z.first->order);
+ if (z.first->order >= 0) // if death outside image
+ d = Death(); // no death occurs outside the image
+ else
+ if (z.first == image_begin_)
+ ++image_begin_;
// 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);
@@ -151,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);
@@ -173,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 Sun Dec 20 07:45:29 2009 -0800
+++ b/include/topology/dynamic-persistence.h Sun Dec 20 08:27:01 2009 -0800
@@ -4,8 +4,6 @@
#include "static-persistence.h"
#include <utilities/types.h>
-#include <boost/progress.hpp>
-
#include <boost/bind.hpp>
#include <boost/lambda/lambda.hpp>
namespace bl = boost::lambda;
@@ -144,11 +142,12 @@
struct PairingTrailsVisitor: public Parent::PairVisitor
{
- PairingTrailsVisitor(Order& order, ConsistencyComparison ccmp):
- order_(order), ccmp_(ccmp) {}
+ PairingTrailsVisitor(Order& order, ConsistencyComparison ccmp, unsigned size):
+ Parent::PairVisitor(size), order_(order), ccmp_(ccmp) {}
void init(iterator i) const { order_.modify(i, boost::bind(&Element::template trail_append<ConsistencyComparison>, bl::_1, &*i, ccmp_)); Count(cTrailLength); } // i->trail_append(&*i, ccmp)
void update(iterator j, iterator i) const { order_.modify(order_.iterator_to(*(i->pair)), boost::bind(&Element::template trail_append<ConsistencyComparison>, bl::_1, &*j, ccmp_)); Count(cTrailLength); } // i->pair->trail_append(&*j, ccmp)
+ void finished(iterator i) const { Parent::PairVisitor::finished(i); }
Order& order_;
ConsistencyComparison ccmp_;
@@ -288,17 +287,14 @@
struct PairingChainsVisitor: public Parent::PairVisitor
{
PairingChainsVisitor(Order& order, ConsistencyComparison ccmp, unsigned size):
- order_(order), ccmp_(ccmp), show_progress(size) {}
+ Parent::PairVisitor(size), order_(order), ccmp_(ccmp) {}
void init(iterator i) const { order_.modify(i, boost::bind(&Element::template chain_append<ConsistencyComparison>, bl::_1, &*i, ccmp_)); } // i->chain_append(&*i, ccmp)
void update(iterator j, iterator i) const { order_.modify(j, boost::bind(&Element::template chain_add<ConsistencyComparison>, bl::_1, i->pair->chain, ccmp_)); } // j->chain.add(i->pair->chain, ccmp_)
- void finished(iterator i) const { CountBy(cChainLength, i->chain.size()); ++show_progress; }
+ void finished(iterator i) const { Parent::PairVisitor::finished(i); CountBy(cChainLength, i->chain.size()); }
Order& order_;
ConsistencyComparison ccmp_;
-
- mutable boost::progress_display
- show_progress;
};
ConsistencyComparison ccmp_;
--- a/include/topology/dynamic-persistence.hpp Sun Dec 20 07:45:29 2009 -0800
+++ b/include/topology/dynamic-persistence.hpp Sun Dec 20 08:27:01 2009 -0800
@@ -34,7 +34,7 @@
DynamicPersistenceTrails<D,CT,OT,E,Cmp,CCmp>::
pair_simplices()
{
- Parent::pair_simplices(begin(), end(), PairingTrailsVisitor(order(), ccmp_));
+ Parent::pair_simplices(begin(), end(), true, PairingTrailsVisitor(order(), ccmp_, size()));
}
template<class D, class CT, class OT, class E, class Cmp, class CCmp>
@@ -279,5 +279,5 @@
DynamicPersistenceChains<D,CT,OT,E,Cmp,CCmp>::
pair_simplices()
{
- Parent::pair_simplices(begin(), end(), PairingChainsVisitor(order(), ccmp_, size()));
+ Parent::pair_simplices(begin(), end(), true, PairingChainsVisitor(order(), ccmp_, size()));
}
--- a/include/topology/field-arithmetic.h Sun Dec 20 07:45:29 2009 -0800
+++ b/include/topology/field-arithmetic.h Sun Dec 20 08:27:01 2009 -0800
@@ -13,6 +13,7 @@
Element id() const { return 1; }
Element zero() const { return 0; }
+ Element init(int a) const { return (a % p_ + p_) % p_; }
Element neg(Element a) const { return p_ - a; }
Element add(Element a, Element b) const { return (a+b) % p_; }
--- a/include/topology/rips.h Sun Dec 20 07:45:29 2009 -0800
+++ b/include/topology/rips.h Sun Dec 20 08:27:01 2009 -0800
@@ -84,7 +84,7 @@
DistanceType distance(const Simplex& s1, const Simplex& s2) const;
- private:
+ protected:
class WithinDistance;
template<class Functor, class NeighborTest>
@@ -96,7 +96,7 @@
const Functor& functor,
bool check_initial = true) const;
- private:
+ protected:
const Distances& distances_;
};
@@ -111,7 +111,7 @@
bool operator()(Vertex u, Vertex v) const { return distances_(u, v) <= max_; }
- private:
+ protected:
const Distances& distances_;
DistanceType max_;
};
@@ -127,7 +127,7 @@
DistanceType operator()(const Simplex& s) const;
- private:
+ protected:
const Distances& distances_;
};
@@ -150,7 +150,7 @@
return e1 < e2;
}
- private:
+ protected:
Evaluator eval_;
};
--- a/include/topology/rips.hpp Sun Dec 20 07:45:29 2009 -0800
+++ b/include/topology/rips.hpp Sun Dec 20 08:27:01 2009 -0800
@@ -136,7 +136,7 @@
functor(s);
}
- if (current.size() == max_dim + 1)
+ if (current.size() == static_cast<size_t>(max_dim) + 1)
return;
rLog(rlRipsDebug, "Traversing %d vertices", candidates.end() - boost::next(excluded));
--- a/include/topology/simplex.h Sun Dec 20 07:45:29 2009 -0800
+++ b/include/topology/simplex.h Sun Dec 20 08:27:01 2009 -0800
@@ -127,12 +127,12 @@
struct DataEvaluator;
struct DimensionExtractor;
- private:
+ protected:
VertexContainer& vertices() { return vdpair_.first(); }
VerticesDataPair vdpair_;
- private:
+ protected:
/* Serialization */
friend class boost::serialization::access;
--- a/include/topology/static-persistence.h Sun Dec 20 07:45:29 2009 -0800
+++ b/include/topology/static-persistence.h Sun Dec 20 08:27:01 2009 -0800
@@ -11,6 +11,8 @@
#include <utilities/types.h>
+#include <boost/progress.hpp>
+
// Element_ should derive from PairCycleData
template<class Data_, class ChainTraits_, class Element_ = use_default>
@@ -88,7 +90,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(), false, PairVisitor(size())); }
// Functions: Accessors
// begin() - returns OrderIndex of the first element
@@ -122,12 +124,13 @@
// Function: pair_simplices(bg, end)
// Compute persistence of the simplices in filtration between bg and end
template<class Visitor>
- void pair_simplices(iterator bg, iterator end, const Visitor& visitor = Visitor());
+ void pair_simplices(iterator bg, iterator end, bool store_negative = false, const Visitor& visitor = Visitor());
// Struct: PairVisitor
// 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`)
@@ -141,7 +144,9 @@
// Function: finished(j)
// Called after the processing of `j` is finished.
- void finished(iterator j) const {}
+ void finished(iterator j) const { ++show_progress; }
+ mutable boost::progress_display
+ show_progress;
};
const Order& order() const { return order_; }
--- a/include/topology/static-persistence.hpp Sun Dec 20 07:45:29 2009 -0800
+++ b/include/topology/static-persistence.hpp Sun Dec 20 08:27:01 2009 -0800
@@ -6,6 +6,8 @@
#include <boost/utility.hpp>
#include <boost/foreach.hpp>
+#include <boost/foreach.hpp>
+
#ifdef LOGGING
static rlog::RLogChannel* rlPersistence = DEF_CHANNEL("topology/persistence", rlog::Log_Debug);
#endif // LOGGING
@@ -43,7 +45,7 @@
template<class Visitor>
void
StaticPersistence<D, CT, OT, E, Cmp>::
-pair_simplices(iterator bg, iterator end, const Visitor& visitor)
+pair_simplices(iterator bg, iterator end, bool store_negative, const Visitor& visitor)
{
#if LOGGING
typename ContainerTraits::OutputMap outmap(order_);
@@ -59,6 +61,17 @@
Cycle z;
swap_cycle(j, z);
rLog(rlPersistence, " has boundary: %s", z.tostring(outmap).c_str());
+
+ // Sparsify the cycle by removing the negative elements
+ if (!store_negative)
+ {
+ typename OrderElement::Cycle zz;
+ BOOST_FOREACH(OrderIndex i, z)
+ if (i->sign()) // positive
+ zz.push_back(i);
+ z.swap(zz);
+ }
+ // --------------------------
CountNum(cPersistencePairBoundaries, z.size());
Count(cPersistencePair);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/include/topology/weighted-rips.h Sun Dec 20 08:27:01 2009 -0800
@@ -0,0 +1,139 @@
+#ifndef __WEIGHTED_RIPS_H__
+#define __WEIGHTED_RIPS_H__
+
+#include <vector>
+#include <string>
+#include "simplex.h"
+#include "rips.h"
+#include <boost/iterator/counting_iterator.hpp>
+
+/**
+ * WeightedRipsSimplex class
+ *
+ * This class sits as an invisible layer between the Simplex datatype passed
+ * to WeightedRips and the class itself. The need for this layer is the need
+ * to store the ``value'' (max inter-vertex distance) of each simplex in the
+ * Weighted Rips complex--something that the user of the class does not need
+ * to be aware of.
+ */
+
+template<class Simplex_, class DistanceType_>
+class WeightedRipsSimplex : public Simplex_
+{
+ public:
+ typedef typename Simplex_::Vertex Vertex;
+ typedef typename Simplex_::VertexContainer VertexContainer;
+ typedef DistanceType_ DistanceType;
+
+ WeightedRipsSimplex(Simplex_ s) : Simplex_(s) { }
+
+ void setSimplexValue(const DistanceType &sv) { simplexValue = sv; }
+ DistanceType getSimplexValue() const { return simplexValue; }
+
+ protected:
+ DistanceType simplexValue;
+};
+
+/**
+ * WeightedRips class
+ *
+ * Class providing basic operations to work with Rips complexes. It implements Bron-Kerbosch algorithm,
+ * and provides simple wrappers for various functions.
+ *
+ * Distances_ is expected to define types IndexType and DistanceType as well as
+ * provide operator()(...) which given two IndexTypes should return
+ * the distance between them. There should be methods begin() and end()
+ * for iterating over IndexTypes as well as a method size().
+ */
+template<class Distances_, class Simplex_ = Simplex<typename Distances_::IndexType> >
+class WeightedRips : public Rips<Distances_, Simplex_>
+{
+ public:
+
+ /* redeclaring the typedefs because they cannot be inherited at compile-time */
+ typedef Distances_ Distances;
+ typedef typename Distances::IndexType IndexType;
+ typedef typename Distances::DistanceType DistanceType;
+
+ typedef WeightedRipsSimplex<Simplex_, DistanceType> Simplex;
+ typedef typename Simplex::Vertex Vertex; // should be the same as IndexType
+ typedef typename Simplex::VertexContainer VertexContainer;
+
+ class Evaluator;
+ class Comparison;
+
+ public:
+ WeightedRips(const Distances& distances):
+ Rips<Distances_, Simplex_>(distances) {}
+
+ template<class Functor>
+ void generate(Dimension k, DistanceType max, const Functor& f) const;
+
+};
+
+/**
+ * DistanceDataStackingFunctor class
+ *
+ * Class providing a functor that is to be called by WeightedRips::generate(). This functor
+ * takes as an argument (to its constructor) the original functor passed by the user to
+ * generate(), and a new ``double'' functor is created. Assuming that the functor acts on
+ * simplices, first the value of the simplex is computed (the radius at which the simplex
+ * appears in the weighted Rips complex), the data field of the simplex is populated with
+ * this value, and then the original functor is called (it has no idea that it was
+ * intercepted).
+ */
+
+template<class Rips_, class Functor_>
+class DistanceDataStackingFunctor
+{
+ public:
+ typedef typename Rips_::Simplex Simplex_;
+
+ DistanceDataStackingFunctor(const Rips_ &r, const Functor_ &f):
+ rips(r), original_functor(f) { }
+
+ void operator()(const Simplex_ &s) const
+ {
+ Simplex_ s_new(s);
+ s_new.setSimplexValue (rips.distance(s_new, s_new));
+ original_functor (s_new);
+ }
+
+ private:
+ const Rips_ &rips;
+ const Functor_ &original_functor;
+};
+
+template<class Distances_, class Simplex_>
+template<class Functor>
+void WeightedRips<Distances_, Simplex_>::generate(Dimension k, DistanceType max, const Functor &f) const
+{
+ Rips<Distances_,Simplex_>::generate(k, max, DistanceDataStackingFunctor<WeightedRips<Distances_, Simplex_>,Functor>(*this, f));
+}
+
+template<class Distances_, class Simplex_>
+class WeightedRips<Distances_, Simplex_>::Evaluator: public Rips<Distances_,Simplex_>::Evaluator
+{
+ public:
+ Evaluator(const Distances& distances):
+ Rips<Distances_, Simplex_>::Evaluator(distances) {}
+
+ DistanceType operator()(const Simplex& s) const { return s.getSimplexValue(); }
+};
+
+template<class Distances_, class Simplex_>
+class WeightedRips<Distances_, Simplex_>::Comparison: public Rips<Distances_,Simplex_>::Comparison
+{
+ public:
+ Comparison(const Distances& distances):
+ Rips<Distances_, Simplex_>::Comparison(distances) {}
+
+ bool operator()(const Simplex& s1, const Simplex& s2) const
+ {
+ if (s1.dimension() != s2.dimension())
+ return s1.dimension() < s2.dimension();
+ return s1.getSimplexValue() < s2.getSimplexValue();
+ }
+};
+
+#endif // __WEIGHTED_RIPS_H__
--- a/include/topology/zigzag-persistence.h Sun Dec 20 07:45:29 2009 -0800
+++ b/include/topology/zigzag-persistence.h Sun Dec 20 08:27:01 2009 -0800
@@ -127,7 +127,7 @@
ZIndex new_z_in_remove(ZigzagPersistence& zz);
- void erasing_z(ZigzagPersistence& zz, ZIndex j) {}
+ void erasing_z(ZigzagPersistence&, ZIndex) {}
Death death(ZigzagPersistence& zz, ZIndex dying_z);
};
--- a/include/topology/zigzag-persistence.hpp Sun Dec 20 07:45:29 2009 -0800
+++ b/include/topology/zigzag-persistence.hpp Sun Dec 20 08:27:01 2009 -0800
@@ -414,7 +414,7 @@
template<class BID, class SD>
bool
ZigzagPersistence<BID,SD>::
-check_consistency(SimplexIndex s_skip, ZIndex z_skip, BIndex b_skip)
+check_consistency(SimplexIndex, ZIndex, BIndex)
{
#ifdef ZIGZAG_CONSISTENCY
#warning "Checking consistency in ZigzagPersistence"
@@ -684,7 +684,7 @@
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)
+new_z_in_add(ZigzagPersistence& zz, const ZColumn&, const BRow&)
{
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()));
@@ -694,7 +694,7 @@
template<class BID, class SD>
typename ZigzagPersistence<BID,SD>::BIndex
ZigzagPersistence<BID,SD>::ZigzagVisitor::
-select_j_in_remove(ZigzagPersistence& zz, const CRow& c_row)
+select_j_in_remove(ZigzagPersistence&, const CRow& c_row)
{
return c_row.front();
}
@@ -712,7 +712,7 @@
template<class BID, class SD>
typename ZigzagPersistence<BID,SD>::Death
ZigzagPersistence<BID,SD>::ZigzagVisitor::
-death(ZigzagPersistence& zz, ZIndex dying_z)
+death(ZigzagPersistence&, ZIndex dying_z)
{
return Death(dying_z->birth);
}
--- a/include/utilities/containers.h Sun Dec 20 07:45:29 2009 -0800
+++ b/include/utilities/containers.h Sun Dec 20 08:27:01 2009 -0800
@@ -188,16 +188,17 @@
typedef vector<T> Container;
typedef SizeStorage<Container> Self;
- SizeStorage(size_t size = 0) {}
+ SizeStorage() {}
+ SizeStorage(size_t) {}
- Self& operator+=(size_t inc) { return *this; }
- Self& operator-=(size_t dec) { return *this; }
+ Self& operator+=(size_t) { return *this; }
+ Self& operator-=(size_t) { 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 swap(SizeStorage&) {}
void clear() {}
};
@@ -209,16 +210,17 @@
typedef deque<T> Container;
typedef SizeStorage<Container> Self;
- SizeStorage(size_t size = 0) {}
+ SizeStorage() {}
+ SizeStorage(size_t) {}
- Self& operator+=(size_t inc) { return *this; }
- Self& operator-=(size_t dec) { return *this; }
+ Self& operator+=(size_t) { return *this; }
+ Self& operator-=(size_t) { 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 swap(SizeStorage&) {}
void clear() {}
};
--- a/include/utilities/counter.h Sun Dec 20 07:45:29 2009 -0800
+++ b/include/utilities/counter.h Sun Dec 20 08:27:01 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 Sun Dec 20 07:45:29 2009 -0800
+++ b/include/utilities/memory.h Sun Dec 20 08:27:01 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
--- a/tools/CMakeLists.txt Sun Dec 20 07:45:29 2009 -0800
+++ b/tools/CMakeLists.txt Sun Dec 20 08:27:01 2009 -0800
@@ -1,15 +1,19 @@
-find_package (Qt4 REQUIRED)
-set (QT_USE_QTOPENGL TRUE)
-set (QT_USE_QTXML TRUE)
-include (${QT_USE_FILE})
+find_package (Qt4)
+if (QT4_FOUND)
+ set (QT_USE_QTOPENGL TRUE)
+ set (QT_USE_QTXML TRUE)
+ include (${QT_USE_FILE})
+
+ add_subdirectory (diagram-viewer)
#find_library (gle_LIBRARY NAMES gle)
#find_library (QGLViewer_LIBRARY NAMES QGLViewer)
#find_path (QGLViewer_INCLUDE_DIR QGLViewer/qglviewer.h)
#include_directories (${QGLViewer_INCLUDE_DIR})
+endif (QT4_FOUND)
+
add_executable (extract-diagram extract-diagram.cpp)
target_link_libraries (extract-diagram ${libraries} ${Boost_SERIALIZATION_LIBRARY})
-add_subdirectory (diagram-viewer)
add_subdirectory (matching)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/draw-diagram/draw.py Sun Dec 20 08:27:01 2009 -0800
@@ -0,0 +1,25 @@
+#!/usr/bin/env python
+
+import pd
+from sys import argv, exit
+
+if len(argv) < 2:
+ print "Usage: %s FILENAME [MULTIPLIER=1] [NOISE=0] [RADIUS=.15] [DIMENSIONS=XMIN,YMIN,XMAX,YMAX]" % argv[0]
+ print " MULTIPLIER - multiply coordinates of each point by this quantity"
+ print " NOISE - filter out points below this persistence"
+ print " RADIUS - radius of a point in the persistence diagram"
+ print " DIMENSIONS - dimensions of the persistence diagram"
+ print
+ print " Example: %s torus.dgm 1 0 .05 -1,-1,10,10" % argv[0]
+ exit()
+
+multiplier = float(argv[2]) if len(argv) > 2 else 1
+noise = float(argv[3]) if len(argv) > 3 else 0
+R = float(argv[4]) if len(argv) > 4 else .15
+dimensions = map(float, argv[5].split(',')) if len(argv) > 5 else None
+
+noise_filter = pd.noise_filter(noise)
+amplify_filter = pd.amplify_filter(multiplier)
+
+dgm = pd.PersistenceDiagram(argv[1], lambda x,y: noise_filter(x,y) and amplify_filter(x,y))
+dgm.savePDF(argv[1] + '.', radius = R, dimensions = dimensions)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/draw-diagram/pd.py Sun Dec 20 08:27:01 2009 -0800
@@ -0,0 +1,113 @@
+#!/usr/bin/python
+
+# Author: Dmitriy Morozov <morozov@cs.duke.edu>
+# (2004-2008) Department of Computer Science, Duke University
+
+import pyx, string, sys, math
+
+class PersistenceDiagram:
+ def drawAxes(self,c,radius,spacing = 1, dimensions = None):
+ if not dimensions:
+ xmax = math.ceil(max(self.xmax,1) + radius + 1)
+ xmin = math.floor(min(self.xmin,-1) - radius - 1)
+ ymax = math.ceil(max(self.ymax,1) + radius + 1)
+ ymin = math.floor(min(self.ymin,-1) - radius - 1)
+ else:
+ xmin,ymin,xmax,ymax = dimensions
+
+ minmax = min([xmax, ymax])
+ maxmin = max([xmin, ymin])
+ for i in xrange(int(math.floor(xmin)),int(math.ceil(xmax)), spacing):
+ c.stroke(pyx.path.line(i,ymin,i,ymax), [pyx.style.linestyle.dashed, pyx.color.gray(0.5)])
+ for j in xrange(int(math.floor(ymin)),int(math.ceil(ymax)), spacing):
+ c.stroke(pyx.path.line(xmin,j,xmax,j), [pyx.style.linestyle.dashed, pyx.color.gray(0.5)])
+ c.stroke(pyx.path.line(xmin,0,xmax,0), [pyx.deco.earrow.normal])
+ c.stroke(pyx.path.line(0,ymin,0,ymax), [pyx.deco.earrow.normal])
+ c.stroke(pyx.path.line(maxmin,maxmin,minmax,minmax))
+
+
+ def drawCanvas(self, c, points, color = 'red', filled = 1, radius = 0.15):
+ for p in points:
+ self.drawPoint(c,p,color,filled,radius)
+
+ def drawPoint(self,c,p, color = 'red', filled = 1, radius = 0.15):
+ if color == 'red':
+ options = [pyx.color.rgb.red]
+ elif color == 'blue':
+ options = [pyx.color.rgb.blue]
+ elif color == 'green':
+ options = [pyx.color.rgb.green]
+ else:
+ options = []
+
+ if filled:
+ draw = c.fill
+ else:
+ draw = c.stroke
+
+ xmax = max(self.xmax,1)
+ xmin = min(self.xmin,-1)
+ ymax = max(self.ymax,1)
+ ymin = min(self.ymin,-1)
+
+ x,y = p
+ if abs(x) == float('inf') or abs(y) == float('inf'): radius = radius * 2
+ if x == float('inf'): x = xmax + 1
+ if x == float('-inf'): x = xmin - 1
+ if y == float('inf'): y = ymax + 1
+ if y == float('-inf'): y = ymin -1
+
+ draw(pyx.path.circle(x,y,radius), options)
+
+ def savePDF(self, filename, color = 'red', filled = 1, radius = 0.15, axes = 1, dimensions = None):
+ for d in self.points.keys():
+ c = pyx.canvas.canvas()
+ self.drawAxes(c, radius, axes, dimensions)
+ self.drawCanvas(c, self.points[d], color, filled, radius)
+ c.writePDFfile(filename + str(d))
+
+ def add(self, d, p, filter):
+ x,y = p
+ p = filter(x,y)
+ if not p: return
+ x,y = p
+
+ if d not in self.points.keys():
+ self.points[d] = []
+ self.points[d] += [(x,y)]
+
+ if abs(x) != float('inf'):
+ self.xmax = max(x,self.xmax)
+ self.xmin = min(x,self.xmin)
+ if abs(y) != float('inf'):
+ self.ymax = max(y,self.ymax)
+ self.ymin = min(y,self.ymin)
+
+ def load(self,filename, filter):
+ self.xmax = self.ymax = 0
+ self.xmin = self.ymin = 0
+ f = file(filename, 'r')
+ for line in f:
+ if line.strip().startswith('#'): continue
+ dim,xstr,ystr = string.split(line)
+ self.add(dim, (float(xstr),float(ystr)), filter)
+
+ def __init__(self, filename, filter = lambda x,y: (x,y)):
+ self.points = {}
+ self.load(filename, filter)
+
+
+def noise_filter(epsilon):
+ def noise(x,y):
+ if y - x <= epsilon: return None
+ return (x,y)
+
+ return noise
+
+def amplify_filter(x_mult, y_mult = None):
+ if not y_mult: y_mult = x_mult
+
+ def amplify(x,y):
+ return (x_mult*x, y_mult*y)
+
+ return amplify