" << std::endl
#ifdef USE_COEFFICIENTS
<< " --modulus compute homology with coefficients in the prime field Z/
Z"
#endif
<< std::endl;
exit(exit_code);
}
int main(int argc, char** argv) {
const char* filename = nullptr;
file_format format = DISTANCE_MATRIX;
index_t dim_max = 1;
value_t threshold = std::numeric_limits::max();
#ifdef USE_COEFFICIENTS
coefficient_t modulus = 2;
#else
const coefficient_t modulus = 2;
#endif
for (index_t i = 1; i < argc; ++i) {
const std::string arg(argv[i]);
if (arg == "--help") {
print_usage_and_exit(0);
} else if (arg == "--dim") {
std::string parameter = std::string(argv[++i]);
size_t next_pos;
dim_max = std::stol(parameter, &next_pos);
if (next_pos != parameter.size()) print_usage_and_exit(-1);
} else if (arg == "--threshold") {
std::string parameter = std::string(argv[++i]);
size_t next_pos;
threshold = std::stof(parameter, &next_pos);
if (next_pos != parameter.size()) print_usage_and_exit(-1);
} else if (arg == "--format") {
std::string parameter = std::string(argv[++i]);
if (parameter == "lower-distance")
format = LOWER_DISTANCE_MATRIX;
else if (parameter == "upper-distance")
format = UPPER_DISTANCE_MATRIX;
else if (parameter == "distance")
format = DISTANCE_MATRIX;
else if (parameter == "point-cloud")
format = POINT_CLOUD;
else if (parameter == "dipha")
format = DIPHA;
else
print_usage_and_exit(-1);
#ifdef USE_COEFFICIENTS
} else if (arg == "--modulus") {
std::string parameter = std::string(argv[++i]);
size_t next_pos;
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];
}
}
std::ifstream file_stream(filename);
if (filename && file_stream.fail()) {
std::cerr << "couldn't open file " << filename << std::endl;
exit(-1);
}
compressed_lower_distance_matrix dist = read_file(filename ? file_stream : std::cin, format);
index_t n = dist.size();
std::cout << "distance matrix with " << n << " points" << std::endl;
auto value_range = std::minmax_element(dist.distances.begin(), dist.distances.end());
std::cout << "value range: [" << *value_range.first << "," << *value_range.second << "]" << std::endl;
dim_max = std::min(dim_max, n - 2);
binomial_coeff_table binomial_coeff(n, dim_max + 2);
std::vector multiplicative_inverse(multiplicative_inverse_vector(modulus));
std::vector columns_to_reduce;
{
union_find dset(n);
std::vector edges;
rips_filtration_comparator comp(dist, 1, binomial_coeff);
for (index_t index = binomial_coeff(n, 2); index-- > 0;) {
value_t diameter = comp.diameter(index);
if (diameter <= threshold) edges.push_back(diameter_index_t(diameter, index));
}
std::sort(edges.rbegin(), edges.rend(), greater_diameter_or_smaller_index());
#ifdef PRINT_PERSISTENCE_PAIRS
std::cout << "persistence intervals in dim 0:" << std::endl;
#endif
std::vector vertices_of_edge(2);
for (auto e : edges) {
vertices_of_edge.clear();
get_simplex_vertices(get_index(e), 1, n, binomial_coeff, 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
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());
#ifdef PRINT_PERSISTENCE_PAIRS
for (index_t i = 0; i < n; ++i)
if (dset.find(i) == i) std::cout << " [0, )" << std::endl << std::flush;
#endif
}
for (index_t dim = 1; dim <= dim_max; ++dim) {
rips_filtration_comparator comp(dist, dim + 1, binomial_coeff);
rips_filtration_comparator comp_prev(dist, dim, binomial_coeff);
hash_map pivot_column_index;
pivot_column_index.reserve(columns_to_reduce.size());
compute_pairs(columns_to_reduce, pivot_column_index, dist, comp, comp_prev, dim, n, threshold, modulus,
multiplicative_inverse, binomial_coeff);
if (dim < dim_max) {
assemble_columns_to_reduce(columns_to_reduce, pivot_column_index, comp, dim, n, threshold, binomial_coeff);
}
}
}