From: Stefan Huber <shuber@sthu.org>
Date: Thu, 7 Nov 2013 18:37:59 +0000 (+0100)
Subject: Add booleanmatrix and simplicialcomplex
X-Git-Tag: v0.1~29
X-Git-Url: https://git.sthu.org/?a=commitdiff_plain;h=267194caeb2fe56530b11c440fefd01c344e5482;p=libstick.git

Add booleanmatrix and simplicialcomplex
---

diff --git a/CMakeLists.txt b/CMakeLists.txt
index c9d4fb2..b5116d9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,30 +1,42 @@
 cmake_minimum_required (VERSION 2.6)
-project (stick) 
+project (libstick)
 
 
 if(NOT CMAKE_SYSTEM_NAME STREQUAL "Windows")
-	set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -ansi")
+	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -ansi -pedantic")
+	# I would like to use -pedantic, but VTK devs are not sufficiently
+	# pedantic
+	set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Werror")
 	set(CMAKE_EXE_LINKER_FLAGS "-Wl,--as-needed" )
 endif()
 
-set(CMAKE_C_FLAGS_RELEASE "-DNDEBUG")
+set(CMAKE_CXX_FLAGS_RELEASE "-DNDEBUG")
 if(CMAKE_BUILD_TYPE STREQUAL "Release")
-	# CMAKE_C_FLAGS_RELEASE is appended, but we need to prepend -O2 to take
+	# CMAKE_CXX_FLAGS_RELEASE is appended, but we need to prepend -O2 to take
 	# effect
-	set(CMAKE_C_FLAGS "-O2 ${CMAKE_C_FLAGS}")
+	set(CMAKE_CXX_FLAGS "-O2 ${CMAKE_CXX_FLAGS}")
 endif()
 
 set(PACKAGE_STRING stick)
 set(PACKAGE_BUGREPORT stefan.huber@ist.ac.at)
 set(PACKAGE_VERSION 0.1)
 
-if(NOT CMAKE_SYSTEM_NAME STREQUAL "Windows")
-	set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Wall -Wextra -std=c99")
-	# I would like to use -pedantic, but VTK devs are not sufficiently
-	# pedantic
-	set(CMAKE_C_FLAGS_DEBUG "${CMAKE_C_FLAGS_DEBUG} -Werror")
-endif()
+option(WITH_UNITTESTS "Enable unit tests" "OFF")
 
 add_subdirectory(include)
 add_subdirectory(lib)
 add_subdirectory(src)
