From d7677b1f3f5ceee9cfdb5275ef00cbf79a714122 Mon Sep 17 00:00:00 2001 From: Bryn Keller Date: Mon, 7 Mar 2016 13:51:48 -0800 Subject: Python interface to PHAT --- .gitignore | 13 ++ README.md | 4 +- include/phat/boundary_matrix.h | 8 +- python/.gitignore | 7 ++ python/MANIFEST.in | 1 + python/README.rst | 110 ++++++++++++++++ python/_phat.cpp | 279 +++++++++++++++++++++++++++++++++++++++++ python/phat.py | 88 +++++++++++++ python/setup.cfg | 5 + python/setup.py | 57 +++++++++ python/src/self_test.py | 168 +++++++++++++++++++++++++ python/src/simple_example.py | 59 +++++++++ src/simple_example.cpp | 2 +- 13 files changed, 796 insertions(+), 5 deletions(-) create mode 100644 .gitignore create mode 100644 python/.gitignore create mode 100644 python/MANIFEST.in create mode 100644 python/README.rst create mode 100644 python/_phat.cpp create mode 100644 python/phat.py create mode 100644 python/setup.cfg create mode 100644 python/setup.py create mode 100644 python/src/self_test.py create mode 100644 python/src/simple_example.py 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 + +//Automatic conversions of stl containers to Python ones +#include + +//Additional support for operators and numpy +#include +#include + +//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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +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 +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 &matrix){ + phat::persistence_pairs pairs; + phat::compute_persistence_pairs(pairs, matrix); + return pairs; + }); + mod.def((std::string("compute_persistence_pairs_dualized_") + suffix).c_str(), + [](phat::boundary_matrix &matrix){ + phat::persistence_pairs pairs; + phat::compute_persistence_pairs_dualized(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 +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 &other) { + return phat::boundary_matrix(other); + }); +} + +// Creates a Python class for a `boundary_matrix`. Boundary matrices are one of two important types +// used by PHAT. +template +void wrap_boundary_matrix(py::module &mod, const std::string &representation_suffix) { + + using mat = phat::boundary_matrix; + + py::class_(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, + "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>(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 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`'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==) + .def("__eq__", &mat::template operator==) + .def("__eq__", &mat::template operator==) + .def("__eq__", &mat::template operator==) + .def("__eq__", &mat::template operator==) + .def("__eq__", &mat::template operator==) + .def("__eq__", &mat::template operator==) + .def("__eq__", &mat::template operator==) + + //#### 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 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(mod, representation_suffix, std::string("sr")); + define_compute_persistence(mod, representation_suffix, std::string("cr")); + define_compute_persistence(mod, representation_suffix, std::string("rr")); + define_compute_persistence(mod, representation_suffix, std::string("tr")); + define_compute_persistence(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(mod, representation_suffix, std::string("btpc")); + define_converter(mod, representation_suffix, std::string("spc")); + define_converter(mod, representation_suffix, std::string("hpc")); + define_converter(mod, representation_suffix, std::string("fpc")); + define_converter(mod, representation_suffix, std::string("vv")); + define_converter(mod, representation_suffix, std::string("vh")); + define_converter(mod, representation_suffix, std:: string("vs")); + define_converter(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(index) >= static_cast(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_(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 &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(m, "btpc"); + wrap_boundary_matrix(m, "spc"); + wrap_boundary_matrix(m, "hpc"); + wrap_boundary_matrix(m, "fpc"); + wrap_boundary_matrix(m, "vv"); + wrap_boundary_matrix(m, "vh"); + wrap_boundary_matrix(m, "vs"); + wrap_boundary_matrix(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 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 ); -- cgit v1.2.3