diff options
author | Gard Spreemann <gard.spreemann@epfl.ch> | 2019-01-28 17:06:07 +0100 |
---|---|---|
committer | Gard Spreemann <gard.spreemann@epfl.ch> | 2019-01-28 17:06:07 +0100 |
commit | 4b0c30171a6fddc8a6dbee5290d7f536cdf978a8 (patch) | |
tree | 31f3cc00f029aea72d6e21d3d9587f5ead835af6 | |
parent | 2ec7e4d5d3dcd3d3970abe48129dc8db1737441d (diff) |
Basic output done, some basic sanity checks done.gspr/representative-cocycles-simpler-output
-rw-r--r-- | ripser.cpp | 162 |
1 files changed, 113 insertions, 49 deletions
@@ -48,10 +48,35 @@ public: template <class Key, class T> class hash_map : public std::unordered_map<Key, T> {}; #endif +static const int64_t DIPHA_MAGIC = 8067171840; +static const int64_t DIPHA_TYPE_DISTANCE_MATRIX = 7; +static const int64_t DIPHA_TYPE_PERSISTENCE_DIAGRAM = 2; + typedef float value_t; typedef int64_t index_t; typedef int16_t coefficient_t; +template <typename T> inline void reverse_endianness(T & x) { + uint8_t * p = reinterpret_cast<uint8_t *>(&x); + std::reverse(p, p + sizeof(T)); +} + +template <typename T> T read(std::istream& s) { + T result; + s.read(reinterpret_cast<char*>(&result), sizeof(T)); + #ifdef BIGENDIAN + reverse_endianness(result); + #endif + return result; +} + +template <typename T> void write(std::ostream& s, T value) { + #ifdef BIGENDIAN + reverse_endianness(value); + #endif + s.write(reinterpret_cast<char*>(&value), sizeof(T)); +} + class binomial_coeff_table { std::vector<std::vector<index_t>> B; @@ -458,7 +483,7 @@ public: hash_map<index_t, index_t>& pivot_column_index, index_t dim); void compute_dim_0_pairs(std::vector<diameter_index_t>& edges, - std::vector<diameter_index_t>& columns_to_reduce) { + std::vector<diameter_index_t>& columns_to_reduce, std::ofstream & outstream_pds, std::ofstream & outstream_reps, int64_t & pair_count) { union_find dset(n); edges = get_edges(); @@ -478,20 +503,23 @@ public: 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) { + write<int64_t>(outstream_pds, 0); + write<double>(outstream_pds, 0); + write<double>(outstream_pds, get_diameter(e)); + ++pair_count; + } #endif - 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 << "]"; + if (nonempty) outstream_reps << " "; nonempty = true; + outstream_reps << "[" << w << "]"; #ifdef USE_COEFFICIENTS - std::cout << ":" << 1; + outstream_reps << ":" << 1; #endif } - std::cout << "}" << std::endl; + outstream_reps << std::endl; dset.link(u, v); @@ -503,16 +531,18 @@ public: #ifdef PRINT_PERSISTENCE_PAIRS for (index_t i = 0; i < n; ++i) if (dset.find(i) == i) { - std::cout << " [0, )" << std::endl << std::flush; - std::cout << "{"; + write<int64_t>(outstream_pds, -0-1); // Indicate essential class. + write<double>(outstream_pds, 0); + write<double>(outstream_pds, NAN); // Undefined/unused for essential classes. + ++pair_count; for (index_t w = 0; w < n; ++w) { - if (w > 0) std::cout << ", "; - std::cout << "[" << w << "]"; + if (w > 0) outstream_reps << " "; + outstream_reps << "[" << w << "]"; #ifdef USE_COEFFICIENTS - std::cout << ":" << 1; + outstream_reps << ":" << 1; #endif } - std::cout << "}" << std::endl; + outstream_reps << std::endl; #endif } } @@ -529,35 +559,34 @@ public: bool& might_be_apparent_pair); template<typename Chain> - void print_chain(Chain& cycle, index_t dim) { + void write_chain(std::ofstream & stream, 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 << "["; + stream << "["; auto it = vertices.begin(); if (it != vertices.end()) { - std::cout << *it++; - while (it != vertices.end()) std::cout << "," << *it++; + stream << *it++; + while (it != vertices.end()) stream << "," << *it++; } - std::cout << "]"; + stream << "]"; #ifdef USE_COEFFICIENTS - std::cout << ":" << normalize(get_coefficient(e), modulus); + stream << ":" << normalize(get_coefficient(e), modulus); #endif cycle.pop(); if (get_index(e = get_pivot(cycle, modulus)) != -1) - std::cout << ", "; + stream << " "; } - std::cout << "}" << std::endl; + stream << 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) { + hash_map<index_t, index_t>& pivot_column_index, index_t dim, std::ofstream & outstream_pds, std::ofstream & outstream_reps, int64_t & pair_count) { #ifdef PRINT_PERSISTENCE_PAIRS std::cout << "persistence intervals in dim " << dim << ":" << std::endl; @@ -658,9 +687,11 @@ public: #ifdef INDICATE_PROGRESS std::cout << "\033[K"; #endif -// std::cout << " [" << diameter << "," << death << ")" << std::endl << std::flush; - std::cout << " [" << diameter << "," << death << "): "; - print_chain(working_reduction_column, dim); + write<int64_t>(outstream_pds, dim); + write<double>(outstream_pds, diameter); + write<double>(outstream_pds, death); + ++pair_count; + write_chain(outstream_reps, working_reduction_column, dim); } #endif pivot_column_index.insert( @@ -703,9 +734,11 @@ public: #ifdef INDICATE_PROGRESS std::cout << "\033[K"; #endif -// std::cout << " [" << diameter << ", )" << std::endl << std::flush; - std::cout << " [" << diameter << ", ): "; - print_chain(working_reduction_column, dim); + write<int64_t>(outstream_pds, -dim - 1); // Essential class in dimension dim + write<double>(outstream_pds, diameter); + write<double>(outstream_pds, NAN); // Undefined/unused for essential classes. + ++pair_count; + write_chain(outstream_reps, working_reduction_column, dim); #endif break; @@ -720,17 +753,17 @@ public: std::vector<diameter_index_t> get_edges(); - void compute_barcodes() { + void compute_barcodes(std::ofstream & outstream_pds, std::ofstream & outstream_reps, int64_t & pair_count) { std::vector<diameter_index_t> simplices, columns_to_reduce; - compute_dim_0_pairs(simplices, columns_to_reduce); + compute_dim_0_pairs(simplices, columns_to_reduce, outstream_pds, outstream_reps, pair_count); for (index_t dim = 1; dim <= dim_max; ++dim) { hash_map<index_t, index_t> pivot_column_index; pivot_column_index.reserve(columns_to_reduce.size()); - compute_pairs(columns_to_reduce, pivot_column_index, dim); + compute_pairs(columns_to_reduce, pivot_column_index, dim, outstream_pds, outstream_reps, pair_count); if (dim < dim_max) { assemble_columns_to_reduce(simplices, columns_to_reduce, pivot_column_index, @@ -1000,12 +1033,6 @@ enum file_format { RIPSER }; -template <typename T> T read(std::istream& s) { - T result; - s.read(reinterpret_cast<char*>(&result), sizeof(T)); - return result; // on little endian: boost::endian::little_to_native(result); -} - compressed_lower_distance_matrix read_point_cloud(std::istream& input_stream) { std::vector<std::vector<value_t>> points; @@ -1154,7 +1181,8 @@ void print_usage_and_exit(int exit_code) { int main(int argc, char** argv) { - const char* filename = nullptr; + const char* ifilename = nullptr; + const char* oprefix = nullptr; file_format format = DISTANCE_MATRIX; @@ -1211,19 +1239,48 @@ int main(int argc, char** argv) { modulus = std::stol(parameter, &next_pos); if (next_pos != parameter.size() || !is_prime(modulus)) print_usage_and_exit(-1); #endif - } else { - if (filename) { print_usage_and_exit(-1); } - filename = argv[i]; + } else if (i == argc - 2) { + ifilename = argv[i]; + } else if (i == argc - 1) { + oprefix = argv[i]; + } else { + print_usage_and_exit(-1); } } - std::ifstream file_stream(filename); - if (filename && file_stream.fail()) { - std::cerr << "couldn't open file " << filename << std::endl; + if (ifilename == nullptr || oprefix == nullptr) + { + std::cerr << "Need input and output file names." << std::endl; + exit(-1); + } + + std::ifstream file_stream(ifilename); + if (ifilename && file_stream.fail()) { + std::cerr << "couldn't open file " << ifilename << std::endl; exit(-1); } - compressed_lower_distance_matrix dist = read_file(filename ? file_stream : std::cin, format); + std::ofstream outstream_pds(std::string(oprefix) + std::string(".pd.dipha"), std::ios::binary); + if (!outstream_pds.is_open()) + { + std::cerr << "Failed to open PD output file." << std::endl; + exit(-1); + } + + std::ofstream outstream_reps(std::string(oprefix) + std::string(".reps.txt")); + if (!outstream_reps.is_open()) + { + std::cerr << "Failed to open representative cocycles output file." << std::endl; + exit(-1); + } + + int64_t pair_count = 0; + + write<int64_t>(outstream_pds, DIPHA_MAGIC); + write<int64_t>(outstream_pds, DIPHA_TYPE_PERSISTENCE_DIAGRAM); + write<int64_t>(outstream_pds, DIPHA_MAGIC); // This is to later be replaced by the number of pairs in the persistence diagram. + + compressed_lower_distance_matrix dist = read_file(ifilename ? file_stream : std::cin, format); value_t min = std::numeric_limits<value_t>::infinity(), max = -std::numeric_limits<value_t>::infinity(), max_finite = max; @@ -1251,13 +1308,20 @@ int main(int argc, char** argv) { std::cout << "distance matrix with " << dist.size() << " points" << std::endl; ripser<compressed_lower_distance_matrix>(std::move(dist), dim_max, threshold, ratio, modulus) - .compute_barcodes(); + .compute_barcodes(outstream_pds, outstream_reps, pair_count); } else { std::cout << "sparse distance matrix with " << dist.size() << " points and " << num_edges << " entries" << std::endl; ripser<sparse_distance_matrix>(sparse_distance_matrix(std::move(dist), threshold), dim_max, threshold, ratio, modulus) - .compute_barcodes(); + .compute_barcodes(outstream_pds, outstream_reps, pair_count); } + + + outstream_pds.seekp(2*sizeof(int64_t), std::ios_base::beg); // Rewind to the spot in the DIPHA file where we write the number of intervals. + write<int64_t>(outstream_pds, pair_count); + + outstream_pds.close(); + outstream_reps.close(); } |