+
+if(WITH_UNITTESTS)
+	add_subdirectory(tests)
+endif()
+
+
+########### Add uninstall target ###############
+configure_file(
+    "${CMAKE_CURRENT_SOURCE_DIR}/cmake/cmake_uninstall.cmake.in"
+    "${CMAKE_CURRENT_BINARY_DIR}/cmake_uninstall.cmake"
+    IMMEDIATE @ONLY)
+
+add_custom_target(uninstall
+    COMMAND ${CMAKE_COMMAND} -P ${CMAKE_CURRENT_BINARY_DIR}/cmake_uninstall.cmake)
diff --git a/cmake/cmake_uninstall.cmake.in b/cmake/cmake_uninstall.cmake.in
new file mode 100644
index 0000000..c6d8094
--- /dev/null
+++ b/cmake/cmake_uninstall.cmake.in
@@ -0,0 +1,22 @@
+if (NOT EXISTS "@CMAKE_CURRENT_BINARY_DIR@/install_manifest.txt")
+    message(FATAL_ERROR "Cannot find install manifest: \"@CMAKE_CURRENT_BINARY_DIR@/install_manifest.txt\"")
+endif(NOT EXISTS "@CMAKE_CURRENT_BINARY_DIR@/install_manifest.txt")
+
+file(READ "@CMAKE_CURRENT_BINARY_DIR@/install_manifest.txt" files)
+string(REGEX REPLACE "\n" ";" files "${files}")
+list(REVERSE files)
+foreach (file ${files})
+    message(STATUS "Uninstalling \"$ENV{DESTDIR}${file}\"")
+    if (EXISTS "$ENV{DESTDIR}${file}")
+        execute_process(
+            COMMAND @CMAKE_COMMAND@ -E remove "$ENV{DESTDIR}${file}"
+            OUTPUT_VARIABLE rm_out
+            RESULT_VARIABLE rm_retval
+        )
+        if(NOT ${rm_retval} EQUAL 0)
+            message(FATAL_ERROR "Problem when removing \"$ENV{DESTDIR}${file}\"")
+        endif (NOT ${rm_retval} EQUAL 0)
+    else (EXISTS "$ENV{DESTDIR}${file}")
+        message(STATUS "File \"$ENV{DESTDIR}${file}\" does not exist.")
+    endif (EXISTS "$ENV{DESTDIR}${file}")
+endforeach(file)
diff --git a/include/CMakeLists.txt b/include/CMakeLists.txt
index 5e88e3f..caf6524 100644
--- a/include/CMakeLists.txt
+++ b/include/CMakeLists.txt
@@ -1,2 +1 @@
-set(INCLUDE_DIR include/${PACKAGE_STRING}-${PACKAGE_VERSION})
-install (DIRECTORY ./ DESTINATION ${INCLUDE_DIR} FILES_MATCHING PATTERN "*.h")
+install (DIRECTORY ./ DESTINATION include FILES_MATCHING PATTERN "*.h")
diff --git a/include/libstick-0.1/booleanmatrix.h b/include/libstick-0.1/booleanmatrix.h
new file mode 100644
index 0000000..ab368d2
--- /dev/null
+++ b/include/libstick-0.1/booleanmatrix.h
@@ -0,0 +1,417 @@
+#ifndef booleanmatrix_h_AechohNgakoghahV
+#define booleanmatrix_h_AechohNgakoghahV
+
+#include <stdint.h>
+#include <cassert>
+#include <cstdlib>
+#include <algorithm>
+#include <vector>
+#include <list>
+#include <set>
+#include <ostream>
+#include <iterator>
+
+
+namespace libstick {
+
+
+/** The base class of BooleanColMatrix and BooleanColRowMatrix which implements
+ * the common logic of both. */
+template<class IT, class D>
+class BooleanColMatrix_base {
+
+    public:
+        typedef IT index_type;
+        typedef std::vector<index_type> column_type;
+        typedef D derived;
+
+        /** Create a matrix with 'width' columns, initalized with zero entries. */
+        BooleanColMatrix_base(size_t width) :
+            cols(width) {
+        }
+
+        /** Get height resp. width of the matrix. */
+        size_t width() const {
+            return cols.size();
+        }
+
+        /** Get the matrix entry at row 'r' and column 'c'. */
+        bool get(index_type r, index_type c) const {
+            assert(c < width());
+            const column_type &col = getColumn(c);
+            return binary_search(col.begin(), col.end(), r);
+        }
+
+        /** Set the matrix entry at row 'r' and column 'c'. */
+        void set(index_type r, index_type c, bool value) {
+            getDerived()->_set(r, c, value);
+        }
+
+        /** Get the c-th column. */
+        const column_type& getColumn(index_type c) const {
+            assert(c < width());
+            return cols[c];
+        }
+
+        /** 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) {
+            assert(c < width());
+
+            // Flip all entries that are set in 'col'.
+            for (typename column_type::const_iterator it = col.begin(); it != col.end(); ++it)
+                set(*it, c, !get(*it, c));
+        }
+
+        /** Two matrices are equal iff they have the same entries */
+        bool operator==(const BooleanColMatrix_base<IT, D> &m) const {
+            return cols == m.cols;
+        }
+
+        /** Two matrices are equal iff they have the same entries */
+        bool operator!=(const BooleanColMatrix_base<IT, D> &m) const {
+            return !(*this == m);
+        }
+
+    protected:
+
+
+        /** Set the matrix entry at row 'r' and column 'c'. */
+        void _set(index_type r, index_type c, bool value) {
+            assert(c < width());
+
+            column_type &col = cols.at(c);
+            // Let us see where to insert the new element
+            typename column_type::iterator it = lower_bound(col.begin(), col.end(), r);
+            bool exists = (it != col.end() && *it == r);
+            assert(get(r,c) == exists);
+
+            // Add 'r' to c-th column
+            if (value) {
+                // r is new, insert it
+                if (!exists)
+                    col.insert(it, r);
+                assert(get(r,c));
+            }
+            // Remove the element
+            else {
+                if (exists)
+                    col.erase(it);
+                assert(!get(r,c));
+            }
+
+#ifndef NDEBUG
+            for (unsigned i=1; i < col.size(); i++)
+                assert(col[i-1] < col[i]);
+#endif
+        }
+
+    private:
+        /** The matrix is the set of columns. */
+        std::vector<column_type> cols;
+
+        /** A cast to the derived type */
+        derived* getDerived() {
+            return static_cast<derived*>(this);
+        }
+};
+
+
+/** This is boolean matrix that is a std::vector of column-vectors and in each
+ * column-vector we save the row-indices of the 1s. It is designed to handle a
+ * few 1s and column-wise operations well. */
+template<class IT=uint32_t>
+class BooleanColMatrix : public BooleanColMatrix_base<IT, BooleanColMatrix<IT> > {
+
+    public:
+        typedef IT index_type;
+        typedef BooleanColMatrix_base<IT, BooleanColMatrix<IT> > base;
+
+        /** Create a matrix with 'width' columns, initalized with zero entries. */
+        BooleanColMatrix(size_t columns) :
+            base(columns) {
+        }
+
+        /** Set the matrix entry at row 'r' and column 'c'. */
+        void _set(index_type r, index_type c, bool value) {
+            base::_set(r, c, value);
+        }
+};
+
+
+/** The base class of BooleanRowMatrix and BooleanColRowMatrix which implements
+ * the common logic of both. */
+template<class IT, class D>
+class BooleanRowMatrix_base {
+
+    public:
+        typedef IT index_type;
+        typedef std::vector<index_type> row_type;
+        typedef D derived;
+
+        /** Create a matrix with 'height' rows, initalized with zero entries.
+         * */
+        BooleanRowMatrix_base(size_t height) :
+            rows(height) {
+        }
+
+        /** Get height resp. width of the matrix. */
+        size_t height() const {
+            return rows.size();
+        }
+
+        /** Get the matrix entry at row 'r' and column 'c'. */
+        bool get(index_type r, index_type c) const {
+            assert(r < height());
+            const row_type &row = getRow(r);
+            return binary_search(row.begin(), row.end(), c);
+        }
+
+        /** Set the matrix entry at row 'r' and column 'c'. */
+        void set(index_type r, index_type c, bool value) {
+            getDerived()->_set(r, c, value);
+        }
+
+        /** Get the r-th row. */
+        const row_type& getRow(index_type r) const {
+            assert(r < height());
+            return rows[r];
+        }
+
+        /** 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) {
+            assert(r < height());
+
+            // Flip all entries that are set in 'row'.
+            for (typename row_type::const_iterator it = row.begin(); it != row.end(); ++it)
+                set(r, *it, !get(r, *it));
+        }
+
+        /** Two matrices are equal iff they have the same entries */
+        bool operator==(const BooleanRowMatrix_base<IT, D> &m) const {
+            return rows == m.rows;
+        }
+
+        /** Two matrices are equal iff they have the same entries */
+        bool operator!=(const BooleanRowMatrix_base<IT, D> &m) const {
+            return !(*this == m);
+        }
+
+    protected:
+        /** Set the matrix entry at row 'r' and column 'c'. */
+        void _set(index_type r, index_type c, bool value) {
+            assert(r < height());
+
+            row_type &row = rows.at(r);
+            // Let us see where to insert/remove the new element
+            typename row_type::iterator it = lower_bound(row.begin(), row.end(), c);
+            bool exists = (it != row.end() && *it == c);
+            assert(get(r,c) == exists);
+
+            // Add 'r' to c-th column
+            if (value) {
+                // r is new, insert it
+                if (!exists)
+                    row.insert(it, c);
+                assert(get(r,c));
+            }
+            // Remove the element
+            else {
+                if (exists)
+                    row.erase(it);
+                assert(!get(r,c));
+            }
+
+#ifndef NDEBUG
+            for (unsigned i=1; i < row.size(); i++)
+                assert(row[i-1] < row[i]);
+#endif
+        }
+
+    private:
+        derived* getDerived() {
+            return static_cast<derived*>(this);
+        }
+
+        /** The matrix is the set of columns. */
+        std::vector<row_type> rows;
+};
+
+
+/** This is boolean matrix that is a std::vector of row-vectors and in each
+ * row-vector we save the column-indices of the 1s. It is designed to handle a
+ * few 1s and row-wise operations well. */
+template<class IT=uint32_t>
+class BooleanRowMatrix : public BooleanRowMatrix_base<IT, BooleanRowMatrix<IT> > {
+
+    public:
+        typedef IT index_type;
+        typedef BooleanRowMatrix_base<IT, BooleanRowMatrix<IT> > base;
+
+        /** Create a matrix with 'height' rows, initalized with zero entries. */
+        BooleanRowMatrix(size_t height) :
+            base(height) {
+        }
+
+        /** Set the matrix entry at row 'r' and column 'c'. */
+        void _set(index_type r, index_type c, bool value) {
+            base::_set(r, c, value);
+        }
+};
+
+
+/** This is a boolean matrix that supports. It is designed to handle a
+ * few 1s and row- and column-wise operations well. */
+template<class IT=uint32_t>
+class BooleanColRowMatrix : public BooleanColMatrix_base<IT, BooleanColRowMatrix<IT> >,
+                            public BooleanRowMatrix_base<IT, BooleanColRowMatrix<IT> > {
+
+    public:
+        typedef BooleanColMatrix_base<IT, BooleanColRowMatrix<IT> > colbase;
+        typedef BooleanRowMatrix_base<IT, BooleanColRowMatrix<IT> > rowbase;
+
+    public:
+        typedef IT index_type;
+
+        /** Create a size x size matrix, initialized with zeros. */
+        BooleanColRowMatrix(size_t size) :
+            colbase(size),
+            rowbase(size) {
+        }
+
+        /** Override implementation. */
+        void _set(index_type r, index_type c, bool value) {
+            colbase::_set(r, c, value);
+            rowbase::_set(r, c, value);
+        }
+
+    public:
+        /** Get height resp. width of the matrix. */
+        size_t size() const {
+            assert(colbase::width() == rowbase::height());
+            return colbase::width();
+        }
+
+        /** Set the matrix entry at row 'r' and column 'c'. */
+        void set(index_type r, index_type c, bool value) {
+            //Calling set() for one base suffices, static polymorphism
+            //calls _set() anyhow.
+            //rowbase::set(r, c, value);
+            colbase::set(r, c, 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));
+            return colbase::get(r, c);
+        }
+
+        /** Two matrices are equal iff they have the same entries */
+        bool operator==(const BooleanColRowMatrix<IT> &m) const {
+            assert(rowbase::operator==(m) == colbase::operator==(m));
+            return colbase::operator==(m);
+        }
+
+        /** Two matrices are equal iff they have the same entries */
+        bool operator!=(const BooleanColRowMatrix<IT> &m) const {
+            return !(*this == m);
+        }
+};
+
+template<class IT>
+std::ostream& operator<<(std::ostream &os, const BooleanColRowMatrix<IT> &mat) {
+    for (unsigned r=0; r < mat.size(); ++r) {
+        for (unsigned c=0; c < mat.size(); ++c)
+            os << (mat.get(r,c) ? "X" : ".");
+
+        if (r < mat.size() - 1)
+            os << std::endl;
+    }
+    return os;
+}
+
+template<class IT>
+std::ostream& operator<<(std::ostream &os, BooleanColRowMatrix<IT> &mat) {
+    const BooleanColRowMatrix<IT> &m = mat;
+    return os << m;
+}
+
+template<class IT, class D>
+std::ostream& operator<<(std::ostream &os, const BooleanRowMatrix_base<IT, D> &mat) {
+    for (unsigned r=0; r < mat.height(); ++r) {
+        // Get the sorted r-th row
+        std::list<IT> row(mat.getRow(r).size());
+        copy(mat.getRow(r).begin(), mat.getRow(r).end(), row.begin());
+
+        os << "|";
+        // Print the elements, and remove them
+        for (unsigned c=0; row.size() > 0; ++c) {
+            if (row.front() == c) {
+                os << "X";
+                row.pop_front();
+            } else
+                os << ".";
+        }
+
+        if (r < mat.height() - 1)
+            os << std::endl;
+    }
+    return os;
+}
+
+template<class IT, class D>
+std::ostream& operator<<(std::ostream &os, BooleanRowMatrix_base<IT, D> &mat) {
+    const BooleanRowMatrix_base<IT, D> &m = mat;
+    return os << m;
+}
+
+template<class IT, class D>
+std::ostream& operator<<(std::ostream &os, const BooleanColMatrix_base<IT, D> &mat) {
+    // The sorted columns
+    std::vector<std::list<IT> > cols;
+    // The max height (max. value) of all columns
+    IT height = 0;
+
+    for (unsigned c=0; c < mat.width(); ++c) {
+        // Get the sorted c-th column
+        std::list<IT> col(mat.getColumn(c).size());
+        copy(mat.getColumn(c).begin(), mat.getColumn(c).end(), col.begin());
+        cols.push_back(col);
+
+        os << "-";
+
+        if (col.size() != 0)
+            height = std::max(height, col.back());
+    }
+
+    os << std::endl;
+
+    for (unsigned r=0; r <= height; ++r) {
+        // Print the elements, and remove them
+        for (unsigned c=0; c < mat.width(); ++c) {
+            std::list<IT> &col = cols.at(c);
+            // This column is already empty
+            if (col.size() == 0)
+                os << " ";
+            else if (col.front() == r) {
+                os << "X";
+                col.pop_front();
+            } else
+                os << ".";
+        }
+
+        os << std::endl;
+    }
+    return os;
+}
+
+template<class IT, class D>
+std::ostream& operator<<(std::ostream &os, BooleanColMatrix_base<IT, D> &mat) {
+    const BooleanColMatrix_base<IT, D> &m = mat;
+    return os << m;
+}
+
+}
+
+#endif
diff --git a/include/libstick-0.1/simplicialcomplex.h b/include/libstick-0.1/simplicialcomplex.h
new file mode 100644
index 0000000..d829721
--- /dev/null
+++ b/include/libstick-0.1/simplicialcomplex.h
@@ -0,0 +1,266 @@
+#ifndef simplicialcomplex_h_nealaezeojeeChuh
+#define simplicialcomplex_h_nealaezeojeeChuh
+
+#include <stdint.h>
+#include <cstdlib>
+#include <cstring>
+#include <algorithm>
+#include <vector>
+
+#include <iostream>
+
+#include <libstick-0.1/booleanmatrix.h>
+
+
+namespace libstick {
+
+/** 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>
+class SimplicialComplex {
+
+    public:
+        /** The type of this class. */
+        typedef SimplicialComplex<MAXDIM, IT, VT> simplcompltype;
+        /** Type of indices of simplices. */
+        typedef IT index_type;
+        /** To every simplex a function value is assigned according to which a
+         * filtration is considered. This is the value type of the function. */
+        typedef VT value_type;
+
+        /** A simplex of the complex. */
+        struct Simplex {
+            /** Dimension of the simplex. */
+            size_t 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);
+
+                Simplex s;
+                s.dim = dim;
+                s.value = value;
+                memcpy(s.faces, faces, face_count_bydim(dim)*sizeof(index_type));
+
+                return s;
+            }
+
+            /** Get number of faces. */
+            size_t face_count() const {
+                return face_count_bydim(dim);
+            }
+
+            /** Get number of faces of a dim-dimensional simplex. */
+            static size_t face_count_bydim(size_t dim) {
+                if (dim == 0)
+                    return 0;
+                return dim + 1;
+            }
+        };
+
+        /** An order of the simplices of complex c. An order can be interpreted
+         * as a permuation of the complex's std::vector of simplices.  */
+        class SimplexOrder {
+
+            public:
+                typedef BooleanColRowMatrix<IT> BoundaryMatrix;
+
+                /** Create a standard order of the complex c, i.e., the identity permutation. */
+                SimplexOrder(const simplcompltype &c) :
+                    c(c)
+                    {
+                    reset();
+                }
+
+                /** Reset order to the identity permutation of the complex's simplices. */
+                void reset() {
+                    order.clear();
+                    for (unsigned i=0; i < c.size(); ++i)
+                        order.push_back(i);
+                    revorder = order;
+                }
+
+                /** Return number of simplices. */
+                size_t size() const {
+                    assert(order.size() == revorder.size());
+                    return order.size();
+                }
+
+                /** Get i-th simplex in the simplex order. */
+                const Simplex& getSimplex(size_t i) const {
+                    assert(i < size());
+                    return c.simplices[order.at(i)];
+                }
+
+                /** Returns true iff the faces of simplex i are before i in this order. */
+                bool isFiltration() const {
+                    assert(size() == c.size());
+
+                    for (unsigned i=0; i < size(); ++i)
+                        for (unsigned f=0; f < getSimplex(i).face_count(); ++f)
+                            if (revorder[getSimplex(i).faces[f]] >= i)
+                                return false;
+
+                    return true;
+                }
+
+                /** Returns true iff isFiltration() gives true and values of simplices
+                 * are monotone w.r.t. this order of simplices. */
+                bool isMonotone() const {
+                    assert(size() == c.size());
+
+                    for (unsigned i=1; i < size(); ++i)
+                        if (getSimplex(i-1).value > getSimplex(i).value)
+                            return false;
+
+                    return isFiltration();
+                }
+
+                /** Sort simplices such that isMonotone() gives true. This
+                 * requires that the complex's isMonotone() gave true
+                 * beforehand.*/
+                void makeMonotoneFiltration() {
+                    assert(c.isMonotone());
+
+                    sort(order.begin(), order.end(), cmpMonotoneFiltration(c));
+                    restoreRevorderFromOrder();
+
+                    assert(c.isMonotone());
+                    assert(isFiltration());
+                    assert(isMonotone());
+                }
+
+                /** Get the boundary matrix of the complex according to this order. */
+                BoundaryMatrix getBoundaryMatrix() const {
+                    BoundaryMatrix mat(size());
+
+                    for (unsigned c=0; c < size(); ++c)
+                        for(unsigned r=0; r < getSimplex(c).face_count(); ++r)
+                            mat.set(revorder[getSimplex(c).faces[r]], c, true);
+
+                    return mat;
+                }
+
+            private:
+                /** Reconstruct 'revorder' by inverting the permutation given by 'order'. */
+                void restoreRevorderFromOrder() {
+                    // Make revorder * order the identity permutation
+                    for (unsigned i=0; i < size(); ++i)
+                        revorder[order[i]] = i;
+                }
+
+                /** The complex of which we consider a simplex order. */
+                const simplcompltype &c;
+
+                /** The i-th simplex in order is the order[i]-th simplex of the
+                 * complex. 'order' can be seen as a permutation of the
+                 * simplices saved in 'c'. */
+                std::vector<index_type> order;
+
+                /** The i-th simplex in the complex is the revorder[i]-th
+                 * simplex in order. 'revorder' can be seen as the inverse
+                 * permutation saved in 'order'. */
+                std::vector<index_type> revorder;
+        };
+
+    public:
+        /** 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) {
+            addSimplex(Simplex::create(dim, faces, value));
+        }
+
+        /** 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(Simplex s) {
+            // Check requirements for faces
+            for (unsigned i=0; i < s.face_count(); ++i) {
+                // Faces are already in complex.
+                assert(s.faces[i] < size());
+                // Faces have dimension dim-1
+                assert(simplices[s.faces[i]].dim == s.dim-1);
+            }
+
+            simplices.push_back(s);
+        }
+
+        /** 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) {
+
+                    if (s.faces[f] >= size())
+                        return false;
+
+                    const Simplex &face = simplices[s.faces[f]];
+                    if (face.dim != s.dim-1)
+                        return false;
+                }
+            }
+            return true;
+        }
+
+        /** Returns true iff simplex's values are monotone w.r.t.
+         * face-inclusion, i.e., for each simplex its value is not smaller than
+         * the values of its faces. Requires that isComplex() gives true. */
+        bool isMonotone() const {
+            assert(isComplex());
+
+            typename std::vector<Simplex>::const_iterator it = ++simplices.begin();
+            for (; it != simplices.end(); ++it)
+                for (unsigned f=0; f < it->face_count(); ++f)
+                    if (simplices[it->faces[f]].value > it->value)
+                        return false;
+
+            return true;
+        }
+
+    private:
+        /** Compares (operator<) two simplices (i.e. indices) in a
+         * SimplexOrder w.r.t. lexicographical order on (value,
+         * dimension)-tuples. */
+        struct cmpMonotoneFiltration {
+            const SimplicialComplex &c;
+
+            cmpMonotoneFiltration(const SimplicialComplex &c) :
+                c(c){
+                }
+
+            bool operator()(index_type i, index_type j) {
+                const Simplex& si = c.simplices[i];
+                const Simplex& sj = c.simplices[j];
+
+                if (si.value < sj.value)
+                    return true;
+                else if (si.value == sj.value)
+                    return si.dim < sj.dim;
+                else
+                    return false;
+            }
+        };
+
+    public:
+        /** A list of simplices */
+        std::vector<Simplex> simplices;
+};
+
+}
+
+
+#endif
diff --git a/include/libstick-0.1/test.h b/include/libstick-0.1/test.h
new file mode 100644
index 0000000..e610c5a
--- /dev/null
+++ b/include/libstick-0.1/test.h
@@ -0,0 +1,6 @@
+#ifndef test_h_yiQuuiphuilutieX
+#define test_h_yiQuuiphuilutieX
+
+extern "C" void test();
+
+#endif
diff --git a/include/stick.h b/include/stick.h
deleted file mode 100644
index 451a77e..0000000
--- a/include/stick.h
+++ /dev/null
@@ -1,6 +0,0 @@
-#ifndef lyxtor_h
-#define lyxtor_h
-
-void test();
-
-#endif
diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt
index 3636f2e..2d03559 100644
--- a/lib/CMakeLists.txt
+++ b/lib/CMakeLists.txt
@@ -1,6 +1,6 @@
-set(stick_SRC stick.c)
+set(stick_SRC test.cc)
 
