- Use colmatrix as boundary matrix.
- Use faster set_symmetric_difference() for colmatrix::add_column().
- Move diagram sanity checks to unit tests.
- Make lowestones a col-index vector indexed by row-indices.
- Test lowestones in unit tests, which are independt of the reduction
algorithm.
namespace libstick {
+template<class T>
+std::ostream& operator<<(std::ostream &, const std::vector<T> &);
+
+
/** The base class of boolean_colmatrix and boolean_colrowmatrix which implements
* the common logic of both. */
template<class IT, class D>
}
public:
- /** Get height resp. width of the matrix. */
+ /** A casting constructor for any colmatrix with same entries type. */
+ template<class D2>
+ boolean_colmatrix_base(const boolean_colmatrix_base<IT, D2> &mat) :
+ cols(mat.get_columns()) {
+ }
+
+ /** Get width of the matrix. */
size_t width() const {
return cols.size();
}
+ /** Get height of the matrix, i.e. maximum row-index + 1 among all columns. */
+ size_t height() const {
+ IT h = 0;
+ for (unsigned c=0; c < width(); ++c)
+ if (cols[c].size() > 0)
+ h = std::max(h, cols[c].back());
+ return h+1;
+ }
+
/** Get the matrix entry at row 'r' and column 'c'. */
bool get(index_type r, index_type c) const {
assert(c < width());
return cols[c];
}
+ /** Get all columns */
+ const std::vector<column_type>& get_columns() const {
+ return cols;
+ }
+
/** Add the column-vector 'col' to the c-th column. Note that 'col'
* actually contains the list of row-indices that are 1. */
void add_column(index_type c, const column_type &col) {
}
/** Two matrices are equal iff they have the same entries */
- bool operator==(const boolean_colmatrix_base<IT, D> &m) const {
- return cols == m.cols;
+ template<class D2>
+ bool operator==(const boolean_colmatrix_base<IT, D2> &m) const {
+ return cols == m.get_columns();
}
/** Two matrices are equal iff they have the same entries */
- bool operator!=(const boolean_colmatrix_base<IT, D> &m) const {
+ template<class D2>
+ bool operator!=(const boolean_colmatrix_base<IT, D2> &m) const {
return !(*this == m);
}
#endif
}
- private:
+ protected:
/** The matrix is the set of columns. */
std::vector<column_type> cols;
public:
typedef IT index_type;
typedef boolean_colmatrix_base<IT, boolean_colmatrix<IT> > base;
+ typedef typename base::column_type column_type;
/** Create a matrix with 'width' columns, initalized with zero entries. */
boolean_colmatrix(size_t columns) :
void _set(index_type r, index_type c, bool value) {
base::_set(r, c, value);
}
+
+ /** A faster implementation of boolean_colmatrix_base::add_column().
+ * Assumes that 'col' is sorted. */
+ void add_column(index_type c, const column_type &col) {
+ assert(c < base::width());
+
+#ifndef NDEBUG
+ for (unsigned i=1; i < col.size(); ++i)
+ assert(col[i-1] < col[i]);
+#endif
+
+ // The original column
+ column_type &orig_col = base::cols[c];
+
+ // Make target column large enough
+ const size_t maxsize = orig_col.size() + col.size();
+ if (tcol.size() < maxsize)
+ tcol.resize(maxsize);
+
+ // Compute symmetric difference
+ typename column_type::iterator it = std::set_symmetric_difference(
+ orig_col.begin(), orig_col.end(), col.begin(), col.end(), tcol.begin());
+
+ // Copy back to the original column
+ orig_col.resize(it - tcol.begin());
+ std::copy(tcol.begin(), it, orig_col.begin());
+ }
+
+ private:
+ /** A temporary container to speed up add_column() */
+ column_type tcol;
};
}
public:
- /** Get height resp. width of the matrix. */
+ /** Get height of the matrix. */
size_t height() const {
return rows.size();
}
#endif
}
- private:
+ protected:
derived* get_derived() {
return static_cast<derived*>(this);
}
rowbase(size) {
}
+ /** Casting a colmatrix into a colrow matrix */
+ template<class D>
+ boolean_colrowmatrix(const boolean_colmatrix_base<IT, D>& mat) :
+ colbase(std::max(mat.width(), mat.height())),
+ rowbase(std::max(mat.width(), mat.height())) {
+ for (unsigned c=0; c < mat.width(); ++c) {
+ const typename colbase::column_type &col = mat.get_column(c);
+ for (unsigned i=0; i < col.size(); ++i)
+ set(col[i], c, true);
+ }
+ for (unsigned r=0; r < size(); ++r) {
+ const typename rowbase::row_type &row = rowbase::get_row(r);
+ for (unsigned i=0; i < row.size(); ++i)
+ assert(get(r, row[i]) == true);
+ }
+ }
+
/** Override implementation. */
void _set(index_type r, index_type c, bool value) {
colbase::_set(r, c, value);
return colbase::operator==(m);
}
+ /** Two matrices are equal iff they have the same entries */
+ template<class D>
+ bool operator==(const boolean_colmatrix_base<IT, D> &m) const {
+ return colbase::operator==(m);
+ }
+
/** Two matrices are equal iff they have the same entries */
bool operator!=(const boolean_colrowmatrix<IT> &m) const {
return !(*this == m);
return os << m;
}
+template<class T>
+std::ostream& operator<<(std::ostream& os, const std::vector<T> &vec) {
+ os << "[";
+
+ typename std::vector<T>::const_iterator it = vec.begin();
+ while ( it != vec.end()) {
+ os << *it;
+ if (++it != vec.end())
+ os << " ";
+ }
+
+ os << "]";
+ return os;
+}
+
}
#endif
namespace libstick {
-template<class T>
-std::ostream& operator<<(std::ostream &, const std::vector<T> &);
-
-
/** Persistence computers various persistent homology information on a
* simplex_order of a simplicial_complex. */
template<int MAXDIM, class IT, class VT>
public:
typedef simplicial_complex<MAXDIM, IT, VT> scomplex;
typedef typename scomplex::simplex_order simplex_order;
- typedef typename simplex_order::boundary_matrix boundary_matrix;
+ typedef boolean_colmatrix<IT> boundary_matrix;
typedef boolean_colmatrix<IT> transformation_matrix;
+ typedef std::vector<IT> lowestones_matrix;
/** A point in a persistence diagram comprises birth- and death-index.
* These are the indices of the simplices of 'order' that give birth
bm = order.get_boundary_matrix();
rbm = bm;
- tm = create_unit_matrix<transformation_matrix>(bm.size());
-
- reduce_boundary_matrix(rbm, tm);
+ tm = create_unit_matrix<transformation_matrix>(bm.width());
+ lowestones = lowestones_matrix(bm.width());
+
+ // Make every column reduced, i.e., it is a zero-column or contains only one 1.
+ for (unsigned c=0; c < rbm.width(); ++c) {
+ //if (c % 100 == 0)
+ //std::cout << "c = " << c << " (" << (100.0*float(c)/rbm.width()) << " %)" << std::endl;
+ // Reduce as long as we need to reduce
+ while (rbm.get_column(c).size() > 0) {
+ // (r, c) is the lowest one of col
+ const IT r = rbm.get_column(c).back();
+ // This column contains its lowest one on the same row
+ const IT c2 = lowestones[r];
+ assert(c2 < c);
+
+ // There is no valid c2, hence we have completely reduced
+ if (c2 == 0)
+ break;
+ // Reduce col by column c2
+ else {
+ assert(rbm.get(r, c));
+ rbm.add_column(c, rbm.get_column(c2));
+ tm.add_column(c, tm.get_column(c2));
+ assert(!rbm.get(r, c));
+ }
+ }
- lowestones = boundary_matrix(rbm.size());
- for (unsigned c=0; c < rbm.size(); ++c)
- if (rbm.get_column(c).size() > 0)
- lowestones.set(rbm.get_column(c).back(), c, true);
+ // A lowest one remained, recall it
+ if (rbm.get_column(c).size() > 0) {
+ const IT r = rbm.get_column(c).back();
+ assert(lowestones[r] == 0);
+ assert(r < c);
+ lowestones[r] = c;
+ }
+ }
}
/** Get the lowest-one matrix of 'order'. */
- const boundary_matrix& get_lowestones_matrix() const {
+ const lowestones_matrix& get_lowestones_matrix() const {
assert(done_matrices);
return lowestones;
}
void interpret_reduction(std::ostream &os) const {
assert(done_matrices);
- for (unsigned c=0; c < bm.size(); ++c) {
+ for (unsigned c=0; c < bm.width(); ++c) {
os << c << ". inserting ";
switch (bm.get_column(c).size()) {
os << std::endl;
os << " boundary: " << bm.get_column(c) << std::endl;
- typename boolean_colrowmatrix<IT>::colbase::column_type col = rbm.get_column(c);
+ typename boolean_colmatrix<IT>::column_type col = rbm.get_column(c);
if (col.size() == 0) {
os << " \033[1;32mbirth\033[0;m of a cycle: " << tm.get_column(c) << std::endl;
} else {
diagrams[d] = diagram();
for (unsigned c=0; c < lowestones.size(); ++c) {
- if (lowestones.get_column(c).size() == 0) {
+ // A cycle was born
+ if (rbm.get_column(c).size() == 0) {
const int dim = order.get_simplex(c).dim;
// Create a diagram point
p.birth = c;
p.death = 0;
// If the class actually dies, get the index of the killing simplex
- if (lowestones.get_row(c).size() > 0)
- p.death = lowestones.get_row(c).back();
+ if (lowestones[c] > 0)
+ p.death = lowestones[c];
_get_persistence_diagram(dim).births.push_back(p);
#ifndef NDEBUG
// Shakespeare principle: A simplex either gives birth or kills, be
// or not be, tertium non datur.
for (unsigned c=0; c < lowestones.size(); ++c) {
- if (lowestones.get_column(c).size() != 0)
+ // A class died with c
+ if (rbm.get_column(c).size() != 0)
assert(c == 0 || deaths.count(c) > 0);
else
assert(births.count(c) > 0);
}
- // Even Jesus died after his birth. Homology classes are nothing better.
- for (int d=-1; d <= MAXDIM; ++d) {
- const diagram &dia = get_persistence_diagram(d);
- for (unsigned i=0; i < dia.births.size(); ++i)
- if (dia.births[i].death > 0)
- assert(dia.births[i].birth < dia.births[i].death);
- }
#endif
}
/** The reduced boundary matrix of 'order' */
boundary_matrix rbm;
/** Only the lowest ones per column of the reduced boundary matrix.
- * Each one at (r, c) corresponds to the death of a class when
- * inserting simplex c, which was born when inserting simplex. */
- boundary_matrix lowestones;
+ * lowestones[i] contains the column-index at which the lowest one at
+ * row i appears, and zero if no such column exists. */
+ lowestones_matrix lowestones;
/** The transformation matrix that transforms bm into rbm, i.e. rbm = bm * tm. */
transformation_matrix tm;
* invariant. Hence, the resulting b is equal to the product of the boundary
* matrix times v. */
template<class IT, class D>
-void reduce_boundary_matrix(boolean_colrowmatrix<IT> &b, boolean_colmatrix_base<IT, D> &v) {
- assert(b.size() == v.width());
+void reduce_boundary_matrix(boolean_colmatrix_base<IT, D> &b, boolean_colmatrix_base<IT, D> &v) {
+ assert(b.width() == v.width());
+
+ // lowestones[i] gives the column whose lowest 1 is at row i. If it contains
+ // zero, there is no such column.
+ std::vector<IT> lowestones(b.width());
// Make every column reduced, i.e., it is a zero-column or contains only one 1.
- for (unsigned c=0; c < b.size(); ++c) {
- const typename boolean_colrowmatrix<IT>::colbase::column_type &col = b.get_column(c);
- if (col.size() == 0)
- continue;
-
- // The column is not a zero-column, take the lowest 1, and reduce the other columns at the row
- IT r = col.back();
- assert(b.get(r, c));
-
- // Get a copy of the r-th row
- typename boolean_colrowmatrix<IT>::rowbase::row_type row(b.get_row(r).size());
- copy(b.get_row(r).begin(), b.get_row(r).end(), row.begin());
- assert(row.size() >= 1);
-
- // Get rid of 1s at that row right of column c.
- typename boolean_colrowmatrix<IT>::row_type::const_iterator it = lower_bound(row.begin(), row.end(), c);
- for(; it != row.end(); ++it) {
- if (*it <= c)
- continue;
-
- assert(b.get(r, *it));
- b.add_column(*it, col);
- v.add_column(*it, v.get_column(c));
- assert(!b.get(r, *it));
+ for (unsigned c=0; c < b.width(); ++c) {
+ //if (c % 100 == 0)
+ //std::cout << "c = " << c << " (" << (100.0*float(c)/b.width()) << " %)" << std::endl;
+ // Reduce as long as we need to reduce
+ while (b.get_column(c).size() > 0) {
+ // (r, c) is the lowest one of col
+ const IT r = b.get_column(c).back();
+ // This column contains its lowest one on the same row
+ const IT c2 = lowestones[r];
+ assert(c2 < c);
+
+ // There is no valid c2, hence we have completely reduced
+ if (c2 == 0)
+ break;
+ // Reduce col by column c2
+ else {
+ assert(b.get(r, c));
+ b.add_column(c, b.get_column(c2));
+ v.add_column(c, v.get_column(c2));
+ assert(!b.get(r, c));
+ }
}
- }
-}
-
-template<class T>
-std::ostream& operator<<(std::ostream& os, const std::vector<T> &vec) {
- os << "[";
- typename std::vector<T>::const_iterator it = vec.begin();
- while ( it != vec.end()) {
- os << *it;
- if (++it != vec.end())
- os << " ";
+ // A lowest one remained, recall it
+ if (b.get_column(c).size() > 0) {
+ const IT r = b.get_column(c).back();
+ assert(lowestones[r] == 0);
+ assert(r < c);
+ lowestones[r] = c;
+ }
}
-
- os << "]";
- return os;
}
class simplex_order {
public:
- typedef boolean_colrowmatrix<IT> boundary_matrix;
+ typedef boolean_colmatrix<IT> boundary_matrix;
/** Create a standard order of the complex c, i.e., the identity permutation. */
simplex_order(const simplcompltype &c) :
for (unsigned c=0; c < size; ++c)
TEST_ASSERT(mat.get_column(c).size() == 0);
+ col.push_back(0);
col.push_back(3);
- col.push_back(14);
col.push_back(8);
- col.push_back(0);
+ col.push_back(14);
// Add column and test for values
mat.add_column(5, col);
sorder so(s);
assert(so.is_filtration());
- typedef sorder::boundary_matrix bmatrix;
+ typedef boolean_colrowmatrix<uint32_t> bmatrix;
bmatrix bm = so.get_boundary_matrix();
// Check for the right vertex incidences
typedef simplicial_complex<3, uint32_t, double> scomplex;
typedef persistence<3, uint32_t, double> pers;
typedef pers::boundary_matrix bm;
+ typedef pers::lowestones_matrix lom;
typedef pers::transformation_matrix tm;
bool setupcalled;
TEST_ADD(persistence_TestSuite::test_matrix_reduction);
TEST_ADD(persistence_TestSuite::test_betti_numbers);
TEST_ADD(persistence_TestSuite::test_lowestones);
+ TEST_ADD(persistence_TestSuite::test_diagram);
}
protected:
const bm &ringbm = ringp.get_boundary_matrix();
const bm &ringrbm = ringp.get_reduced_boundary_matrix();
const tm &ringtm = ringp.get_transformation_matrix();
- TEST_ASSERT(ringrbm == ringbm * ringtm);
+ TEST_ASSERT(ringrbm == boolean_colrowmatrix<uint32_t>(ringbm) * ringtm);
+
pers torusp(otorus);
torusp.compute_matrices();
const bm &torusbm = torusp.get_boundary_matrix();
const bm &torusrbm = torusp.get_reduced_boundary_matrix();
const tm &torustm = torusp.get_transformation_matrix();
- TEST_ASSERT(torusrbm == torusbm * torustm);
-
- uint32_t torusrbme_coords[][2] = {
- {0, 1}, {2, 7}, {3, 7}, {1, 8}, {2, 8}, {4, 10}, {5, 10}, {4, 11}, {6, 11}, {1, 13},
- {4, 13}, {10, 15}, {13, 15}, {14, 15}, {9, 19}, {10, 19}, {13, 19}, {17, 19}, {7, 20},
- {16, 20}, {18, 20}, {7, 21}, {9, 21}, {10, 21}, {11, 21}, {13, 21}, {16, 21}, {12, 23},
- {13, 23}, {22, 23}, {7, 24}, {8, 24}, {9, 24}, {10, 24}, {11, 24}, {12, 24}, {25, 28},
- {26, 28}, {25, 29}, {27, 29}, {1, 31}, {25, 31}, {28, 33}, {31, 33}, {32, 33}, {30, 35},
- {31, 35}, {34, 35}, {8, 37}, {30, 37}, {31, 37}, {36, 37}, {9, 40}, {28, 40}, {31, 40},
- {39, 40}, {9, 41}, {28, 41}, {29, 41}, {31, 41}, {38, 41}, {7, 42}, {8, 42}, {9, 42},
- {28, 42}, {29, 42}, {30, 42}, {7, 49}, {8, 49}, {9, 49}, {28, 49}, {29, 49}, {43, 49},
- {47, 49}, {10, 50}, {11, 50}, {28, 50}, {29, 50}, {43, 50}, {45, 50}, {29, 51}, {44, 51},
- {48, 51}, {10, 52}, {28, 52}, {43, 52}, {44, 52}, {28, 53}, {43, 53}, {46, 53 }
- };
- bm torusrbme(torusbm.size());
- torusrbme.set_all(torusrbme_coords, sizeof(torusrbme_coords)/(2*sizeof(uint32_t)), true);
- TEST_ASSERT(torusrbme == torusrbm);
+ TEST_ASSERT(torusrbm == boolean_colrowmatrix<uint32_t>(torusbm) * torustm);
+
+ const lom &toruslo = torusp.get_lowestones_matrix();
+ uint32_t torusloe_coords[] = {1, 0, 8, 7, 13, 10, 11, 0, 0, 0, 0,
+ 0, 24, 0, 15, 0, 21, 19, 20, 0, 0, 0, 23, 0, 0, 31, 28, 29, 0,
+ 0, 42, 0, 33, 0, 35, 0, 37, 0, 41, 40, 0, 0, 0, 0, 52, 50, 53,
+ 49, 51, 0, 0, 0, 0, 0, 0};
+ lom torusloe(torusloe_coords, torusloe_coords + sizeof(torusloe_coords)/sizeof(uint32_t));
+ TEST_ASSERT(torusloe == toruslo);
+
pers ballp(oball);
ballp.compute_matrices();
const bm &ballbm = ballp.get_boundary_matrix();
const bm &ballrbm = ballp.get_reduced_boundary_matrix();
const tm &balltm = ballp.get_transformation_matrix();
- TEST_ASSERT(ballrbm == ballbm * balltm);
-
- uint32_t ballrbme_coords[][2] = {
- {0, 1}, {1, 6}, {2, 6}, {1, 7}, {3, 7}, {1, 8}, {4, 8}, {1, 10}, {5, 10}, {6, 14},
- {10, 14}, {11, 14}, {6, 15}, {7, 15}, {10, 15}, {12, 15}, {6, 16}, {7, 16}, {8, 16},
- {10, 16}, {13, 16}, {6, 17}, {7, 17}, {8, 17}, {9, 17}, {6, 19}, {7, 19}, {18, 19},
- {14, 22}, {15, 22}, {19, 22}, {21, 22}, {14, 23}, {15, 23}, {16, 23}, {17, 23}, {19, 23},
- {20, 23 }
- };
- bm ballrbme(ballrbm.size());
- ballrbme.set_all(ballrbme_coords, sizeof(ballrbme_coords)/(2*sizeof(uint32_t)), true);
- TEST_ASSERT(ballrbme == ballrbm);
+ TEST_ASSERT(ballrbm == boolean_colrowmatrix<uint32_t>(ballbm) * balltm);
+
+ const lom &balllo = ballp.get_lowestones_matrix();
+ uint32_t ballloe_coords[] = {1, 0, 6, 7, 8, 10, 0, 0, 0, 17, 0, 14,
+ 15, 16, 0, 0, 0, 0, 19, 0, 23, 22, 0, 0};
+ lom ballloe(ballloe_coords, ballloe_coords + sizeof(ballloe_coords)/sizeof(uint32_t));
+ TEST_ASSERT(ballloe == balllo);
//std::cout << std::endl;
//ballp.interpret_reduction(std::cout);
void test_lowestones_impl(const scomplex::simplex_order& so) {
pers p(so);
p.compute_matrices();
- const bm &lowestones = p.get_lowestones_matrix();
-
- for (unsigned c=0; c < lowestones.size(); ++c)
- assert(lowestones.get_column(c).size() <= 1);
- for (unsigned r=0; r < lowestones.size(); ++r)
- assert(lowestones.get_row(r).size() <= 1);
- for (unsigned c=0; c < lowestones.size(); ++c) {
+ const lom &lowestones = p.get_lowestones_matrix();
+
+ // Check that lowestones contains every column-index only once,
+ // except the zeros.
+ lom cpy = lowestones;
+ sort(cpy.begin(), cpy.end());
+ lom::iterator it = cpy.begin();
+ // Skip the zeros
+ while (it != cpy.end() && *it == 0)
+ ++it;
+ // From now on, no column index appears twice
+ for (; it+1 != cpy.end(); ++it)
+ TEST_ASSERT(*it < *(it+1));
+
+ for (unsigned r=0; r < lowestones.size(); ++r) {
// If (r,c) is a one then
// (i) a cycle dies with c --> row c is zero
// (ii) a cycle is born with r --> column r is zero
//Hence
// (i) the column r is a zero-column
// (i) the row c is a zero-column
- if (lowestones.get_column(c).size() > 0) {
- const uint32_t r = lowestones.get_column(c).back();
- assert(lowestones.get_column(r).size() == 0);
- assert(lowestones.get_row(c).size() == 0);
+
+ const uint32_t c = lowestones[r];
+ if (c > 0) {
+ TEST_ASSERT(p.get_reduced_boundary_matrix().get(r,c) == true);
+ TEST_ASSERT(p.get_reduced_boundary_matrix().get_column(r).size() == 0);
+ TEST_ASSERT(lowestones[c] == 0);
}
}
}
+
+ void test_diagram() {
+ test_diagram_impl(oring);
+ test_diagram_impl(otorus);
+ test_diagram_impl(oball);
+ }
+
+ void test_diagram_impl(const scomplex::simplex_order& so) {
+ pers p(so);
+ p.compute_diagrams();
+
+ // Even Jesus died after his birth. Homology classes are nothing better.
+ for (int d=-1; d <= 3; ++d) {
+ const pers::diagram &dia = p.get_persistence_diagram(d);
+ for (unsigned i=0; i < dia.births.size(); ++i)
+ if (dia.births[i].death > 0)
+ assert(dia.births[i].birth < dia.births[i].death);
+ }
+
+ }
};
#endif
projects / libstick.git / commitdiff