diff options
| -rw-r--r-- | Makefile | 24 | ||||
| -rw-r--r-- | ripser.cpp | 118 | ||||
| -rw-r--r-- | ripser.xcodeproj/project.pbxproj | 1 |
3 files changed, 100 insertions, 43 deletions
@@ -1,24 +1,8 @@ -build: ripser +build: ripser-representatives -all: ripser ripser-coeff ripser-reduction ripser-coeff-reduction ripser-debug - - -ripser: ripser.cpp - c++ -std=c++11 ripser.cpp -o ripser -Ofast -D NDEBUG - -ripser-coeff: ripser.cpp - c++ -std=c++11 ripser.cpp -o ripser-coeff -Ofast -D NDEBUG -D USE_COEFFICIENTS - -ripser-reduction: ripser.cpp - c++ -std=c++11 ripser.cpp -o ripser-reduction -Ofast -D NDEBUG -D ASSEMBLE_REDUCTION_MATRIX - -ripser-coeff-reduction: ripser.cpp - c++ -std=c++11 ripser.cpp -o ripser-coeff-reduction -Ofast -D NDEBUG -D USE_COEFFICIENTS -D ASSEMBLE_REDUCTION_MATRIX - -ripser-debug: ripser.cpp - c++ -std=c++11 ripser.cpp -o ripser-debug -g - +ripser-representatives: ripser.cpp + c++ -std=c++11 ripser.cpp -o ripser-representatives -Ofast -D NDEBUG -D _USE_COEFFICIENTS -D ASSEMBLE_REDUCTION_MATRIX -D USE_GOOGLE_HASHMAP clean: - rm -f ripser ripser-coeff ripser-reduction ripser-coeff-reduction ripser-debug + rm -f ripser-representatives @@ -36,6 +36,8 @@ with this program. If not, see <http://www.gnu.org/licenses/>. #include <sstream> #include <unordered_map> +#include "prettyprint.hpp" + #ifdef USE_GOOGLE_HASHMAP #include <sparsehash/sparse_hash_map> template <class Key, class T> class hash_map : public google::sparse_hash_map<Key, T> { @@ -78,6 +80,10 @@ bool is_prime(const coefficient_t n) { return true; } +coefficient_t normalize(const coefficient_t n, const coefficient_t modulus) { + return n > modulus / 2 ? n - modulus : n; +} + std::vector<coefficient_t> multiplicative_inverse_vector(const coefficient_t m) { std::vector<coefficient_t> inverse(m); inverse[1] = 1; @@ -432,6 +438,12 @@ public: return out; } + std::vector<index_t> vertices_of_simplex(const index_t simplex_index, const index_t dim) { + std::vector<index_t> vertices(dim + 1); + get_simplex_vertices(simplex_index, dim, n, vertices.begin()); + return vertices; + } + value_t compute_diameter(const index_t index, index_t dim) const { value_t diam = -std::numeric_limits<value_t>::infinity(); @@ -461,31 +473,56 @@ public: greater_diameter_or_smaller_index<diameter_index_t>()); #ifdef PRINT_PERSISTENCE_PAIRS - std::cout << "persistence intervals in dim 0:" << std::endl; + std::cout << "persistence intervals in dim 0:" << std::endl; #endif - - std::vector<index_t> vertices_of_edge(2); - for (auto e : edges) { - vertices_of_edge.clear(); - get_simplex_vertices(get_index(e), 1, n, std::back_inserter(vertices_of_edge)); - index_t u = dset.find(vertices_of_edge[0]), v = dset.find(vertices_of_edge[1]); - - if (u != v) { + + std::vector<index_t> vertices_of_edge(2); + for (auto e : edges) { + vertices_of_edge.clear(); + get_simplex_vertices(get_index(e), 1, n, std::back_inserter(vertices_of_edge)); + index_t u = dset.find(vertices_of_edge[0]), v = dset.find(vertices_of_edge[1]); + + if (u != v) { #ifdef PRINT_PERSISTENCE_PAIRS - if (get_diameter(e) != 0) - std::cout << " [0," << get_diameter(e) << ")" << std::endl; + if (get_diameter(e) != 0) + std::cout << " [0," << get_diameter(e) << ")" << std::endl; #endif - dset.link(u, v); - } else - columns_to_reduce.push_back(e); - } - std::reverse(columns_to_reduce.begin(), columns_to_reduce.end()); - + std::cout << "{"; + bool nonempty = false; + for (index_t w = 0; w < n; ++w) + if (dset.find(w) == u) { + if (nonempty) std::cout << ", "; nonempty = true; + std::cout << "[" << w << "]"; +#ifdef USE_COEFFICIENTS + std::cout << ":" << 1; +#endif + } + std::cout << "}" << std::endl; + + + dset.link(u, v); + } else + columns_to_reduce.push_back(e); + } + std::reverse(columns_to_reduce.begin(), columns_to_reduce.end()); + #ifdef PRINT_PERSISTENCE_PAIRS - for (index_t i = 0; i < n; ++i) - if (dset.find(i) == i) std::cout << " [0, )" << std::endl << std::flush; + for (index_t i = 0; i < n; ++i) + if (dset.find(i) == i) { + std::cout << " [0, )" << std::endl << std::flush; + std::cout << "{"; + for (index_t w = 0; w < n; ++w) { + if (w > 0) std::cout << ", "; + std::cout << "[" << w << "]"; +#ifdef USE_COEFFICIENTS + std::cout << ":" << 1; #endif - } + } + std::cout << "}" << std::endl; +#endif + } + } + template <typename Column, typename Iterator> diameter_entry_t add_coboundary_and_get_pivot(Iterator column_begin, Iterator column_end, @@ -524,6 +561,34 @@ public: return get_pivot(working_coboundary, modulus); } + + template<typename Chain> + void print_chain(Chain& cycle, index_t dim) { + diameter_entry_t e; + + std::cout << "{"; + while (get_index(e = get_pivot(cycle, modulus)) != -1) { + vertices.resize(dim + 1); + get_simplex_vertices(get_index(e), dim, n, + vertices.rbegin()); + + std::cout << "["; + auto it = vertices.begin(); + if (it != vertices.end()) { + std::cout << *it++; + while (it != vertices.end()) std::cout << "," << *it++; + } + std::cout << "]"; +#ifdef USE_COEFFICIENTS + std::cout << ":" << normalize(get_coefficient(e), modulus); +#endif + cycle.pop(); + if (get_index(e = get_pivot(cycle, modulus)) != -1) + std::cout << ", "; + } + std::cout << "}" << std::endl; + + } void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce, hash_map<index_t, index_t>& pivot_column_index, index_t dim) { @@ -627,8 +692,9 @@ public: #ifdef INDICATE_PROGRESS std::cout << "\033[K"; #endif - std::cout << " [" << diameter << "," << death << ")" << std::endl - << std::flush; +// std::cout << " [" << diameter << "," << death << ")" << std::endl << std::flush; + std::cout << " [" << diameter << "," << death << "): "; + print_chain(working_reduction_column, dim); } #endif pivot_column_index.insert( @@ -668,7 +734,13 @@ public: } } else { #ifdef PRINT_PERSISTENCE_PAIRS - std::cout << " [" << diameter << ", )" << std::endl << std::flush; +#ifdef INDICATE_PROGRESS + std::cout << "\033[K"; +#endif +// std::cout << " [" << diameter << ", )" << std::endl << std::flush; + std::cout << " [" << diameter << ", ): "; + print_chain(working_reduction_column, dim); + #endif break; } diff --git a/ripser.xcodeproj/project.pbxproj b/ripser.xcodeproj/project.pbxproj index b5f6881..c2667d1 100644 --- a/ripser.xcodeproj/project.pbxproj +++ b/ripser.xcodeproj/project.pbxproj @@ -229,6 +229,7 @@ _FILE_FORMAT_LOWER_DISTANCE_MATRIX, _FILE_FORMAT_POINT_CLOUD, _USE_COEFFICIENTS, + ASSEMBLE_REDUCTION_MATRIX, INDICATE_PROGRESS, ); PRODUCT_NAME = "$(TARGET_NAME)"; |