- Move matrix_reduction() to a persistence class.
- Add persistence test suite.
- Add a few helper functions to boolean_matrix classes.
typedef std::vector<index_type> column_type;
typedef D derived;
+ protected:
/** Create a matrix with 'width' columns, initalized with zero entries. */
boolean_colmatrix_base(size_t width) :
cols(width) {
}
+ public:
/** Get height resp. width of the matrix. */
size_t width() const {
return cols.size();
get_derived()->_set(r, c, value);
}
+ /** For each of the 'count'-many (row-index, column-pair) pair in 'indices', set the specific value. */
+ void set_all(index_type indices[][2], size_t count, bool value) {
+ for (unsigned i=0; i < count; ++i)
+ set(indices[i][0], indices[i][1], value);
+ }
+
/** Get the c-th column. */
const column_type& get_column(index_type c) const {
assert(c < width());
return !(*this == m);
}
+ /** Print index pairs of the 1s as { {r, c}, {r, c}, {r, c} } */
+ void print_indexpairs(std::ostream &os) const {
+ bool first=true;
+
+ os << "{";
+ for (unsigned c=0; c < width(); ++c) {
+ const column_type &col = get_column(c);
+ for (typename column_type::const_iterator it = col.begin(); it != col.end(); ++it) {
+ if (first)
+ first = false;
+ else
+ os << ",";
+ os << " {" << *it << ", " << c << "}";
+ }
+ }
+ os << " }";
+ }
+
protected:
/** Set the matrix entry at row 'r' and column 'c'. */
void _set(index_type r, index_type c, bool value) {
typedef std::vector<index_type> row_type;
typedef D derived;
+ protected:
/** Create a matrix with 'height' rows, initalized with zero entries.
* */
boolean_rowmatrix_base(size_t height) :
rows(height) {
}
+ public:
/** Get height resp. width of the matrix. */
size_t height() const {
return rows.size();
get_derived()->_set(r, c, value);
}
+ /** For each of the 'count'-many (row-index, column-pair) pair in 'indices', set the specific value. */
+ void set_all(index_type indices[][2], size_t count, bool value) {
+ for (unsigned i=0; i < count; ++i)
+ set(indices[i][0], indices[i][1], value);
+ }
+
/** Get the r-th row. */
const row_type& get_row(index_type r) const {
assert(r < height());
colbase::set(r, c, value);
}
+ /** For each of the 'count'-many (row-index, column-pair) pair in 'indices', set the specific value. */
+ void set_all(index_type indices[][2], size_t count, bool value) {
+ colbase::set_all(indices, count, value);
+ }
+
/** Get the matrix entry at row 'r' and column 'c'. */
bool get(index_type r, index_type c) const {
assert(colbase::get(r, c) == rowbase::get(r, c));
/** Multiply with matrix b from the right */
template<class D>
- boolean_colrowmatrix<IT> operator*(const boolean_colmatrix_base<IT, D> &b) {
+ boolean_colrowmatrix<IT> operator*(const boolean_colmatrix_base<IT, D> &b) const {
boolean_colrowmatrix<IT> c(size());
multiply_matrix(c, *this, b);
return c;
}
/** 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 boolean_rowmatrix_base<IT, D> &a, const boolean_colmatrix_base<IT, D> &b) {
+template<class IT, class D1, class D2, class RT>
+void multiply_matrix(RT &result, const boolean_rowmatrix_base<IT, D1> &a, const boolean_colmatrix_base<IT, D2> &b) {
assert(a.height() == b.width());
for (unsigned r=0; r < a.height(); ++r) {
- const typename boolean_rowmatrix_base<IT, D>::row_type &row = a.get_row(r);
+ const typename boolean_rowmatrix_base<IT, D1>::row_type &row = a.get_row(r);
for (unsigned c=0; c < b.width(); ++c) {
- const typename boolean_colmatrix_base<IT, D>::column_type &col = b.get_column(c);
+ const typename boolean_colmatrix_base<IT, D2>::column_type &col = b.get_column(c);
if (count_set_intersection(row.begin(), row.end(), col.begin(), col.end()) % 2 == 1)
result.set(r, c, true);
}
template<class IT>
std::ostream& operator<<(std::ostream &os, const boolean_colrowmatrix<IT> &mat) {
+
+ os << " ";
+ for (unsigned c=0; c < mat.size(); ++c) {
+ const unsigned d = (c % 10);
+ if (d > 0)
+ os << d;
+ else
+ os << " ";
+ }
+ os << std::endl;
+
for (unsigned r=0; r < mat.size(); ++r) {
+ const unsigned d = (r % 10);
+ if (d > 0)
+ os << d;
+ else
+ os << " ";
+
for (unsigned c=0; c < mat.size(); ++c)
os << (mat.get(r,c) ? "X" : ".");
+++ /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 reduce_boundary_matrix(boolean_colrowmatrix<IT> &b, boolean_colmatrix_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 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));
- }
- }
-}
-
-}
-
-#endif
--- /dev/null
+#ifndef persistence_h_aPeinoXeiFeethuz
+#define persistence_h_aPeinoXeiFeethuz
+
+#include <iostream>
+#include <ostream>
+#include <vector>
+
+#include "simplicialcomplex.h"
+#include "booleanmatrix.h"
+
+
+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>
+class persistence {
+
+ public:
+ typedef typename simplicial_complex<MAXDIM, IT, VT>::simplex_order simplex_order;
+ typedef typename simplex_order::boundary_matrix boundary_matrix;
+ typedef boolean_colmatrix<IT> transformation_matrix;
+
+ /** A point in a diagram comprises birth- and death-index. These are the
+ * indices of the simplices of 'order' that give birth and death to the
+ * class. */
+ struct diagram_point {
+ IT birth;
+ IT death;
+ };
+
+
+ public:
+ /** Create a new peristence object on the given simplex_order */
+ persistence(const simplex_order &order) :
+ order(order),
+ bm(0),
+ rbm(0),
+ tm(0)
+ {
+ reset();
+ }
+
+ /** Reset all results gained from 'order'. */
+ void reset() {
+ done_matrices = false;
+ }
+
+ /** Get boundary matrix 'bm' of 'order', compute reduces boundary
+ * matrix 'rbm', and the transformation matrix 'tm' such that rbm = bm *
+ * tm. */
+ void compute_matrices() {
+ if (done_matrices)
+ return;
+ done_matrices = true;
+
+ bm = order.get_boundary_matrix();
+ rbm = bm;
+ tm = create_unit_matrix<transformation_matrix>(bm.size());
+ reduce_boundary_matrix(rbm, tm);
+
+ assert(rbm == bm * tm);
+ }
+
+ /** Get the boundary matrix of 'order'. */
+ const boundary_matrix& get_boundary_matrix() {
+ compute_matrices();
+ return bm;
+ }
+
+ /** Get the reduced boundary matrix of 'order'. */
+ const boundary_matrix& get_reduced_boundary_matrix() {
+ compute_matrices();
+ return rbm;
+ }
+
+ /** Get the transformation matrix of 'order', which transforms the
+ * boundary matrix to the reduced boundary matrix, when multiplied from
+ * the right. */
+ const transformation_matrix& get_transformation_matrix() {
+ compute_matrices();
+ return tm;
+ }
+
+
+ /** Print the sequence of births and deaths of cycles (homology classes). */
+ void interpret_reduction(std::ostream &os) {
+ compute_matrices();
+ for (unsigned c=0; c < bm.size(); ++c) {
+ os << c << ". inserting ";
+
+ switch (bm.get_column(c).size()) {
+ case 0:
+ os << "dummy vertex of dim -1";
+ break;
+ case 1:
+ os << "vertex";
+ break;
+ case 2:
+ os << "edge";
+ break;
+ case 3:
+ os << "triangle";
+ break;
+ case 4:
+ os << "tetrahedron";
+ break;
+ default:
+ os << "simplex";
+ break;
+ }
+
+ os << std::endl;
+ os << " boundary: " << bm.get_column(c) << std::endl;
+
+ typename boolean_colrowmatrix<IT>::colbase::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 {
+ os << " \033[1;31mdeath\033[0;m of a cycle: " << rbm.get_column(c) << std::endl;
+ os << " boundary of: " << tm.get_column(c) << std::endl;
+ }
+ os << std::endl;
+ }
+ }
+
+ private:
+
+ /** The underlying simplex order of this diagram. */
+ const simplex_order ℴ
+
+ bool done_matrices;
+
+ boundary_matrix bm;
+ boundary_matrix rbm;
+ transformation_matrix tm;
+};
+
+
+/** 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 reduce_boundary_matrix(boolean_colrowmatrix<IT> &b, boolean_colmatrix_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 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));
+ }
+ }
+}
+
+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;
+}
+
+
+} // namespace libstick
+
+#endif
simplices.push_back(s);
}
+ /** Add an array of simplices */
+ void add_simplices(Simplex* sarray, size_t count) {
+ for (unsigned i=0; i < count; ++i)
+ add_simplex(sarray[i]);
+ }
+
/** 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 is_complex() const {
#include "booleanmatrix.h"
#include "simplicialcomplex.h"
+#include "persistence.h"
using namespace std;
ts.add(auto_ptr<Test::Suite>(new boolean_matrix_TestSuite));
ts.add(auto_ptr<Test::Suite>(new simplicial_complex_TestSuite));
+ ts.add(auto_ptr<Test::Suite>(new persistence_TestSuite));
Test::TextOutput output(Test::TextOutput::Verbose);
return ts.run(output);
--- /dev/null
+#ifndef persistence_h_Fahlaewahgaebaqu
+#define persistence_h_Fahlaewahgaebaqu
+
+#include <cpptest.h>
+#include <cpptest-suite.h>
+
+#include <libstick-0.1/persistence.h>
+
+using namespace libstick;
+
+
+class persistence_TestSuite: public Test::Suite {
+
+ private:
+ typedef simplicial_complex<3, uint32_t, double> scomplex;
+ typedef persistence<3, uint32_t, double> pers;
+ typedef pers::boundary_matrix bm;
+ typedef pers::transformation_matrix tm;
+
+ bool setupcalled;
+ scomplex ball, ring, torus;
+ scomplex::simplex_order oball, oring, otorus;
+
+ public:
+ persistence_TestSuite() :
+ setupcalled(false),
+ oball(ball),
+ oring(ring),
+ otorus(torus)
+ {
+ TEST_ADD(persistence_TestSuite::test_matrix_reduction);
+ }
+
+ protected:
+ virtual void setup() {
+ if (setupcalled)
+ return;
+ setupcalled = true;
+
+ // 7 ------------------------------------- 8
+ // \ /
+ // \ 4 ------------------------- 5 /
+ // \ \ / /
+ // \ \ 1 ------- 2 / /
+ // \ \ \ / / /
+ // \ \ \ / / /
+ // \ \ \ / / /
+ // \ \ \ / / /
+ // \ \ 3 / /
+ // \ \ / /
+ // \ \ / /
+ // \ \ / /
+ // \ \ / /
+ // \ \ / /
+ // \ \ / /
+ // \ 6 /
+ // \ /
+ // \ /
+ // \ /
+ // 9
+
+ scomplex::Simplex ssring[] = {
+ // dimension, faces, value...
+ /* 1 */ {0, {}, 1},
+ /* 2 */ {0, {}, 2},
+ /* 3 */ {0, {}, 3},
+ /* 4 */ {0, {}, 4},
+ /* 5 */ {0, {}, 5},
+ /* 6 */ {0, {}, 6},
+ /* 7 */ {1, {2, 3}, 6.01},
+ /* 8 */ {1, {3, 1}, 6.02},
+ /* 9 */ {1, {1, 2}, 6.03},
+ /* 10 */ {1, {4, 5}, 6.04},
+ /* 11 */ {1, {5, 6}, 6.05},
+ /* 12 */ {1, {6, 4}, 6.06},
+ /* 13 */ {1, {1, 4}, 6.07},
+ /* 14 */ {1, {1, 5}, 6.08},
+ /* 15 */ {2, {13, 14, 10}, 6.0801},
+ /* 16 */ {1, {3, 6}, 6.09},
+ /* 17 */ {1, {2, 5}, 6.10},
+ /* 18 */ {1, {2, 6}, 6.11},
+ /* 19 */ {2, {9, 14, 17}, 6.1101},
+ /* 20 */ {2, {7, 16, 18}, 6.1102},
+ /* 21 */ {2, {11, 17, 18}, 6.1103},
+ /* 22 */ {1, {1, 6}, 6.12},
+ /* 23 */ {2, {12, 13, 22}, 6.1201},
+ /* 24 */ {2, {8, 16, 22}, 6.1202}
+ };
+ const size_t cntssring = sizeof(ssring)/sizeof(scomplex::Simplex);
+ ring.add_simplices(ssring, cntssring);
+ oring.reset();
+
+ scomplex::Simplex sstorus[] = {
+ // dimension, faces, value...
+ /* 25 */ {0, {}, 7},
+ /* 26 */ {0, {}, 8},
+ /* 27 */ {0, {}, 9},
+ /* 28 */ {1, {25, 26}, 9.01},
+ /* 29 */ {1, {26, 27}, 9.02},
+ /* 30 */ {1, {25, 27}, 9.03},
+ /* 31 */ {1, {25, 1}, 9.04},
+ /* 32 */ {1, {26, 1}, 9.05},
+ /* 33 */ {2, {28, 31, 32}, 9.0501},
+ /* 34 */ {1, {27, 1}, 9.06},
+ /* 35 */ {2, {30, 31, 34}, 9.0601},
+ /* 36 */ {1, {27, 3}, 9.07},
+ /* 37 */ {2, {36, 34, 8}, 9.0701},
+ /* 38 */ {1, {27, 2}, 9.08},
+ /* 39 */ {1, {26, 2}, 9.09},
+ /* 40 */ {2, {9, 32, 39}, 9.0901},
+ /* 41 */ {2, {29, 38, 39}, 9.0902},
+ /* 42 */ {2, {7, 38, 36}, 9.0903},
+ /* 43 */ {1, {4, 25}, 9.10},
+ /* 44 */ {1, {5, 26}, 9.11},
+ /* 45 */ {1, {6, 27}, 9.12},
+ /* 46 */ {1, {4, 26}, 9.13},
+ /* 47 */ {1, {4, 27}, 9.14},
+ /* 48 */ {1, {5, 27}, 9.15},
+ /* 49 */ {2, {43, 47, 30}, 9.1501},
+ /* 50 */ {2, {12, 45, 47}, 9.1502},
+ /* 51 */ {2, {29, 44, 48}, 9.1503},
+ /* 52 */ {2, {48, 11, 45}, 9.1504},
+ /* 53 */ {2, {10, 46, 44}, 9.1505},
+ /* 54 */ {2, {43, 46, 28}, 9.1506},
+ };
+ const size_t cntsstorus = sizeof(sstorus)/sizeof(scomplex::Simplex);
+ torus = ring;
+ torus.add_simplices(sstorus, cntsstorus);
+ otorus.reset();
+
+ scomplex::Simplex ssball[] = {
+ // dimension, faces, value...
+ {0, {}, 1},
+ {0, {}, 2},
+ {0, {}, 3},
+ {0, {}, 4},
+ {0, {}, 5},
+ {1, {1, 2}, 6},
+ {1, {2, 3}, 7},
+ {1, {3, 4}, 8},
+ {1, {4, 1}, 9},
+ {1, {1, 5}, 10},
+ {1, {2, 5}, 11},
+ {1, {3, 5}, 12},
+ {1, {4, 5}, 13},
+ {2, {6, 10, 11}, 14},
+ {2, {7, 11, 12}, 15},
+ {2, {8, 12, 13}, 16},
+ {2, {9, 13, 10}, 17},
+ {1, {1, 3}, 18},
+ {2, {18, 6, 7}, 19},
+ {2, {18, 8, 9}, 20},
+ {2, {18, 10, 12}, 21},
+ {3, {21, 14, 15, 19}, 22},
+ {3, {21, 16, 17, 20}, 23},
+ };
+ const size_t cntssball = sizeof(ssball)/sizeof(scomplex::Simplex);
+ ball.add_simplices(ssball, cntssball);
+ oball.reset();
+ }
+
+ virtual void tear_down() {
+ }
+
+ void test_matrix_reduction() {
+ pers ringp(oring);
+ 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);
+
+ pers torusp(otorus);
+ 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);
+
+
+ pers ballp(oball);
+ 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);
+
+ //std::cout << std::endl;
+ //ballp.interpret_reduction(std::cout);
+ }
+};
+
+#endif
#include <cpptest-suite.h>
#include <libstick-0.1/simplicialcomplex.h>
-#include <libstick-0.1/matrixreduction.h>
using namespace libstick;
class simplicial_complex_TestSuite: public Test::Suite {
private:
- typedef simplicial_complex<2, uint32_t, double> scomplex;
+ typedef simplicial_complex<3, uint32_t, double> scomplex;
typedef scomplex::simplex_order::boundary_matrix bm;
bool setupcalled;
TEST_ADD(simplicial_complex_TestSuite::test_is_order_filtration);
TEST_ADD(simplicial_complex_TestSuite::test_is_order_monotone);
TEST_ADD(simplicial_complex_TestSuite::test_boundary_matrix);
- TEST_ADD(simplicial_complex_TestSuite::test_matrix_reduction);
}
protected:
return;
setupcalled = true;
- const unsigned num = 11;
- scomplex::Simplex ss[num] = {
+ scomplex::Simplex ss[] = {
// dimension, faces, value...
{0, {}, 1},
{0, {}, 2},
{2, {7, 8, 9}, 10},
{2, {5, 6, 9}, 11}
};
+ const size_t cntss = sizeof(ss)/sizeof(scomplex::Simplex);
+ c1.add_simplices(ss, cntss);
// This is o1 This is o2 This is o3(b)
// (value = index) (values) (values)
// | \| | \| | \|
// 4 ----7---- 3 4 ----7---- 3 4 ----7---- 3
- // Build the complex
- for (unsigned i=0; i<num; ++i)
- c1.add_simplex(ss[i]);
-
c2 = c1;
c2.simplices[6].value = 12;
void test_boundary_matrix() {
bm mat1 = o1.get_boundary_matrix();
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, 5, true);
- mat1e.set(2, 5, true);
- mat1e.set(2, 6, true);
- mat1e.set(3, 6, true);
- mat1e.set(3, 7, 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);
+ uint32_t mat1e_coords[][2] = {
+ {0, 1}, {0, 2}, {0, 3}, {0, 4}, {1, 5}, {2, 5}, {2, 6}, {3, 6}, {3, 7}, {4, 7}, {1, 8}, {4, 8},
+ {1, 9}, {3, 9}, {7, 10}, {8, 10}, {9, 10}, {5, 11}, {6, 11}, {9, 11}
+ };
+ mat1e.set_all(mat1e_coords, sizeof(mat1e_coords)/(2*sizeof(uint32_t)), true);
TEST_ASSERT(mat1 == mat1e);
+
bm mat3b = o3b.get_boundary_matrix();
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, 5, true);
- mat3be.set(2, 5, true);
- mat3be.set(3, 6, 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);
+ uint32_t mat3be_coords[][2] = {
+ {0, 1}, {0, 2}, {0, 3}, {0, 4}, {1, 5}, {2, 5}, {3, 6}, {4, 6}, {1, 7}, {4, 7}, {1, 8}, {3, 8},
+ {6, 9}, {7, 9}, {8, 9}, {2, 10}, {3, 10}, {5, 11}, {8, 11}, {10, 11}
+ };
+ mat3be.set_all(mat3be_coords, sizeof(mat3be_coords)/(2*sizeof(uint32_t)), true);
TEST_ASSERT(mat3b == mat3be);
//std::cout << mat3b << std::endl << std::endl;
//std::cout << ((bm::colbase) mat3b) << std::endl << std::endl;
//std::cout << ((bm::rowbase) mat3b) << std::endl << std::endl;
}
-
- void test_matrix_reduction() {
- bm b = o3b.get_boundary_matrix();
- bm v = create_unit_matrix<bm>(b.size());
-
- bm b_orig = b;
- reduce_boundary_matrix(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;
- }
};
#endif
projects / libstick.git / commitdiff