Skip to content

Commit

Permalink
cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
@ubauer
ubauer committed Aug 21, 2016
1 parent cb53a8a commit a76f7e8
Show file tree
Hide file tree
Showing 2 changed files with 47 additions and 86 deletions.
2 changes: 1 addition & 1 deletion .clang-format
View file
@@ -1,7 +1,7 @@
BasedOnStyle: LLVM
IndentWidth: 4
TabWidth: 4
ColumnLimit: 100
ColumnLimit: 120
AccessModifierOffset: -4
AllowShortIfStatementsOnASingleLine: true
AllowShortLoopsOnASingleLine: true
Expand Down
131 changes: 46 additions & 85 deletions ripser.cpp
View file
Expand Up @@ -96,8 +96,7 @@ std::vector<coefficient_t> multiplicative_inverse_vector(const coefficient_t m)

template <typename OutputIterator>
OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t n,
const binomial_coeff_table& binomial_coeff,
OutputIterator out) {
const binomial_coeff_table& binomial_coeff, OutputIterator out) {
--n;

for (index_t k = dim + 1; k > 0; --k) {
Expand All @@ -124,8 +123,7 @@ OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t n,
return out;
}

std::vector<index_t> vertices_of_simplex(const index_t simplex_index, const index_t dim,
const index_t n,
std::vector<index_t> vertices_of_simplex(const index_t simplex_index, const index_t dim, const index_t n,
const binomial_coeff_table& binomial_coeff) {
std::vector<index_t> vertices;
get_simplex_vertices(simplex_index, dim, n, binomial_coeff, std::back_inserter(vertices));
Expand All @@ -136,15 +134,12 @@ std::vector<index_t> vertices_of_simplex(const index_t simplex_index, const inde
struct entry_t {
index_t index;
coefficient_t coefficient;
entry_t(index_t _index, coefficient_t _coefficient)
: index(_index), coefficient(_coefficient) {}
entry_t(index_t _index, coefficient_t _coefficient) : index(_index), coefficient(_coefficient) {}
entry_t(index_t _index) : index(_index), coefficient(1) {}
entry_t() : index(0), coefficient(1) {}
};

entry_t make_entry(index_t _index, coefficient_t _coefficient) {
return entry_t(_index, _coefficient);
}
entry_t make_entry(index_t _index, coefficient_t _coefficient) { return entry_t(_index, _coefficient); }
index_t get_index(entry_t e) { return e.index; }
index_t get_coefficient(entry_t e) { return e.coefficient; }
void set_coefficient(entry_t& e, const coefficient_t c) { e.coefficient = c; }
Expand Down Expand Up @@ -188,20 +183,14 @@ class diameter_entry_t : public std::pair<value_t, entry_t> {
const entry_t& get_entry(const diameter_entry_t& p) { return p.second; }
entry_t& get_entry(diameter_entry_t& p) { return p.second; }
const index_t get_index(const diameter_entry_t& p) { return get_index(get_entry(p)); }
const coefficient_t get_coefficient(const diameter_entry_t& p) {
return get_coefficient(get_entry(p));
}
const coefficient_t get_coefficient(const diameter_entry_t& p) { return get_coefficient(get_entry(p)); }
const value_t& get_diameter(const diameter_entry_t& p) { return p.first; }
void set_coefficient(diameter_entry_t& p, const coefficient_t c) {
set_coefficient(get_entry(p), c);
}
diameter_entry_t make_diameter_entry(value_t _diameter, index_t _index,
coefficient_t _coefficient) {
void set_coefficient(diameter_entry_t& p, const coefficient_t c) { set_coefficient(get_entry(p), c); }
diameter_entry_t make_diameter_entry(value_t _diameter, index_t _index, coefficient_t _coefficient) {
return std::make_pair(_diameter, make_entry(_index, _coefficient));
}
diameter_entry_t make_diameter_entry(diameter_index_t _diameter_index, coefficient_t _coefficient) {
return std::make_pair(get_diameter(_diameter_index),
make_entry(get_index(_diameter_index), _coefficient));
return std::make_pair(get_diameter(_diameter_index), make_entry(get_index(_diameter_index), _coefficient));
}

template <typename Entry> struct greater_diameter_or_smaller_index {
Expand Down Expand Up @@ -230,18 +219,16 @@ template <typename DistanceMatrix> class rips_filtration_comparator {
get_simplex_vertices(index, dim, dist.size(), binomial_coeff, vertices.begin());

for (index_t i = 0; i <= dim; ++i)
for (index_t j = 0; j < i; ++j) {
diam = std::max(diam, dist(vertices[i], vertices[j]));
}
for (index_t j = 0; j < i; ++j) { diam = std::max(diam, dist(vertices[i], vertices[j])); }
return diam;
}

bool operator()(const index_t a, const index_t b) const {
assert(a < binomial_coeff(dist.size(), dim + 1));
assert(b < binomial_coeff(dist.size(), dim + 1));

return greater_diameter_or_smaller_index<diameter_index_t>()(
diameter_index_t(diameter(a), a), diameter_index_t(diameter(b), b));
return greater_diameter_or_smaller_index<diameter_index_t>()(diameter_index_t(diameter(a), a),
diameter_index_t(diameter(b), b));
}

template <typename Entry> bool operator()(const Entry& a, const Entry& b) const {
Expand All @@ -255,8 +242,7 @@ class simplex_coboundary_enumerator {
const binomial_coeff_table& binomial_coeff;

public:
simplex_coboundary_enumerator(index_t _idx, index_t _dim, index_t _n,
const binomial_coeff_table& _binomial_coeff)
simplex_coboundary_enumerator(index_t _idx, index_t _dim, index_t _n, const binomial_coeff_table& _binomial_coeff)
: idx(_idx), modified_idx(_idx), v(_n - 1), k(_dim + 1), binomial_coeff(_binomial_coeff) {}

bool has_next() {
Expand All @@ -271,8 +257,7 @@ class simplex_coboundary_enumerator {
}

std::pair<entry_t, index_t> next() {
auto result =
std::make_pair(make_entry(modified_idx + binomial_coeff(v, k + 1), k & 1 ? -1 : 1), v);
auto result = std::make_pair(make_entry(modified_idx + binomial_coeff(v, k + 1), k & 1 ? -1 : 1), v);
--v;
return result;
}
Expand Down Expand Up @@ -323,16 +308,14 @@ template <> void compressed_distance_matrix<UPPER_TRIANGULAR>::init_rows() {
}
}

template <>
value_t compressed_distance_matrix<UPPER_TRIANGULAR>::operator()(const index_t i,
const index_t j) const {
return i == j ? 0 : rows[std::min(i, j)][std::max(i, j)];
template <> value_t compressed_distance_matrix<UPPER_TRIANGULAR>::operator()(index_t i, index_t j) const {
std::tie(i, j) = std::minmax(i, j);
return i == j ? 0 : rows[i][j];
}

template <>
value_t compressed_distance_matrix<LOWER_TRIANGULAR>::operator()(const index_t i,
const index_t j) const {
return i == j ? 0 : rows[std::max(i, j)][std::min(i, j)];
template <> value_t compressed_distance_matrix<LOWER_TRIANGULAR>::operator()(index_t i, index_t j) const {
std::tie(i, j) = std::minmax(i, j);
return i == j ? 0 : rows[j][i];
}

typedef compressed_distance_matrix<LOWER_TRIANGULAR> compressed_lower_distance_matrix;
Expand All @@ -345,9 +328,9 @@ class euclidean_distance_matrix {
euclidean_distance_matrix(std::vector<std::vector<value_t>>&& _points) : points(_points) {}

value_t operator()(const index_t i, const index_t j) const {
return std::sqrt(std::inner_product(
points[i].begin(), points[i].end(), points[j].begin(), value_t(), std::plus<value_t>(),
[](value_t u, value_t v) { return (u - v) * (u - v); }));
return std::sqrt(std::inner_product(points[i].begin(), points[i].end(), points[j].begin(), value_t(),
std::plus<value_t>(),
[](value_t u, value_t v) { return (u - v) * (u - v); }));
}

size_t size() const { return points.size(); }
Expand Down Expand Up @@ -436,17 +419,15 @@ template <typename ValueType> class compressed_sparse_matrix {
}
};

template <typename Heap>
void push_entry(Heap& column, index_t i, coefficient_t c, value_t diameter) {
template <typename Heap> void push_entry(Heap& column, index_t i, coefficient_t c, value_t diameter) {
entry_t e = make_entry(i, c);
column.push(std::make_pair(diameter, e));
}

template <typename Comparator>
void assemble_columns_to_reduce(std::vector<diameter_index_t>& columns_to_reduce,
hash_map<index_t, index_t>& pivot_column_index,
const Comparator& comp, index_t dim, index_t n, value_t threshold,
const binomial_coeff_table& binomial_coeff) {
hash_map<index_t, index_t>& pivot_column_index, const Comparator& comp, index_t dim,
index_t n, value_t threshold, const binomial_coeff_table& binomial_coeff) {
index_t num_simplices = binomial_coeff(n, dim + 2);

columns_to_reduce.clear();
Expand Down Expand Up @@ -476,9 +457,8 @@ void assemble_columns_to_reduce(std::vector<diameter_index_t>& columns_to_reduce
}

template <typename DistanceMatrix, typename ComparatorCofaces, typename Comparator>
void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
hash_map<index_t, index_t>& pivot_column_index, const DistanceMatrix& dist,
const ComparatorCofaces& comp, const Comparator& comp_prev, index_t dim,
void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce, hash_map<index_t, index_t>& pivot_column_index,
const DistanceMatrix& dist, const ComparatorCofaces& comp, const Comparator& comp_prev, index_t dim,
index_t n, value_t threshold, coefficient_t modulus,
const std::vector<coefficient_t>& multiplicative_inverse,
const binomial_coeff_table& binomial_coeff) {
Expand All @@ -502,8 +482,7 @@ void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
auto column_to_reduce = columns_to_reduce[i];

#ifdef ASSEMBLE_REDUCTION_MATRIX
std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>,
smaller_index<diameter_entry_t>>
std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, smaller_index<diameter_entry_t>>
reduction_column;
#endif

Expand All @@ -516,8 +495,8 @@ void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
#ifdef INDICATE_PROGRESS
if ((i + 1) % 1000 == 0)
std::cout << "\033[K"
<< "reducing column " << i + 1 << "/" << columns_to_reduce.size()
<< " (diameter " << diameter << ")" << std::flush << "\r";
<< "reducing column " << i + 1 << "/" << columns_to_reduce.size() << " (diameter " << diameter
<< ")" << std::flush << "\r";
#endif

index_t j = i;
Expand Down Expand Up @@ -565,13 +544,11 @@ void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
assert(simplex_diameter == comp_prev.diameter(get_index(simplex)));

#ifdef ASSEMBLE_REDUCTION_MATRIX
reduction_column.push(
make_diameter_entry(simplex_diameter, get_index(simplex), simplex_coefficient));
reduction_column.push(make_diameter_entry(simplex_diameter, get_index(simplex), simplex_coefficient));
#endif

vertices.clear();
get_simplex_vertices(get_index(simplex), dim, n, binomial_coeff,
std::back_inserter(vertices));
get_simplex_vertices(get_index(simplex), dim, n, binomial_coeff, std::back_inserter(vertices));

coface_entries.clear();
simplex_coboundary_enumerator cofaces(get_index(simplex), dim, n, binomial_coeff);
Expand All @@ -581,14 +558,11 @@ void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
index_t covertex = coface_descriptor.second;
index_t coface_index = get_index(coface);
value_t coface_diameter = simplex_diameter;
for (index_t v : vertices) {
coface_diameter = std::max(coface_diameter, dist(v, covertex));
}
for (index_t v : vertices) { coface_diameter = std::max(coface_diameter, dist(v, covertex)); }
assert(comp.diameter(coface_index) == coface_diameter);

if (coface_diameter <= threshold) {
coefficient_t coface_coefficient =
get_coefficient(coface) * get_coefficient(simplex) * factor;
coefficient_t coface_coefficient = get_coefficient(coface) * get_coefficient(simplex) * factor;
coface_coefficient %= modulus;
if (coface_coefficient < 0) coface_coefficient += modulus;
assert(coface_coefficient >= 0);
Expand Down Expand Up @@ -660,8 +634,7 @@ void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
#else
const coefficient_t coefficient = 1;
#endif
reduction_matrix.push_back(
make_diameter_entry(get_diameter(e), index, coefficient));
reduction_matrix.push_back(make_diameter_entry(get_diameter(e), index, coefficient));
}
#else
#ifdef USE_COEFFICIENTS
Expand All @@ -678,13 +651,7 @@ void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
#endif
}

enum file_format {
LOWER_DISTANCE_MATRIX,
UPPER_DISTANCE_MATRIX,
DISTANCE_MATRIX,
POINT_CLOUD,
DIPHA
};
enum file_format { LOWER_DISTANCE_MATRIX, UPPER_DISTANCE_MATRIX, DISTANCE_MATRIX, POINT_CLOUD, DIPHA };

template <typename T> T read(std::istream& s) {
T result;
Expand Down Expand Up @@ -712,8 +679,7 @@ compressed_lower_distance_matrix read_point_cloud(std::istream& input_stream) {

index_t n = eucl_dist.size();

std::cout << "point cloud with " << n << " points in dimension "
<< eucl_dist.points.front().size() << std::endl;
std::cout << "point cloud with " << n << " points in dimension " << eucl_dist.points.front().size() << std::endl;

std::vector<value_t> distances;

Expand Down Expand Up @@ -810,16 +776,12 @@ void print_usage_and_exit(int exit_code) {
<< "Options:" << std::endl
<< std::endl
<< " --help print this screen" << std::endl
<< " --format use the specified file format for the input. Options are:"
<< std::endl
<< " lower-distance (lower triangular distance matrix; default)"
<< std::endl
<< " upper-distance (upper triangular distance matrix)"
<< std::endl
<< " --format use the specified file format for the input. Options are:" << std::endl
<< " lower-distance (lower triangular distance matrix; default)" << std::endl
<< " upper-distance (upper triangular distance matrix)" << std::endl
<< " distance (full distance matrix)" << std::endl
<< " point-cloud (point cloud in Euclidean space)" << std::endl
<< " dipha (distance matrix in DIPHA file format)"
<< std::endl
<< " dipha (distance matrix in DIPHA file format)" << std::endl
<< " --dim <k> compute persistent homology up to dimension <k>" << std::endl
<< " --threshold <t> compute Rips complexes up to diameter <t>" << std::endl
#ifdef USE_COEFFICIENTS
Expand Down Expand Up @@ -898,8 +860,8 @@ int main(int argc, char** argv) {

std::cout << "distance matrix with " << n << " points" << std::endl;

auto result = std::minmax_element(dist.distances.begin(), dist.distances.end());
std::cout << "value range: [" << *result.first << "," << *result.second << "]" << 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);

Expand All @@ -917,11 +879,10 @@ int main(int argc, char** argv) {
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, dist, comp, comp_prev, dim, n,
threshold, modulus, multiplicative_inverse, binomial_coeff);
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);
assemble_columns_to_reduce(columns_to_reduce, pivot_column_index, comp, dim, n, threshold, binomial_coeff);
}
}

0 comments on commit a76f7e8

Please sign in to comment.