if(NOT CMAKE_SYSTEM_NAME STREQUAL "Windows")
- set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -ansi -pedantic")
+ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -std=c++98 -pedantic")
# I would like to use -pedantic, but VTK devs are not sufficiently
# pedantic
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Werror")
/** 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) {
+ void addColumn(index_type c, const column_type &col) {
assert(c < width());
// Flip all entries that are set in 'col'.
}
#ifndef NDEBUG
+ // C++11 would have is_sorted
for (unsigned i=1; i < col.size(); i++)
assert(col[i-1] < col[i]);
#endif
/** Add the row-vector 'row' to the r-th row. Note that 'row'
* actually contains the list of column-indices that are 1. */
- void add_row(index_type r, const row_type &row) {
+ void addRow(index_type r, const row_type &row) {
assert(r < height());
// Flip all entries that are set in 'row'.
}
#ifndef NDEBUG
+ // C++11 would have is_sorted
for (unsigned i=1; i < row.size(); i++)
assert(row[i-1] < row[i]);
#endif
bool operator!=(const BooleanColRowMatrix<IT> &m) const {
return !(*this == m);
}
+
+ /** Multiply with matrix b from the right */
+ template<class D>
+ BooleanColRowMatrix<IT> operator*(const BooleanColMatrix_base<IT, D> &b) {
+ BooleanColRowMatrix<IT> c(size());
+ multiply_matrix(c, *this, b);
+ return c;
+ }
};
+
+/** Counts the number of common elements in two sorted vectors. Equal to
+ * counting the number of elements given by std::set_intersection. */
+template <class InputIterator1, class InputIterator2>
+size_t count_set_intersection (InputIterator1 first1, InputIterator1 last1, InputIterator2 first2, InputIterator2 last2)
+{
+ size_t count = 0;
+
+ // As long as we did not either end, look for common elements
+ while (first1!=last1 && first2!=last2)
+ {
+ if (*first1 < *first2)
+ ++first1;
+ else if (*first2 < *first1)
+ ++first2;
+ // Element is common to both
+ else {
+ ++count;
+ ++first1;
+ ++first2;
+ }
+ }
+ return count;
+}
+
+/** Multiply a*b and save the product in 'result'. It is assumed that 'result' is intially empty and has appropriate size. */
+template<class IT, class D, class RT>
+void multiply_matrix(RT &result, const BooleanRowMatrix_base<IT, D> &a, const BooleanColMatrix_base<IT, D> &b) {
+ assert(a.height() == b.width());
+
+ for (unsigned r=0; r < a.height(); ++r) {
+ const typename BooleanRowMatrix_base<IT, D>::row_type &row = a.getRow(r);
+ for (unsigned c=0; c < b.width(); ++c) {
+ const typename BooleanColMatrix_base<IT, D>::column_type &col = b.getColumn(c);
+ if (count_set_intersection(row.begin(), row.end(), col.begin(), col.end()) % 2 == 1)
+ result.set(r, c, true);
+ }
+ }
+}
+
+template<class MT>
+MT unitMatrix(size_t size) {
+ MT mat(size);
+ for (unsigned i=0; i < size; ++i)
+ mat.set(i, i, true);
+ return mat;
+}
+
template<class IT>
std::ostream& operator<<(std::ostream &os, const BooleanColRowMatrix<IT> &mat) {
for (unsigned r=0; r < mat.size(); ++r) {
--- /dev/null
+#ifndef matrixreduction_h_aPeinoXeiFeethuz
+#define matrixreduction_h_aPeinoXeiFeethuz
+
+#include <iostream>
+
+#include "booleanmatrix.h"
+
+
+namespace libstick {
+
+/** Takes the boundary matrix b and transforms it inplace into a reduced matrix
+ * b. The matrix v is required to be a unit matrix having the same dimensions
+ * as b. It is maintained in order to leave the product b*inverse(v) as an
+ * invariant. Hence, the resulting b is equal to the product of the boundary
+ * matrix times v. */
+template<class IT, class D>
+void reduceBoundaryMatrix(BooleanColRowMatrix<IT> &b, BooleanColMatrix_base<IT, D> &v) {
+ assert(b.size() == v.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 BooleanColRowMatrix<IT>::colbase::column_type &col = b.getColumn(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 BooleanColRowMatrix<IT>::rowbase::row_type row(b.getRow(r).size());
+ copy(b.getRow(r).begin(), b.getRow(r).end(), row.begin());
+ assert(row.size() >= 1);
+
+ // Get rid of 1s at that row right of column c.
+ typename BooleanColRowMatrix<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.addColumn(*it, col);
+ v.addColumn(*it, v.getColumn(c));
+ assert(!b.get(r, *it));
+ }
+ }
+}
+
+}
+
+#endif
#include <cstring>
#include <algorithm>
#include <vector>
+#include <limits>
#include <iostream>
/** A simplicial complex is a std::vector of simplices such that each face is
* also part of the complex. Every simplex has dimension at most MAXDIM. The
* indices of simplices resp. their faces are of type IT. To each simplex a
- * value is assigend, which is of type VT. */
-template<size_t MAXDIM, class IT=uint32_t, class VT=double>
+ * value is assigend, which is of type VT. When a SimplicialComplex is
+ * instantiated, a single (-1) dimensional simplex is automatically created.
+ * Each 0-dimensional simplex automatically has this simplex as its face.
+ * Consequently, the innner class SimplexOrder gives the extended boundary
+ * matrix. */
+template<int MAXDIM, class IT=uint32_t, class VT=double>
class SimplicialComplex {
public:
/** A simplex of the complex. */
struct Simplex {
/** Dimension of the simplex. */
- size_t dim;
+ int dim;
/** The indices of the faces of the simplex. */
index_type faces[MAXDIM+1];
/** The value of the simplex. */
value_type value;
/** Create a new simplex with dimension 'dim', (dim+1)-faces and
- * its value. */
- static Simplex create(size_t dim, index_type* faces, value_type value) {
- assert(dim <= MAXDIM);
+ * its value. If simpley is 0-dimensional, its face is
+ * automatically set to one (-1)-dimensional simplex. */
+ static Simplex create(int dim, index_type* faces, value_type value) {
+ assert(0 <= dim && dim <= MAXDIM);
Simplex s;
s.dim = dim;
s.value = value;
- memcpy(s.faces, faces, face_count_bydim(dim)*sizeof(index_type));
+
+ if (dim > 0)
+ memcpy(s.faces, faces, faceCountByDim(dim)*sizeof(index_type));
+ else
+ s.faces[0] = 0;
+
+ return s;
+ }
+
+ /** Create a (-1)-dimensional simplex. It has the lowest possible value. */
+ static Simplex create_minusonedim_simplex() {
+ Simplex s;
+
+ s.dim = -1;
+ s.faces[0] = 0;
+ s.value = std::numeric_limits<VT>::has_infinity
+ ? -std::numeric_limits<VT>::infinity()
+ : std::numeric_limits<VT>::min();
return s;
}
/** Get number of faces. */
- size_t face_count() const {
- return face_count_bydim(dim);
+ size_t faceCount() const {
+ return faceCountByDim(dim);
}
/** Get number of faces of a dim-dimensional simplex. */
- static size_t face_count_bydim(size_t dim) {
- if (dim == 0)
- return 0;
+ static size_t faceCountByDim(int dim) {
+ assert(-1 <= dim && dim <= MAXDIM);
return dim + 1;
}
};
assert(size() == c.size());
for (unsigned i=0; i < size(); ++i)
- for (unsigned f=0; f < getSimplex(i).face_count(); ++f)
+ for (unsigned f=0; f < getSimplex(i).faceCount(); ++f)
if (revorder[getSimplex(i).faces[f]] >= i)
return false;
BoundaryMatrix mat(size());
for (unsigned c=0; c < size(); ++c)
- for(unsigned r=0; r < getSimplex(c).face_count(); ++r)
+ for(unsigned r=0; r < getSimplex(c).faceCount(); ++r)
mat.set(revorder[getSimplex(c).faces[r]], c, true);
return mat;
};
public:
+ SimplicialComplex() {
+ // Add the one minus-one dimensional simplex
+ addSimplex(Simplex::create_minusonedim_simplex());
+ }
+
/** Return number of simplices. */
size_t size() const {
return simplices.size();
/** Add a simplex to the complex. The dimension of the faces must be
* dim-1, and they must already be part of the complex. */
- void addSimplex(size_t dim, index_type* faces, value_type value) {
+ void addSimplex(int dim, index_type* faces, value_type value) {
addSimplex(Simplex::create(dim, faces, value));
}
* dim-1, and they must already be part of the complex. */
void addSimplex(Simplex s) {
// Check requirements for faces
- for (unsigned i=0; i < s.face_count(); ++i) {
+ for (unsigned i=0; i < s.faceCount(); ++i) {
// Faces are already in complex.
assert(s.faces[i] < size());
// Faces have dimension dim-1
/** Return true iff for each simplex i with dimension dim it holds that
* the faces of i are contained in the complex and have dimension dim-1. */
bool isComplex() const {
- typename std::vector<Simplex>::const_iterator it = ++simplices.begin();
for (unsigned i=0; i < size(); ++i) {
const Simplex &s = simplices[i];
- for (unsigned f=0; f < s.face_count(); ++f) {
+ for (unsigned f=0; f < s.faceCount(); ++f) {
if (s.faces[f] >= size())
return false;
typename std::vector<Simplex>::const_iterator it = ++simplices.begin();
for (; it != simplices.end(); ++it)
- for (unsigned f=0; f < it->face_count(); ++f)
+ for (unsigned f=0; f < it->faceCount(); ++f)
if (simplices[it->faces[f]].value > it->value)
return false;
col.push_back(0);
// Add column and test for values
- mat.add_column(5, col);
+ mat.addColumn(5, col);
for (unsigned c=0; c < size; ++c) {
if (c==5) {
for (unsigned r=0; r < size; ++r )
row.push_back(3);
// Add row and test for values
- mat.add_row(5, row);
+ mat.addRow(5, row);
for (unsigned r=0; r < size; ++r) {
if (r==5) {
for (unsigned c=0; c < size; ++c )
#include <cpptest-suite.h>
#include <libstick-0.1/simplicialcomplex.h>
+#include <libstick-0.1/matrixreduction.h>
using namespace libstick;
TEST_ADD(SimplicialComplexTestSuite::test_isOrderFiltration);
TEST_ADD(SimplicialComplexTestSuite::test_isOrderMonotone);
TEST_ADD(SimplicialComplexTestSuite::test_boundaryMatrix);
+ TEST_ADD(SimplicialComplexTestSuite::test_matrixreduction);
}
protected:
const unsigned num = 11;
scomplex::Simplex ss[num] = {
// dimension, faces, value...
- {0, {0, 0, 0}, 0},
- {0, {0, 0, 0}, 1},
- {0, {0, 0, 0}, 2},
- {0, {0, 0, 0}, 3},
- {1, {0, 1, 0}, 4},
- {1, {1, 2, 0}, 5},
- {1, {2, 3, 0}, 6},
- {1, {3, 0, 0}, 7},
- {1, {0, 2, 0}, 8},
- {2, {6, 7, 8}, 9},
- {2, {4, 5, 8}, 10}
+ {0, {}, 1},
+ {0, {}, 2},
+ {0, {}, 3},
+ {0, {}, 4},
+ {1, {1, 2}, 5},
+ {1, {2, 3}, 6},
+ {1, {3, 4}, 7},
+ {1, {4, 1}, 8},
+ {1, {1, 3}, 9},
+ {2, {7, 8, 9}, 10},
+ {2, {5, 6, 9}, 11}
};
// This is o1 This is o2 This is o3(b)
// (value = index) (values) (values)
//
- // 0 ----4---- 1 0 ----4---- 1 0 ----4---- 1
- // |\ | |\ | |\ |
- // | \ 10 | | \ 10 | | \ 12 |
- // | \ | | \ | | \ |
- // | \ | | \ | | \ |
- // 7 8 5 7 8 11 7 8 11
- // | \ | | \ | | \ |
- // | 9 \ | | 9 \ | | 9 \ |
- // | \ | | \ | | \ |
- // | \| | \| | \|
- // 3 ----6---- 2 3 ----6---- 2 3 ----6---- 2
+ // 1 ----5---- 2 1 ----5---- 2 1 ----5---- 2
+ // |\ | |\ | |\ |
+ // | \ 11 | | \ 11 | | \ 13 |
+ // | \ | | \ | | \ |
+ // | \ | | \ | | \ |
+ // 8 9 6 8 9 12 8 9 12
+ // | \ | | \ | | \ |
+ // | 10 \ | | 10 \ | | 10 \ |
+ // | \ | | \ | | \ |
+ // | \| | \| | \|
+ // 4 ----7---- 3 4 ----7---- 3 4 ----7---- 3
// Build the complex
for (unsigned i=0; i<num; ++i)
c1.addSimplex(ss[i]);
c2 = c1;
- c2.simplices[5].value = 11;
+ c2.simplices[6].value = 12;
c3 = c2;
- c3.simplices[10].value = 12;
+ c3.simplices[11].value = 13;
o1.reset();
o2.reset();
void test_boundaryMatrix() {
bm mat1 = o1.getBoundaryMatrix();
bm mat1e(c1.size());
+ mat1e.set(0, 1, true);
+ mat1e.set(0, 2, true);
+ mat1e.set(0, 3, true);
mat1e.set(0, 4, true);
- mat1e.set(1, 4, true);
mat1e.set(1, 5, true);
mat1e.set(2, 5, true);
mat1e.set(2, 6, true);
mat1e.set(3, 6, true);
mat1e.set(3, 7, true);
- mat1e.set(0, 7, true);
- mat1e.set(0, 8, true);
- mat1e.set(2, 8, true);
- mat1e.set(6, 9, true);
- mat1e.set(7, 9, true);
- mat1e.set(8, 9, true);
- mat1e.set(4, 10, true);
- mat1e.set(5, 10, true);
+ mat1e.set(4, 7, true);
+ mat1e.set(4, 8, true);
+ mat1e.set(1, 8, true);
+ mat1e.set(1, 9, true);
+ mat1e.set(3, 9, true);
+ mat1e.set(7, 10, true);
mat1e.set(8, 10, true);
+ mat1e.set(9, 10, true);
+ mat1e.set(5, 11, true);
+ mat1e.set(6, 11, true);
+ mat1e.set(9, 11, true);
TEST_ASSERT(mat1 == mat1e);
bm mat3b = o3b.getBoundaryMatrix();
bm mat3be(c1.size());
+ mat3be.set(0, 1, true);
+ mat3be.set(0, 2, true);
+ mat3be.set(0, 3, true);
mat3be.set(0, 4, true);
- mat3be.set(1, 4, true);
+ mat3be.set(1, 5, true);
mat3be.set(2, 5, true);
- mat3be.set(3, 5, true);
mat3be.set(3, 6, true);
- mat3be.set(0, 6, true);
- mat3be.set(0, 7, true);
- mat3be.set(2, 7, true);
- mat3be.set(5, 8, true);
- mat3be.set(6, 8, true);
- mat3be.set(7, 8, true);
- mat3be.set(1, 9, true);
- mat3be.set(2, 9, true);
- mat3be.set(4, 10, true);
- mat3be.set(7, 10, true);
- mat3be.set(9, 10, true);
+ mat3be.set(4, 6, true);
+ mat3be.set(4, 7, true);
+ mat3be.set(1, 7, true);
+ mat3be.set(1, 8, true);
+ mat3be.set(3, 8, true);
+ mat3be.set(6, 9, true);
+ mat3be.set(7, 9, true);
+ mat3be.set(8, 9, true);
+ mat3be.set(2, 10, true);
+ mat3be.set(3, 10, true);
+ mat3be.set(5, 11, true);
+ mat3be.set(8, 11, true);
+ mat3be.set(10, 11, true);
+ TEST_ASSERT(mat3b == mat3be);
//std::cout << mat3b << std::endl << std::endl;
+ //std::cout << mat3be << std::endl << std::endl;
//std::cout << ((bm::colbase) mat3b) << std::endl << std::endl;
//std::cout << ((bm::rowbase) mat3b) << std::endl << std::endl;
+ }
- TEST_ASSERT(mat3b == mat3be);
+ void test_matrixreduction() {
+ bm b = o3b.getBoundaryMatrix();
+ bm v = unitMatrix<bm>(b.size());
+
+ bm b_orig = b;
+ reduceBoundaryMatrix(b, v);
+ TEST_ASSERT(b == b_orig*v);
+
+ //std::cout << std::endl << "after reduce: " << std::endl;
+ //std::cout << std::endl << b << std::endl;
+ //std::cout << std::endl << v << std::endl;
+ //std::cout << std::endl << (b_orig*v) << std::endl;
}
};
projects / libstick.git / commitdiff