-include_directories(${stick_SOURCE_DIR}/include)
+include_directories(${libstick_SOURCE_DIR}/include)
 
 configure_file(${CMAKE_CURRENT_SOURCE_DIR}/libstick.pc.cmake ${CMAKE_CURRENT_BINARY_DIR}/libstick.pc @ONLY)
 
diff --git a/lib/stick.c b/lib/stick.c
deleted file mode 100644
index d6adeff..0000000
--- a/lib/stick.c
+++ /dev/null
@@ -1,8 +0,0 @@
-#include <stdio.h>
-
-#include "stick.h"
-
-
-void test() {
-	printf("Hello world.\n");
-}
diff --git a/lib/test.cc b/lib/test.cc
new file mode 100644
index 0000000..75a9499
--- /dev/null
+++ b/lib/test.cc
@@ -0,0 +1,8 @@
+#include <stdio.h>
+
+#include "libstick-0.1/test.h"
+
+
+void test() {
+    printf("Hello world.\n");
+}
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 66012c5..48fffed 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,9 +1,9 @@
-set(stick_cli_SRC main.c)
+set(stick_cli_SRC main.cc)
 
 configure_file(${CMAKE_CURRENT_SOURCE_DIR}/config.h.cmake ${CMAKE_CURRENT_BINARY_DIR}/config.h)
 include_directories(${CMAKE_CURRENT_BINARY_DIR})
 
