summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorBryn Keller <xoltar@xoltar.org>2016-03-07 13:51:48 -0800
committerBryn Keller <xoltar@xoltar.org>2016-03-07 13:51:48 -0800
commitd7677b1f3f5ceee9cfdb5275ef00cbf79a714122 (patch)
treefb6e9356e5093b0e077979faacf9de6f215152c2
parent41a3fef94446ec6ce115c65d92890ae9b888690c (diff)
Python interface to PHAT
-rw-r--r--.gitignore13
-rw-r--r--README.md4
-rw-r--r--include/phat/boundary_matrix.h8
-rw-r--r--python/.gitignore7
-rw-r--r--python/MANIFEST.in1
-rw-r--r--python/README.rst110
-rw-r--r--python/_phat.cpp279
-rw-r--r--python/phat.py88
-rw-r--r--python/setup.cfg5
-rw-r--r--python/setup.py57
-rw-r--r--python/src/self_test.py168
-rw-r--r--python/src/simple_example.py59
-rw-r--r--src/simple_example.cpp2
13 files changed, 796 insertions, 5 deletions
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..e05a4ac
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,13 @@
+GPATH
+GRTAGS
+GTAGS
+CMakeCache.txt
+self_test
+simple_example
+phat
+info
+convert
+cmake_install.cmake
+benchmark
+Makefile
+CMakeFiles \ No newline at end of file
diff --git a/README.md b/README.md
index b953afa..9a30ffa 100644
--- a/README.md
+++ b/README.md
@@ -7,7 +7,7 @@ Ulrich Bauer, Michael Kerber, Jan Reininghaus
## Contributors: ##
-Hubert Wagner
+Hubert Wagner, Bryn Keller
## Downloads: ##
* [PHAT, v1.4.0](https://drive.google.com/uc?id=0B7Yz6TPEpiGEWmNyeVFsNXgtUGc&export=download)
@@ -138,4 +138,4 @@ References:
1. H.Edelsbrunner, J.Harer: Computational Topology, An Introduction. American Mathematical Society, 2010, ISBN 0-8218-4925-5
2. V.de Silva, D.Morozov, M.Vejdemo-Johansson: Dualities in persistent (co)homology. Inverse Problems 27, 2011
3. C.Chen, M.Kerber: Persistent Homology Computation With a Twist. 27th European Workshop on Computational Geometry, 2011.
-4. U.Bauer, M.Kerber, J.Reininghaus: Clear and Compress: Computing Persistent Homology in Chunks. [http://arxiv.org/pdf/1303.0477.pdf arXiv:1303.0477] \ No newline at end of file
+4. U.Bauer, M.Kerber, J.Reininghaus: Clear and Compress: Computing Persistent Homology in Chunks. [http://arxiv.org/pdf/1303.0477.pdf arXiv:1303.0477]
diff --git a/include/phat/boundary_matrix.h b/include/phat/boundary_matrix.h
index 10c66cc..f864dee 100644
--- a/include/phat/boundary_matrix.h
+++ b/include/phat/boundary_matrix.h
@@ -68,7 +68,7 @@ namespace phat {
// finalizes given column
void finalize( index idx ) { rep._finalize( idx ); }
- // syncronizes all internal data structures -- has to be called before and after any multithreaded access!
+ // synchronizes all internal data structures -- has to be called before and after any multithreaded access!
void sync() { rep._sync(); }
// info functions -- independent of chosen 'Representation'
@@ -283,12 +283,16 @@ namespace phat {
// Loads boundary_matrix from given file
// Format: nr_columns % dim1 % N1 % row1 row2 % ...% rowN1 % dim2 % N2 % ...
- bool load_binary( std::string filename )
+ bool load_binary(const std::string &filename )
{
std::ifstream input_stream( filename.c_str( ), std::ios_base::binary | std::ios_base::in );
if( input_stream.fail( ) )
return false;
+ //TODO Note that if you read a file that isn't actually a data file, you may get
+ //a number of columns that is bigger than the available memory, which leads to crashes
+ //with deeply confusing error messages. Consider ways to prevent this.
+ //Magic number in the file header? Check for more columns than bytes in the file?
int64_t nr_columns;
input_stream.read( (char*)&nr_columns, sizeof( int64_t ) );
this->set_num_cols( (index)nr_columns );
diff --git a/python/.gitignore b/python/.gitignore
new file mode 100644
index 0000000..b576037
--- /dev/null
+++ b/python/.gitignore
@@ -0,0 +1,7 @@
+__pycache__
+*.so
+build
+dist
+old
+phat.egg-info
+pybind11
diff --git a/python/MANIFEST.in b/python/MANIFEST.in
new file mode 100644
index 0000000..dba0deb
--- /dev/null
+++ b/python/MANIFEST.in
@@ -0,0 +1 @@
+recursive-include include *.h \ No newline at end of file
diff --git a/python/README.rst b/python/README.rst
new file mode 100644
index 0000000..b67ff47
--- /dev/null
+++ b/python/README.rst
@@ -0,0 +1,110 @@
+Persistent Homology Algorithm Toolkit (PHAT)
+============================================
+
+This is a Python interface for the `Persistent Homology Algorithm Toolkit`_, a software library
+that contains methods for computing the persistence pairs of a
+filtered cell complex represented by an ordered boundary matrix with Z\ :sub:`2` coefficients.
+
+For an introduction to persistent homology, see the textbook [1]_. This software package
+contains code for several algorithmic variants:
+
+* The "standard" algorithm (see [1]_, p.153)
+* The "row" algorithm from [2]_ (called pHrow in that paper)
+* The "twist" algorithm, as described in [3]_ (default algorithm)
+* The "chunk" algorithm presented in [4]_
+* The "spectral sequence" algorithm (see [1]_, p.166)
+
+All but the standard algorithm exploit the special structure of the boundary matrix
+to take shortcuts in the computation. The chunk and the spectral sequence algorithms
+make use of multiple CPU cores if PHAT is compiled with OpenMP support.
+
+All algorithms are implemented as function objects that manipulate a given
+``boundary_matrix`` (to be defined below) object to reduced form.
+From this reduced form one can then easily extract the persistence pairs.
+Alternatively, one can use the ``compute_persistence_pairs function`` which takes an
+algorithm as a parameter, reduces the given ``boundary_matrix`` and stores the
+resulting pairs in a given ``persistence_pairs`` object.
+
+The ``boundary_matrix`` class takes a "Representation" class as a parameter.
+This representation defines how columns of the matrix are represented and how
+low-level operations (e.g., column additions) are performed. The right choice of the
+representation class can be as important for the performance of the program as choosing
+the algorithm. We provide the following choices of representation classes:
+
+* ``vector_vector``: Each column is represented as a sorted ``std::vector`` of integers, containing the indices of the non-zero entries of the column. The matrix itself is a ``std::vector`` of such columns.
+* ``vector_heap``: Each column is represented as a heapified ``std::vector`` of integers, containing the indices of the non-zero entries of the column. The matrix itself is a ``std::vector`` of such columns.
+* ``vector_set``: Each column is a ``std::set`` of integers, with the same meaning as above. The matrix is stored as a ``std::vector`` of such columns.
+* ``vector_list``: Each column is a sorted ``std::list`` of integers, with the same meaning as above. The matrix is stored as a ``std::vector`` of such columns.
+* ``sparse_pivot_column``: The matrix is stored as in the vector_vector representation. However, when a column is manipulated, it is first converted into a ``std::set``, using an extra data field called the "pivot column". When another column is manipulated later, the pivot column is converted back to the ``std::vector`` representation. This can lead to significant speed improvements when many columns are added to a given pivot column consecutively. In a multicore setup, there is one pivot column per thread.
+* ``heap_pivot_column``: The same idea as in the sparse version. Instead of a ``std::set``, the pivot column is represented by a ``std::priority_queue``.
+* ``full_pivot_column``: The same idea as in the sparse version. However, instead of a ``std::set``, the pivot column is expanded into a bit vector of size n (the dimension of the matrix). To avoid costly initializations, the class remembers which entries have been manipulated for a pivot column and updates only those entries when another column becomes the pivot.
+* ``bit_tree_pivot_column`` (default representation): Similar to the ``full_pivot_column`` but the implementation is more efficient. Internally it is a bit-set with fast iteration over nonzero elements, and fast access to the maximal element.
+
+Sample usage
+------------
+
+From src/simple_example.py in the source::
+
+ print("""
+ we will build an ordered boundary matrix of this simplicial complex consisting of a single triangle:
+
+ 3
+ |\\
+ | \\
+ | \\
+ | \\ 4
+ 5| \\
+ | \\
+ | 6 \\
+ | \\
+ |________\\
+ 0 2 1
+
+ """)
+
+ import phat
+
+ # set the dimension of the cell that each column represents:
+ dimensions = [0, 0, 1, 0, 1, 1, 2]
+
+ # define a boundary matrix with the chosen internal representation
+ boundary_matrix = phat.boundary_matrix(representation = phat.representations.vector_vector)
+
+ # set the respective columns -- the columns entries have to be sorted
+ boundary_matrix.set_dims(dimensions)
+ boundary_matrix.set_col(0, [])
+ boundary_matrix.set_col(1, [])
+ boundary_matrix.set_col(2, [0,1])
+ boundary_matrix.set_col(3, [])
+ boundary_matrix.set_col(4, [1,3])
+ boundary_matrix.set_col(5, [0,3])
+ boundary_matrix.set_col(6, [2,4,5])
+
+ # print some information of the boundary matrix:
+ print()
+ print("the boundary matrix has %d columns:" % boundary_matrix.get_num_cols())
+ for col_idx in range(boundary_matrix.get_num_cols()):
+ s = "column %d represents a cell of dimension %d." % (col_idx, boundary_matrix.get_dim(col_idx))
+ if (not boundary_matrix.is_empty(col_idx)):
+ s = s + " its boundary consists of the cells " + " ".join([str(c) for c in boundary_matrix.get_col(col_idx)])
+ print(s)
+ print("overall, the boundary matrix has %d entries." % boundary_matrix.get_num_entries())
+
+ pairs = phat.compute_persistence_pairs(boundary_matrix)
+
+ pairs.sort()
+
+ print()
+
+ print("there are %d persistence pairs: " % len(pairs))
+ for pair in pairs:
+ print("birth: %d, death: %d" % pair)
+
+References:
+
+.. [1] H.Edelsbrunner, J.Harer: Computational Topology, An Introduction. American Mathematical Society, 2010, ISBN 0-8218-4925-5
+.. [2] V.de Silva, D.Morozov, M.Vejdemo-Johansson: Dualities in persistent (co)homology. Inverse Problems 27, 2011
+.. [3] C.Chen, M.Kerber: Persistent Homology Computation With a Twist. 27th European Workshop on Computational Geometry, 2011.
+.. [4] U.Bauer, M.Kerber, J.Reininghaus: Clear and Compress: Computing Persistent Homology in Chunks. arXiv:1303.0477_
+.. _arXiv:1303.0477: http://arxiv.org/pdf/1303.0477.pdf
+.. _`Persistent Homology Algorithm Toolkit`: https://bitbucket.org/phat/phat-code
diff --git a/python/_phat.cpp b/python/_phat.cpp
new file mode 100644
index 0000000..c43a3f5
--- /dev/null
+++ b/python/_phat.cpp
@@ -0,0 +1,279 @@
+//Required header for using pybind11
+#include <pybind11/pybind11.h>
+
+//Automatic conversions of stl containers to Python ones
+#include <pybind11/stl.h>
+
+//Additional support for operators and numpy
+#include <pybind11/operators.h>
+#include <pybind11/numpy.h>
+
+//All the things we're going to wrap
+#include "phat/persistence_pairs.h"
+#include "phat/compute_persistence_pairs.h"
+#include "phat/boundary_matrix.h"
+#include "phat/representations/abstract_pivot_column.h"
+#include <phat/representations/vector_vector.h>
+#include <phat/representations/vector_heap.h>
+#include <phat/representations/vector_set.h>
+#include <phat/representations/vector_list.h>
+#include <phat/representations/sparse_pivot_column.h>
+#include <phat/representations/heap_pivot_column.h>
+#include <phat/representations/full_pivot_column.h>
+#include <phat/representations/bit_tree_pivot_column.h>
+#include <phat/algorithms/twist_reduction.h>
+#include <phat/algorithms/standard_reduction.h>
+#include <phat/algorithms/row_reduction.h>
+#include <phat/algorithms/chunk_reduction.h>
+#include <phat/algorithms/spectral_sequence_reduction.h>
+
+namespace py = pybind11;
+
+//## Some template functions we'll need later
+
+// This function defines two Python functions in the extension module, that are named
+// `compute_persistence_pairs_${rep}_${reduction}`
+// `compute_persistence_pairs_dualized_${rep}_${reductionx}`.
+//
+// The Python user will never see these, since we will use (in phat.py) the type of the
+// boundary matrix and the requested reduction to dispatch to the correct function
+// required.
+//
+// These functions are the main operations of PHAT. In the Python version, they take
+// a boundary matrix, and return a persistence_pairs object.
+//
+// `Reduction` will be an algorithm, `Representation` is a type that controls
+// how the boundary matrix stores its internal state.
+//
+// We will be using this function to define these two functions for every combination
+// of `Representation` and `Reduction` that PHAT supports.
+template <typename Reduction, typename Representation>
+void define_compute_persistence(py::module &mod,
+ const std::string &representation_suffix,
+ const std::string &reduction_suffix) {
+
+ auto suffix = representation_suffix + std::string("_") + reduction_suffix;
+
+ //We don't annotate these with doc strings or py::args because
+ //they are only used internally by code in phat.py
+ mod.def((std::string("compute_persistence_pairs_") + suffix).c_str(),
+ [](phat::boundary_matrix<Representation> &matrix){
+ phat::persistence_pairs pairs;
+ phat::compute_persistence_pairs<Reduction>(pairs, matrix);
+ return pairs;
+ });
+ mod.def((std::string("compute_persistence_pairs_dualized_") + suffix).c_str(),
+ [](phat::boundary_matrix<Representation> &matrix){
+ phat::persistence_pairs pairs;
+ phat::compute_persistence_pairs_dualized<Reduction>(pairs, matrix);
+ return pairs;
+ });
+}
+
+// Define a function to convert a `boundary_matrix` with one internal representation to a
+// `boundary_matrix` with a different internal representation. Like with define_compute_persistence,
+// the user will never see this function, but it is used internally by phat.py.
+template <typename SelfRep, typename OtherRep>
+void define_converter(py::module &mod, const std::string &self_suffix, const std::string &other_suffix) {
+ //We don't annotate these with doc strings or py::args because
+ //they are only used internally by code in phat.py
+ mod.def((std::string("convert_") + other_suffix + "_to_" + self_suffix).c_str(),
+ [](phat::boundary_matrix<OtherRep> &other) {
+ return phat::boundary_matrix<SelfRep>(other);
+ });
+}
+
+// Creates a Python class for a `boundary_matrix<T>`. Boundary matrices are one of two important types
+// used by PHAT.
+template<class T>
+void wrap_boundary_matrix(py::module &mod, const std::string &representation_suffix) {
+
+ using mat = phat::boundary_matrix<T>;
+
+ py::class_<mat>(mod, (std::string("boundary_matrix_") + representation_suffix).c_str())
+ //Default no-args constructor
+ .def(py::init())
+ //#### Loading and extracting data from Python lists
+ //Note we can use references to member functions (even templated ones) directly in many cases.
+ .def("load_vector_vector",
+ &mat::template load_vector_vector<phat::index, phat::dimension>,
+ "Load this instance with the given columns and dimensions",
+ py::arg("columns"), py::arg("dimensions"))
+ .def("get_vector_vector", [](mat &m) {
+ std::vector< std::vector< int > > vector_vector_matrix;
+ std::vector< int > vector_dims;
+ m.save_vector_vector( vector_vector_matrix, vector_dims );
+ return std::tuple<std::vector<std::vector<int>>, std::vector<int>>(vector_vector_matrix, vector_dims);
+ },
+ "Extract the data in the boundary matrix into a list of columns, and a list of dimensions that correspond to the columns")
+ //#### Loading and saving files
+ .def("load_binary", &mat::load_binary,
+ "Load this instance with data from a binary file")
+ .def("save_binary", &mat::save_binary,
+ "Save this instance to a binary file")
+ .def("load_ascii", &mat::load_ascii,
+ "Load this instance with data from a text file")
+ .def("save_ascii", &mat::save_ascii,
+ "Save this instance to a text file")
+ //#### Getting and setting dimensions
+ //Note that boundary_matrix dimensions are not normal matrix dimensions,
+ //They refer to the dimension of the simplex stored in the given column.
+ .def("get_dim", &mat::get_dim,
+ "Get the dimension list for this boundary matrix. "
+ "The dimension list shows the dimension for each column in the matrix")
+ .def("set_dim", &mat::set_dim, "Set the dimension for a single column",
+ py::arg("index"), py::arg("dimension"))
+ //The `set_dims` method is an example of making the data structure easier to use
+ //from Python. This is a method that doesn't exist in the C++ class, but we add it
+ //using a C++ lambda. This ability to enhance the binding using lambdas
+ //is an *extremely* handy tool.
+ .def("set_dims", [](mat &m, std::vector<phat::index> dims) {
+ m.set_num_cols(dims.size());
+ for(size_t i = 0; i < dims.size(); i++) {
+ m.set_dim(i, dims[i]);
+ }
+ },
+ "Set the dimension list for this boundary matrix",
+ py::arg("dimensions"))
+
+ //#### \__eq__
+ //The `boundary_matrix<T>`'s `operator==` is templated, which could make a Python wrapper
+ //very tricky indeed. Luckily, when we define multiple
+ //methods with the same name but different C++ types, pybind11 will create a Python method
+ //that chooses between them based on type tags that it manages. This is *also* extremely handy.
+ .def("__eq__", &mat::template operator==<phat::bit_tree_pivot_column>)
+ .def("__eq__", &mat::template operator==<phat::sparse_pivot_column>)
+ .def("__eq__", &mat::template operator==<phat::heap_pivot_column>)
+ .def("__eq__", &mat::template operator==<phat::full_pivot_column>)
+ .def("__eq__", &mat::template operator==<phat::vector_vector>)
+ .def("__eq__", &mat::template operator==<phat::vector_heap>)
+ .def("__eq__", &mat::template operator==<phat::vector_set>)
+ .def("__eq__", &mat::template operator==<phat::vector_list>)
+
+ //#### Data access
+
+ // In `get_col`, since Python is garbage collected, the C++ idiom of passing in a collection
+ // to load doesn't make much sense. We can simply allocate an STL vector and
+ // return it. The pybind11 framework will take ownership and hook it into the
+ // Python reference counting system.
+ .def("get_col", [](mat &m, phat::index col_index) {
+ std::vector<phat::index> col;
+ m.get_col(col_index, col);
+ return col;
+ },
+ "Extract a single column as a list",
+ py::arg("index"))
+ .def("set_col", &mat::set_col,
+ "Set the values for a given column",
+ py::arg("index"), py::arg("column"))
+ .def("get_num_cols", &mat::get_num_cols)
+ .def("is_empty", &mat::is_empty)
+ .def("get_num_entries", &mat::get_num_entries);
+
+ //#### Compute persistence
+ // Define compute_persistence(_dualized) for all possible reductions.
+ define_compute_persistence<phat::standard_reduction, T>(mod, representation_suffix, std::string("sr"));
+ define_compute_persistence<phat::chunk_reduction, T>(mod, representation_suffix, std::string("cr"));
+ define_compute_persistence<phat::row_reduction, T>(mod, representation_suffix, std::string("rr"));
+ define_compute_persistence<phat::twist_reduction, T>(mod, representation_suffix, std::string("tr"));
+ define_compute_persistence<phat::spectral_sequence_reduction, T>(mod, representation_suffix, std::string("ssr"));
+ //#### Converters
+ //Define functions to convert from this kind of `boundary_matrix` to any of the other types
+ define_converter<T, phat::bit_tree_pivot_column>(mod, representation_suffix, std::string("btpc"));
+ define_converter<T, phat::sparse_pivot_column>(mod, representation_suffix, std::string("spc"));
+ define_converter<T, phat::heap_pivot_column>(mod, representation_suffix, std::string("hpc"));
+ define_converter<T, phat::full_pivot_column>(mod, representation_suffix, std::string("fpc"));
+ define_converter<T, phat::vector_vector>(mod, representation_suffix, std::string("vv"));
+ define_converter<T, phat::vector_heap>(mod, representation_suffix, std::string("vh"));
+ define_converter<T, phat::vector_set>(mod, representation_suffix, std:: string("vs"));
+ define_converter<T, phat::vector_list>(mod, representation_suffix, std::string("vl"));
+}
+//fix_index checks for out-of-bounds indexes, and converts negative indices to positive ones
+//e.g. pairs[-1] => pairs[len(pairs) - 1]
+phat::index fix_index(const phat::persistence_pairs &p, int index) {
+ //Note get_num_pairs returns type index, which is not unsigned, though it comes from
+ //std::vector.size, which is size_t.
+ phat::index pairs = p.get_num_pairs();
+ assert(pairs > 0);
+ if (index < 0) {
+ index = pairs + index;
+ }
+ if ((index < 0) || static_cast<size_t>(index) >= static_cast<size_t>(pairs)) {
+ //pybind11 helpfully converts C++ exceptions into Python ones
+ throw py::index_error();
+ }
+ return index;
+}
+
+//Here we define the wrapper for the persistence_pairs class. Unlike `boundary_matrix`, this
+//class is not templated, so is simpler to wrap.
+void wrap_persistence_pairs(py::module &m) {
+ py::class_<phat::persistence_pairs>(m, "persistence_pairs")
+ //No-args constructor
+ .def(py::init())
+
+ //This is a method that takes two ints
+ .def("append_pair",
+ &phat::persistence_pairs::append_pair,
+ "Appends a single (birth, death) pair",
+ py::arg("birth"), py::arg("death"))
+
+ //This is a method that takes two ints
+ .def("set_pair",
+ &phat::persistence_pairs::set_pair,
+ "Sets the (birth, death) pair at a given index",
+ py::arg("index"), py::arg("birth"), py::arg("death"))
+
+ //#### Python collection support
+ .def("__len__", &phat::persistence_pairs::get_num_pairs)
+ // Unlike set_pair, this takes a Python 2-tuple
+ .def("__setitem__",
+ [](phat::persistence_pairs &p, int index, std::pair<phat::index,phat::index> &pair) {
+ phat::index idx = fix_index(p, index);
+ p.set_pair(idx, pair.first, pair.second);
+ })
+ // \__len\__ and \__getitem\__ together serve to make this a Python iterable
+ // so you can do `for i in pairs: blah`. A nicer way is to support \__iter\__,
+ // which we leave for future work.
+ .def("__getitem__", [](const phat::persistence_pairs &p, int index) {
+ phat::index idx = fix_index(p, index);
+ return p.get_pair(idx);
+ })
+ .def("clear", &phat::persistence_pairs::clear, "Empties the collection")
+ .def("sort", &phat::persistence_pairs::sort, "Sort in place")
+ .def("__eq__", &phat::persistence_pairs::operator==)
+ //#### File operations
+ .def("load_ascii", &phat::persistence_pairs::load_ascii,
+ "Load the contents of a text file into this instance")
+ .def("save_ascii", &phat::persistence_pairs::save_ascii,
+ "Save this instance to a text file")
+ .def("save_binary", &phat::persistence_pairs::save_binary,
+ "Save the contents of this instance to a binary file")
+ .def("load_binary", &phat::persistence_pairs::load_binary,
+ "Load the contents of a binary file into this instance");
+}
+
+//## Define the module
+//This is where we actually define the `_phat` module. We'll also have a `phat` module that's written
+//in Python, which will use `_phat` as an implementation detail.
+PYBIND11_PLUGIN(_phat) {
+ //Create the module object. First arg is the name, second is the module docstring.
+ py::module m("_phat", "C++ wrapper for PHAT. Please use the phat module, not the _phat module");
+
+ //Wrap the `persistence_pairs` class
+ wrap_persistence_pairs(m);
+
+ //#### Generate all the different representations of `boundary_matrix`
+ wrap_boundary_matrix<phat::bit_tree_pivot_column>(m, "btpc");
+ wrap_boundary_matrix<phat::sparse_pivot_column>(m, "spc");
+ wrap_boundary_matrix<phat::heap_pivot_column>(m, "hpc");
+ wrap_boundary_matrix<phat::full_pivot_column>(m, "fpc");
+ wrap_boundary_matrix<phat::vector_vector>(m, "vv");
+ wrap_boundary_matrix<phat::vector_heap>(m, "vh");
+ wrap_boundary_matrix<phat::vector_set>(m, "vs");
+ wrap_boundary_matrix<phat::vector_list>(m, "vl");
+
+ //We're done!
+ return m.ptr();
+
+}
diff --git a/python/phat.py b/python/phat.py
new file mode 100644
index 0000000..7cc0ae2
--- /dev/null
+++ b/python/phat.py
@@ -0,0 +1,88 @@
+import _phat
+import enum
+
+from _phat import persistence_pairs
+
+__all__ = ['boundary_matrix',
+ 'persistence_pairs',
+ 'compute_persistence_pairs',
+ 'compute_persistence_pairs_dualized']
+
+
+"""Bindings for the Persistent Homology Algorithm Toolbox
+
+
+Please see https://bitbucket.org/phat-code/phat for more information.
+"""
+
+
+class representations(enum.Enum):
+ """Available representations for internal storage of columns in
+ a `boundary_matrix`
+ """
+ bit_tree_pivot_column = 1
+ sparse_pivot_column = 2
+ full_pivot_column = 3
+ vector_vector = 4
+ vector_heap = 5
+ vector_set = 6
+ vector_list = 7
+
+class reductions(enum.Enum):
+ "Available reduction algorithms"
+ twist_reduction = 1
+ chunk_reduction = 2
+ standard_reduction = 3
+ row_reduction = 4
+ spectral_sequence_reduction = 5
+
+def __short_name(name):
+ return "".join([n[0] for n in name.split("_")])
+
+def convert(source, to_representation):
+ """Internal - function to convert from one `boundary_matrix` implementation to another"""
+ class_name = source.__class__.__name__
+ source_rep_short_name = class_name[len('boundary_matrix_'):]
+ to_rep_short_name = __short_name(to_representation.name)
+ function = getattr(_phat, "convert_%s_to_%s" % (source_rep_short_name, to_rep_short_name))
+ return function(source)
+
+def boundary_matrix(representation = representations.bit_tree_pivot_column, source = None):
+ """Returns an instance of a `boundary_matrix` class.
+ The boundary matrix will use the specified implementation for storing its column data.
+ If the `source` parameter is specified, it will be assumed to be another boundary matrix,
+ whose data should be copied into the new matrix.
+ """
+ if source:
+ return convert(source, representation)
+ else:
+ class_name = representation.name
+ short_name = __short_name(class_name)
+ function = getattr(_phat, "boundary_matrix_" + short_name)
+ return function()
+
+def compute_persistence_pairs(matrix,
+ reduction = reductions.twist_reduction):
+ """Computes persistence pairs (birth, death) for the given boundary matrix."""
+ class_name = matrix.__class__.__name__
+ representation_short_name = class_name[len('boundary_matrix_'):]
+ algo_name = reduction.name
+ algo_short_name = __short_name(algo_name)
+ function = getattr(_phat, "compute_persistence_pairs_" + representation_short_name + "_" + algo_short_name)
+ return function(matrix)
+
+def compute_persistence_pairs_dualized(matrix,
+ reduction = reductions.twist_reduction):
+ """Computes persistence pairs (birth, death) from the dualized form of the given boundary matrix."""
+ class_name = matrix.__class__.__name__
+ representation_short_name = class_name[len('boundary_matrix_'):]
+ algo_name = reduction.name
+ algo_short_name = __short_name(algo_name)
+ function = getattr(_phat, "compute_persistence_pairs_dualized_" + representation_short_name + "_" + algo_short_name)
+ return function(matrix)
+
+
+
+
+
+
diff --git a/python/setup.cfg b/python/setup.cfg
new file mode 100644
index 0000000..458847f
--- /dev/null
+++ b/python/setup.cfg
@@ -0,0 +1,5 @@
+# [bdist_wheel]
+# This flag says that the code is written to work on both Python 2 and Python
+# 3. If at all possible, it is good practice to do this. If you cannot, you
+# will need to generate wheels for each Python version that you support.
+# universal=1 \ No newline at end of file
diff --git a/python/setup.py b/python/setup.py
new file mode 100644
index 0000000..5e71cc7
--- /dev/null
+++ b/python/setup.py
@@ -0,0 +1,57 @@
+from setuptools import setup, Extension, find_packages
+from setuptools.command.build_ext import build_ext
+import sys
+import os.path
+ext_modules = [
+ Extension(
+ '_phat',
+ ['_phat.cpp'],
+ include_dirs=['include', '../include'],
+ language='c++',
+ ),
+]
+
+here = os.path.abspath(os.path.dirname(__file__))
+
+# Get the long description from the README file
+with open(os.path.join(here, 'README.rst'), encoding='utf-8') as f:
+ long_description = f.read()
+
+class BuildExt(build_ext):
+ """A custom build extension for adding compiler-specific options."""
+ c_opts = {
+ 'msvc': ['/EHsc'],
+ 'unix': ['-std=c++11'],
+ }
+
+ if sys.platform == 'darwin':
+ c_opts['unix'] += ['-stdlib=libc++', '-mmacosx-version-min=10.7']
+
+ def build_extensions(self):
+ ct = self.compiler.compiler_type
+ opts = self.c_opts.get(ct, [])
+ for ext in self.extensions:
+ ext.extra_compile_args = opts
+ build_ext.build_extensions(self)
+
+setup(
+ name='phat',
+ version='0.0.1',
+ author='Bryn Keller',
+ author_email='bryn.keller@intel.com',
+ url='https://bitbucket.org/phat-code/phat',
+ description='Python bindings for PHAT',
+ license = 'LGPL',
+ keywords='algebraic-topology PHAT distributed topology persistent-homology',
+ long_description=long_description,
+ ext_modules=ext_modules,
+ install_requires=['pybind11'],
+ cmdclass={'build_ext': BuildExt},
+ py_modules = ['phat'],
+ # packages = find_packages(exclude = ['doc', 'test'])
+ )
+
+
+
+
+
diff --git a/python/src/self_test.py b/python/src/self_test.py
new file mode 100644
index 0000000..c8174fb
--- /dev/null
+++ b/python/src/self_test.py
@@ -0,0 +1,168 @@
+import sys
+import phat
+
+if __name__=='__main__':
+ test_data = (sys.argv[1:] and sys.argv[1]) or "../../examples/torus.bin"
+
+ print("Reading test data %s in binary format ..." % test_data)
+
+ boundary_matrix = phat.boundary_matrix()
+ # This is broken for some reason
+ if not boundary_matrix.load_binary(test_data):
+ # if not boundary_matrix.load_ascii(test_data):
+ print("Error: test data %s not found!" % test_data)
+ sys.exit(1)
+
+ error = False
+
+ def compute_chunked(mat):
+ return phat.compute_persistence_pairs(mat, phat.reductions.chunk_reduction)
+
+ print("Comparing representations using Chunk algorithm ...")
+ print("Running Chunk - Sparse ...")
+ sparse_boundary_matrix = phat.boundary_matrix(phat.representations.sparse_pivot_column, boundary_matrix)
+ sparse_pairs = compute_chunked(sparse_boundary_matrix)
+
+ print("Running Chunk - Heap ...")
+ heap_boundary_matrix = phat.boundary_matrix(phat.representations.vector_heap, boundary_matrix)
+ heap_pairs = compute_chunked(heap_boundary_matrix)
+
+ print("Running Chunk - Full ...")
+ full_boundary_matrix = phat.boundary_matrix(phat.representations.full_pivot_column, boundary_matrix)
+ full_pairs = compute_chunked(full_boundary_matrix)
+
+ print("Running Chunk - BitTree ...")
+ bit_tree_boundary_matrix = phat.boundary_matrix(phat.representations.bit_tree_pivot_column, boundary_matrix)
+ bit_tree_pairs = compute_chunked(bit_tree_boundary_matrix)
+
+ print("Running Chunk - Vec_vec ...")
+ vec_vec_boundary_matrix = phat.boundary_matrix(phat.representations.vector_vector, boundary_matrix)
+ vec_vec_pairs = compute_chunked(vec_vec_boundary_matrix)
+
+ print("Running Chunk - Vec_heap ...")
+ vec_heap_boundary_matrix = phat.boundary_matrix(phat.representations.vector_heap, boundary_matrix)
+ vec_heap_pairs = compute_chunked(vec_heap_boundary_matrix)
+
+ print("Running Chunk - Vec_set ...")
+ vec_set_boundary_matrix = phat.boundary_matrix(phat.representations.vector_set, boundary_matrix)
+ vec_set_pairs = compute_chunked(vec_set_boundary_matrix)
+
+ print("Running Chunk - Vec_list ...")
+ vec_list_boundary_matrix = phat.boundary_matrix(phat.representations.vector_list, boundary_matrix)
+ vec_list_pairs = compute_chunked(vec_list_boundary_matrix)
+
+ if sparse_pairs != heap_pairs:
+ print("Error: sparse and heap differ!", file=sys.stderr)
+ error = True
+ if heap_pairs != full_pairs:
+ print("Error: heap and full differ!", file=sys.stderr)
+ error = True
+ if full_pairs != vec_vec_pairs:
+ print("Error: full and vec_vec differ!", file=sys.stderr)
+ error = True
+ if vec_vec_pairs != vec_heap_pairs:
+ print("Error: vec_vec and vec_heap differ!", file=sys.stderr)
+ error = True
+ if vec_heap_pairs != vec_set_pairs:
+ print("Error: vec_heap and vec_set differ!", file=sys.stderr)
+ error = True
+ if vec_set_pairs != bit_tree_pairs:
+ print("Error: vec_set and bit_tree differ!", file=sys.stderr)
+ error = True
+ if bit_tree_pairs != vec_list_pairs:
+ print("Error: bit_tree and vec_list differ!", file=sys.stderr)
+ error = True
+ if vec_list_pairs != sparse_pairs:
+ print("Error: vec_list and sparse differ!", file=sys.stderr)
+ error = True
+ if error:
+ sys.exit(1)
+ else:
+ print("All results are identical (as they should be)")
+
+ print("Comparing algorithms using BitTree representation ...")
+ print("Running Twist - BitTree ...")
+
+ def bit_tree_mat():
+ return phat.boundary_matrix(phat.representations.bit_tree_pivot_column, boundary_matrix)
+
+ reps = phat.representations
+ reds = phat.reductions
+ pairs = phat.compute_persistence_pairs
+
+ twist_boundary_matrix = bit_tree_mat()
+ twist_pairs = pairs(twist_boundary_matrix, reds.twist_reduction)
+
+ print("Running Standard - BitTree ...")
+ std_boundary_matrix = bit_tree_mat()
+ std_pairs = pairs(std_boundary_matrix, reds.standard_reduction)
+
+ print("Running Chunk - BitTree ...")
+ chunk_boundary_matrix = bit_tree_mat()
+ chunk_pairs = pairs(chunk_boundary_matrix, reds.chunk_reduction)
+
+ print("Running Row - BitTree ...")
+ row_boundary_matrix = bit_tree_mat()
+ row_pairs = pairs(row_boundary_matrix, reds.row_reduction)
+
+ print("Running Spectral sequence - BitTree ...")
+ ss_boundary_matrix = bit_tree_mat()
+ ss_pairs = pairs(ss_boundary_matrix, reds.spectral_sequence_reduction)
+
+ if twist_pairs != std_pairs:
+ print("Error: twist and standard differ!", file=sys.stderr)
+ error = True
+ if std_pairs != chunk_pairs:
+ print("Error: standard and chunk differ!", file=sys.stderr)
+ error = True
+ if chunk_pairs != row_pairs:
+ print("Error: chunk and row differ!", file=sys.stderr)
+ error = True
+ if row_pairs != ss_pairs:
+ print("Error: row and spectral sequence differ!", file=sys.stderr)
+ error = True
+ if ss_pairs != twist_pairs:
+ print("Error: spectral sequence and twist differ!", file=sys.stderr)
+ error = True
+ if error:
+ sys.exit(1)
+ else:
+ print("All results are identical (as they should be)")
+
+ print("Comparing primal and dual approach using Chunk - Full ...")
+
+ primal_boundary_matrix = phat.boundary_matrix(reps.full_pivot_column, boundary_matrix)
+ primal_pairs = phat.compute_persistence_pairs(primal_boundary_matrix, reds.chunk_reduction)
+
+ dual_boundary_matrix = phat.boundary_matrix(reps.full_pivot_column, boundary_matrix)
+ dual_pairs = phat.compute_persistence_pairs_dualized(dual_boundary_matrix)
+
+ if primal_pairs != dual_pairs:
+ print("Error: primal and dual differ!", file=sys.stderr)
+ error = True
+
+
+ if error:
+ sys.exit(1)
+ else:
+ print("All results are identical (as they should be)")
+
+ print("Testing vector<vector> interface ...")
+
+ (vector_vector_matrix, vector_dims) = boundary_matrix.get_vector_vector()
+
+ vector_vector_boundary_matrix = phat.boundary_matrix(phat.representations.bit_tree_pivot_column)
+
+ vector_vector_boundary_matrix.load_vector_vector(vector_vector_matrix, vector_dims)
+
+ if vector_vector_boundary_matrix != boundary_matrix:
+ print("Error: [load|save]_vector_vector bug", file=sys.stderr)
+ error = True
+
+
+ if error:
+ sys.exit(1)
+ else:
+ print("Test passed!")
+
+
diff --git a/python/src/simple_example.py b/python/src/simple_example.py
new file mode 100644
index 0000000..82cf6be
--- /dev/null
+++ b/python/src/simple_example.py
@@ -0,0 +1,59 @@
+"""This file contains a simple example that demonstrates the usage of the library interface"""
+
+if __name__ == "__main__":
+
+ print("""
+ we will build an ordered boundary matrix of this simplicial complex consisting of a single triangle:
+
+ 3
+ |\\
+ | \\
+ | \\
+ | \\ 4
+ 5| \\
+ | \\
+ | 6 \\
+ | \\
+ |________\\
+ 0 2 1
+
+""")
+
+ import phat
+
+ # set the dimension of the cell that each column represents:
+ dimensions = [0, 0, 1, 0, 1, 1, 2]
+
+ # define a boundary matrix with the chosen internal representation
+ boundary_matrix = phat.boundary_matrix(representation = phat.representations.vector_vector)
+
+ # set the respective columns -- the columns entries have to be sorted
+ boundary_matrix.set_dims(dimensions)
+ boundary_matrix.set_col(0, [])
+ boundary_matrix.set_col(1, [])
+ boundary_matrix.set_col(2, [0,1])
+ boundary_matrix.set_col(3, [])
+ boundary_matrix.set_col(4, [1,3])
+ boundary_matrix.set_col(5, [0,3])
+ boundary_matrix.set_col(6, [2,4,5])
+
+ # print some information of the boundary matrix:
+ print()
+ print("The boundary matrix has %d columns:" % boundary_matrix.get_num_cols())
+ for col_idx in range(boundary_matrix.get_num_cols()):
+ s = "Column %d represents a cell of dimension %d." % (col_idx, boundary_matrix.get_dim(col_idx))
+ if (not boundary_matrix.is_empty(col_idx)):
+ s = s + " Its boundary consists of the cells " + " ".join([str(c) for c in boundary_matrix.get_col(col_idx)])
+ print(s)
+ print("Overall, the boundary matrix has %d entries." % boundary_matrix.get_num_entries())
+
+ pairs = phat.compute_persistence_pairs(boundary_matrix)
+
+ pairs.sort()
+
+ print()
+
+ print("There are %d persistence pairs: " % len(pairs))
+ for pair in pairs:
+ print("Birth: %d, Death: %d" % pair)
+
diff --git a/src/simple_example.cpp b/src/simple_example.cpp
index aff3ff9..ba8aafe 100644
--- a/src/simple_example.cpp
+++ b/src/simple_example.cpp
@@ -96,7 +96,7 @@ int main( int argc, char** argv )
std::cout << std::endl;
std::cout << "The boundary matrix has " << boundary_matrix.get_num_cols() << " columns: " << std::endl;
for( phat::index col_idx = 0; col_idx < boundary_matrix.get_num_cols(); col_idx++ ) {
- std::cout << "Colum " << col_idx << " represents a cell of dimension " << (int)boundary_matrix.get_dim( col_idx ) << ". ";
+ std::cout << "Column " << col_idx << " represents a cell of dimension " << (int)boundary_matrix.get_dim( col_idx ) << ". ";
if( !boundary_matrix.is_empty( col_idx ) ) {
std::vector< phat::index > temp_col;
boundary_matrix.get_col( col_idx, temp_col );