-include_directories(${stick_SOURCE_DIR}/include)
+include_directories(${libstick_SOURCE_DIR}/include)
 
 find_package(VTK REQUIRED)
 INCLUDE(${VTK_USE_FILE})
diff --git a/src/main.c b/src/main.c
deleted file mode 100644
index 2c3c4fd..0000000
--- a/src/main.c
+++ /dev/null
@@ -1,13 +0,0 @@
-#include <stdio.h>
-#include <stdlib.h>
-
-#include <stick.h>
-
-int main(int argc, char* argv[]) {
-	(void) argc;
-	(void) argv;
-
-	test();
-	return EXIT_SUCCESS;
-}
-
diff --git a/src/main.cc b/src/main.cc
new file mode 100644
index 0000000..ce34f6f
--- /dev/null
+++ b/src/main.cc
@@ -0,0 +1,13 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+#include <libstick-0.1/test.h>
+
+int main(int argc, char* argv[]) {
+    (void) argc;
+    (void) argv;
+
+    test();
+    return EXIT_SUCCESS;
+}
+
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
new file mode 100644
index 0000000..2ef2b65
--- /dev/null
+++ b/tests/CMakeLists.txt
@@ -0,0 +1,13 @@
+set(tests_SRC main.cc)
+
+INCLUDE (FindPkgConfig)
+
+pkg_search_module(libcpptest REQUIRED libcpptest)
+include_directories(${libcpptest_INCLUDE_DIRS})
+
+include_directories(${libstick_SOURCE_DIR}/include)
+
+add_executable(tests ${tests_SRC})
+target_link_libraries(tests
+	stick
+	${libcpptest_LIBRARIES})
diff --git a/tests/booleanmatrix.h b/tests/booleanmatrix.h
new file mode 100644
index 0000000..8554bf3
--- /dev/null
+++ b/tests/booleanmatrix.h
@@ -0,0 +1,111 @@
+#ifndef booleanmatrix_h_ohYeiveiKanashie
+#define booleanmatrix_h_ohYeiveiKanashie
+
+#include <cpptest.h>
+#include <cpptest-suite.h>
+
+#include <libstick-0.1/booleanmatrix.h>
+
+using namespace libstick;
+
+
+class BooleanmatrixTestSuite: public Test::Suite {
+    public:
+        BooleanmatrixTestSuite() {
+            TEST_ADD(BooleanmatrixTestSuite::test_getsetsize<BooleanRowMatrix<> >);
+            TEST_ADD(BooleanmatrixTestSuite::test_getsetsize<BooleanColMatrix<> >);
+            TEST_ADD(BooleanmatrixTestSuite::test_getsetsize<BooleanColRowMatrix<> >);
+
+            TEST_ADD(BooleanmatrixTestSuite::test_addgetcolumn<BooleanColMatrix<> >);
+            TEST_ADD(BooleanmatrixTestSuite::test_addgetcolumn<BooleanColRowMatrix<> >);
+
+            TEST_ADD(BooleanmatrixTestSuite::test_addgetrow<BooleanRowMatrix<> >);
+            TEST_ADD(BooleanmatrixTestSuite::test_addgetrow<BooleanColRowMatrix<> >);
+        }
+
+    protected:
+        virtual void setup() {
+        }
+
+        virtual void tear_down() {
+        }
+
+    private:
+        template<typename T>
+        void test_getsetsize() {
+            const unsigned size = 17;
+            T mat(size);
+
+            // Must all be zero
+            for (unsigned r=0; r < size; ++r)
+                for (unsigned c=0; c < size; ++c)
+                    TEST_ASSERT(mat.get(r, c) == false);
+
+            // Set a few of them one
+            for (unsigned r=0; r < size; ++r)
+                for (unsigned c=0; c < size; ++c)
+                    if ((257*c + 65537*r + 7)%11 == 0)
+                        mat.set(r, c, true);
+
+            // Check if they are one
+            for (unsigned r=0; r < size; ++r)
+                for (unsigned c=0; c < size; ++c)
+                    TEST_ASSERT(((257*c + 65537*r + 7)%11 == 0) == mat.get(r, c))
+        }
+
+        template<typename T>
+        void test_addgetcolumn() {
+            const unsigned size = 17;
+            T mat(size);
+            typename T::column_type col;
+
+            // Check if columns are empty
+            for (unsigned c=0; c < size; ++c)
+                TEST_ASSERT(mat.getColumn(c).size() == 0);
+
+            col.push_back(3);
+            col.push_back(14);
+            col.push_back(8);
+            col.push_back(0);
+
+            // Add column and test for values
+            mat.add_column(5, col);
+            for (unsigned c=0; c < size; ++c) {
+                if (c==5) {
+                    for (unsigned r=0; r < size; ++r )
+                        TEST_ASSERT(mat.get(r,c) == (r==0 || r==3 || r==8 || r==14));
+                } else {
+                    TEST_ASSERT(mat.getColumn(c).size() == 0);
+                }
+            }
+        }
+
+        template<typename T>
+        void test_addgetrow() {
+            const unsigned size = 17;
+            T mat(size);
+            typename T::row_type row;
+
+            // Check if rows are empty
+            for (unsigned r=0; r < size; ++r)
+                TEST_ASSERT(mat.getRow(r).size() == 0);
+
+            row.push_back(0);
+            row.push_back(8);
+            row.push_back(14);
+            row.push_back(3);
+
+            // Add row and test for values
+            mat.add_row(5, row);
+            for (unsigned r=0; r < size; ++r) {
+                if (r==5) {
+                    for (unsigned c=0; c < size; ++c )
+                        TEST_ASSERT(mat.get(r,c) == (c==0 || c==3 || c==8 || c==14));
+                } else {
+                    TEST_ASSERT(mat.getRow(r).size() == 0);
+                }
+            }
+        }
+};
+
+#endif
diff --git a/tests/main.cc b/tests/main.cc
new file mode 100644
index 0000000..6dca333
--- /dev/null
+++ b/tests/main.cc
@@ -0,0 +1,22 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "booleanmatrix.h"
+#include "simplicialcomplex.h"
+
+using namespace std;
+
+
+int main(int argc, char* argv[]) {
+    (void) argc;
+    (void) argv;
+
+    Test::Suite ts;
+
+    ts.add(auto_ptr<Test::Suite>(new BooleanmatrixTestSuite));
+    ts.add(auto_ptr<Test::Suite>(new SimplicialComplexTestSuite));
+
+    Test::TextOutput output(Test::TextOutput::Verbose);
+    return ts.run(output);
+}
+
diff --git a/tests/simplicialcomplex.h b/tests/simplicialcomplex.h
new file mode 100644
index 0000000..fa829e1
--- /dev/null
+++ b/tests/simplicialcomplex.h
@@ -0,0 +1,168 @@
+#ifndef simplicialcomplex_h_ooDeimaexieghaev
+#define simplicialcomplex_h_ooDeimaexieghaev
+
+#include <cpptest.h>
+#include <cpptest-suite.h>
+
+#include <libstick-0.1/simplicialcomplex.h>
+
+using namespace libstick;
+
+
+class SimplicialComplexTestSuite: public Test::Suite {
+
+    private:
+        typedef SimplicialComplex<2, uint32_t, double> scomplex;
+        typedef scomplex::SimplexOrder::BoundaryMatrix bm;
+
+        bool setupcalled;
+        scomplex c1, c2, c3;
+        scomplex::SimplexOrder o1, o2, o3, o3b;
+
+    public:
+        SimplicialComplexTestSuite() :
+            setupcalled(false),
+            o1(c1),
+            o2(c2),
+            o3(c3),
+            o3b(c3)
+            {
+            TEST_ADD(SimplicialComplexTestSuite::test_isComplex);
+            TEST_ADD(SimplicialComplexTestSuite::test_isMonotoneComplex);
+            TEST_ADD(SimplicialComplexTestSuite::test_isOrderFiltration);
+            TEST_ADD(SimplicialComplexTestSuite::test_isOrderMonotone);
+            TEST_ADD(SimplicialComplexTestSuite::test_boundaryMatrix);
+        }
+
+    protected:
+        virtual void setup() {
+            if (setupcalled)
+                return;
+            setupcalled = true;
+
+            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}
+            };
+
+            //  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
+
+            // Build the complex
+            for (unsigned i=0; i<num; ++i)
+                c1.addSimplex(ss[i]);
+
+            c2 = c1;
+            c2.simplices[5].value = 11;
+
+            c3 = c2;
+            c3.simplices[10].value = 12;
+
+            o1.reset();
+            o2.reset();
+            o3.reset();
+            o3b.reset();
+            o3b.makeMonotoneFiltration();
+        }
+
+        virtual void tear_down() {
+        }
+
+        void test_isComplex() {
+            TEST_ASSERT(c1.isComplex());
+            TEST_ASSERT(c2.isComplex());
+            TEST_ASSERT(c3.isComplex());
+        }
+
+        void test_isMonotoneComplex() {
+            TEST_ASSERT(c1.isMonotone());
+            TEST_ASSERT(!c2.isMonotone());
+            TEST_ASSERT(c3.isMonotone());
+        }
+
+        void test_isOrderFiltration() {
+            TEST_ASSERT(o1.isFiltration());
+            TEST_ASSERT(o2.isFiltration());
+            TEST_ASSERT(o3.isFiltration());
+            TEST_ASSERT(o3b.isFiltration());
+        }
+
+        void test_isOrderMonotone() {
+            TEST_ASSERT(o1.isMonotone());
+            TEST_ASSERT(!o2.isMonotone());
+            TEST_ASSERT(!o3.isMonotone());
+            TEST_ASSERT(o3b.isMonotone());
+        }
+
+        void test_boundaryMatrix() {
+            bm mat1 = o1.getBoundaryMatrix();
+            bm mat1e(c1.size());
+            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(8, 10, true);
+            TEST_ASSERT(mat1 == mat1e);
+
+            bm mat3b = o3b.getBoundaryMatrix();
+            bm mat3be(c1.size());
+            mat3be.set(0, 4, true);
+            mat3be.set(1, 4, 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);
+
+            //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;
+
+            TEST_ASSERT(mat3b == mat3be);
+        }
+};
+
+#endif