Add doc, examples and unitary tests

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/get_persistence@1213 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: dc7980ef3f1baef17f2e42770da146d6c961e672

7 files changed, 422 insertions, 61 deletions

diff --git a/src/Persistent_cohomology/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt
index 1f32da17..186a6c33 100644
--- a/src/Persistent_cohomology/example/CMakeLists.txt
+++ b/src/Persistent_cohomology/example/CMakeLists.txt
@@ -74,12 +74,17 @@ if(GMPXX_FOUND AND GMP_FOUND)
         add_executable(periodic_alpha_complex_3d_persistence periodic_alpha_complex_3d_persistence.cpp)
         target_link_libraries(periodic_alpha_complex_3d_persistence ${Boost_SYSTEM_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES} ${CGAL_LIBRARY})
 
+        add_executable(custom_persistence_sort custom_persistence_sort.cpp)
+        target_link_libraries(custom_persistence_sort ${Boost_SYSTEM_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES} ${CGAL_LIBRARY})
+
         if (TBB_FOUND)
           target_link_libraries(alpha_complex_persistence ${TBB_RELEASE_LIBRARY})
           target_link_libraries(periodic_alpha_complex_3d_persistence ${TBB_RELEASE_LIBRARY})
+          target_link_libraries(custom_persistence_sort ${TBB_RELEASE_LIBRARY})
         endif()
         add_test(alpha_complex_persistence_2_0_45 ${CMAKE_CURRENT_BINARY_DIR}/alpha_complex_persistence ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off -m 0.45 -p 2)
         add_test(periodic_alpha_complex_3d_persistence_2_0 ${CMAKE_CURRENT_BINARY_DIR}/periodic_alpha_complex_3d_persistence ${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off 2 0)
+        add_test(custom_persistence_sort ${CMAKE_CURRENT_BINARY_DIR}/custom_persistence_sort)
 
       else()
         message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.")
diff --git a/src/Persistent_cohomology/example/custom_persistence_sort.cpp b/src/Persistent_cohomology/example/custom_persistence_sort.cpp
new file mode 100644
index 00000000..4dab4560
--- /dev/null
+++ b/src/Persistent_cohomology/example/custom_persistence_sort.cpp
@@ -0,0 +1,115 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2014  INRIA Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <CGAL/Epick_d.h>
+#include <CGAL/point_generators_d.h>
+#include <CGAL/algorithm.h>
+#include <CGAL/assertions.h>
+
+#include <gudhi/Alpha_complex.h>
+#include <gudhi/Persistent_cohomology.h>
+
+#include <iostream>
+#include <iterator>
+#include <vector>
+#include <fstream>  // for std::ofstream
+
+
+using Kernel = CGAL::Epick_d< CGAL::Dimension_tag<3> >;
+using Point = Kernel::Point_d;
+using Alpha_complex = Gudhi::alphacomplex::Alpha_complex<Kernel>;
+
+std::vector<Point> random_points() {
+  // Instanciate a random point generator
+  CGAL::Random rng(0);
+
+  // Generate "points_number" random points in a vector
+  std::vector<Point> points;
+
+  // Generates 1000 random 3D points on a sphere of radius 4.0
+  CGAL::Random_points_on_sphere_d<Point> rand_outside(3, 4.0, rng);
+  CGAL::cpp11::copy_n(rand_outside, 1000, std::back_inserter(points));
+  // Generates 2000 random 3D points in a sphere of radius 3.0
+  CGAL::Random_points_in_ball_d<Point> rand_inside(3, 3.0, rng);
+  CGAL::cpp11::copy_n(rand_inside, 2000, std::back_inserter(points));
+
+  return points;
+}
+
+/*
+ * Compare two intervals by dimension, then by length.
+ */
+struct cmp_intervals_by_dim_then_length {
+  explicit cmp_intervals_by_dim_then_length(Alpha_complex * sc)
+      : sc_(sc) { }
+
+  template<typename Persistent_interval>
+  bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) {
+    if (sc_->dimension(get < 0 > (p1)) == sc_->dimension(get < 0 > (p2)))
+      return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1))
+              > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2)));
+    else
+      return (sc_->dimension(get < 0 > (p1)) > sc_->dimension(get < 0 > (p2)));
+  }
+  Alpha_complex* sc_;
+};
+
+int main(int argc, char **argv) {
+  std::vector<Point> points = random_points();
+
+  // Alpha complex persistence computation from generated points
+  Alpha_complex alpha_complex_from_points(points, 0.6);
+
+  using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology< Alpha_complex,
+      Gudhi::persistent_cohomology::Field_Zp >;
+  Persistent_cohomology pcoh(alpha_complex_from_points);
+
+  // initializes the coefficient field for homology - Z/3Z
+  pcoh.init_coefficients(3);
+  pcoh.compute_persistent_cohomology(0.2);
+
+  // Custom sort and output persistence
+  cmp_intervals_by_dim_then_length cmp(&alpha_complex_from_points);
+  auto persistent_pairs = pcoh.get_persistent_pairs();
+  std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp);
+  for (auto pair : persistent_pairs) {
+    std::cout << alpha_complex_from_points.dimension(get<0>(pair)) << " "
+          << alpha_complex_from_points.filtration(get<0>(pair)) << " "
+          << alpha_complex_from_points.filtration(get<1>(pair)) << std::endl;
+  }
+
+  // Persistent Betti numbers
+  std::cout << "The persistent Betti numbers in interval [0.40, 0.41] are : ";
+  for (int dim = 0; dim < alpha_complex_from_points.dimension(); dim++)
+    std::cout << "b" << dim << " = " << pcoh.persistent_betti_number(dim, 0.40, 0.41) << " ; ";
+  std::cout << std::endl;
+
+  // Betti numbers
+  std::vector<int> betti_numbers = pcoh.betti_numbers();
+  std::cout << "The Betti numbers are : ";
+  for (std::size_t i = 0; i < betti_numbers.size(); i++)
+    std::cout << "b" << i << " = " << betti_numbers[i] << " ; ";
+  std::cout << std::endl;
+
+  return 0;
+}
+
diff --git a/src/Persistent_cohomology/example/plain_homology.cpp b/src/Persistent_cohomology/example/plain_homology.cpp
index 9e5adb5d..5afce9e2 100644
--- a/src/Persistent_cohomology/example/plain_homology.cpp
+++ b/src/Persistent_cohomology/example/plain_homology.cpp
@@ -42,24 +42,6 @@ struct MyOptions : Simplex_tree_options_full_featured {
 };
 typedef Simplex_tree<MyOptions> ST;
 
-  /*
-   * Compare two intervals by dimension, then by length.
-   */
-  struct cmp_intervals_by_dim_then_length {
-    explicit cmp_intervals_by_dim_then_length(ST * sc)
-        : sc_(sc) {
-    }
-    template<typename Persistent_interval>
-    bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) {
-      if (sc_->dimension(get < 0 > (p1)) == sc_->dimension(get < 0 > (p2)))
-        return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1))
-            > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2)));
-      else
-        return (sc_->dimension(get < 0 > (p1)) > sc_->dimension(get < 0 > (p2)));
-    }
-    ST* sc_;
-  };
-
 int main() {
   ST st;
 
@@ -102,34 +84,10 @@ int main() {
   pcoh.output_diagram();
   
   
-  // ********************************************************
-  // get_persistence
-  // ********************************************************
-  std::cout << std::endl;
+  // Print the Betti numbers are b0=2 and b1=1.
   std::cout << std::endl;
-
-  // First version
-  std::vector<int> betti_numbers = pcoh.betti_numbers();
-  std::cout << "The Betti numbers are : ";
-  for (std::size_t i = 0; i < betti_numbers.size(); i++)
-    std::cout << "b" << i << " = " << betti_numbers[i] << " ; ";
-  std::cout << std::endl;
-
-  // Second version
   std::cout << "The Betti numbers are : ";
   for (int i = 0; i < st.dimension(); i++)
     std::cout << "b" << i << " = " << pcoh.betti_number(i) << " ; ";
   std::cout << std::endl;
-  
-  // Get persistence
-  cmp_intervals_by_dim_then_length cmp(&st);
-  auto persistent_pairs = pcoh.get_persistent_pairs();
-  std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp);
-  for (auto pair : persistent_pairs) {
-    std::cout << st.dimension(get<0>(pair)) << " "
-          << st.filtration(get<0>(pair)) << " "
-          << st.filtration(get<1>(pair)) << std::endl;
-  }
-
-
 }
diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h
index f0cb159e..5b4c215a 100644
--- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h
+++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h
@@ -88,7 +88,7 @@ class Persistent_cohomology {
   /** \brief Initializes the Persistent_cohomology class.
    *
    * @param[in] cpx Complex for which the persistent homology is computed.
-   cpx is a model of FilteredComplex
+   * cpx is a model of FilteredComplex
    */
   explicit Persistent_cohomology(Complex_ds& cpx)
       : cpx_(&cpx),
@@ -117,7 +117,7 @@ class Persistent_cohomology {
   /** \brief Initializes the Persistent_cohomology class.
    *
    * @param[in] cpx Complex for which the persistent homology is compiuted.
-   cpx is a model of FilteredComplex
+   * cpx is a model of FilteredComplex
    *
    * @param[in] persistence_dim_max if true, the persistent homology for the maximal dimension in the
    *                                complex is computed. If false, it is ignored. Default is false.
@@ -554,14 +554,6 @@ class Persistent_cohomology {
     Complex_ds * sc_;
   };
 
-  /*
-   * Returns true when Filtration_value type accepts infinity values and the given value is equal to infinity.
-   */
-  bool is_infinity(Filtration_value value) const {
-    bool has_infinity = std::numeric_limits<Filtration_value>::has_infinity;
-    return (has_infinity && value == std::numeric_limits<Filtration_value>::infinity());
-  }
-
  public:
   /** \brief Output the persistence diagram in ostream.
    *
@@ -576,13 +568,14 @@ class Persistent_cohomology {
   void output_diagram(std::ostream& ostream = std::cout) {
     cmp_intervals_by_length cmp(cpx_);
     std::sort(std::begin(persistent_pairs_), std::end(persistent_pairs_), cmp);
+    bool has_infinity = std::numeric_limits<Filtration_value>::has_infinity;
     for (auto pair : persistent_pairs_) {
       // Special case on windows, inf is "1.#INF" (cf. unitary tests and R package TDA)
-      if (is_infinity(cpx_->filtration(get<1>(pair)))) {
-        ostream << /*get<2>(pair) <<*/ "  " << cpx_->dimension(get<0>(pair)) << " "
+      if (has_infinity && cpx_->filtration(get<1>(pair)) == std::numeric_limits<Filtration_value>::infinity()) {
+        ostream << get<2>(pair) << "  " << cpx_->dimension(get<0>(pair)) << " "
           << cpx_->filtration(get<0>(pair)) << " inf " << std::endl;
       } else {
-        ostream << /*get<2>(pair) <<*/ "  " << cpx_->dimension(get<0>(pair)) << " "
+        ostream << get<2>(pair) << "  " << cpx_->dimension(get<0>(pair)) << " "
           << cpx_->filtration(get<0>(pair)) << " "
           << cpx_->filtration(get<1>(pair)) << " " << std::endl;
       }
@@ -601,14 +594,15 @@ class Persistent_cohomology {
   }
 
   /** @brief Returns Betti numbers.
-   *
+   * @return A vector of persistent Betti numbers.
    */
   std::vector<int> betti_numbers() const {
     // Init Betti numbers vector with zeros until Simplicial complex dimension
     std::vector<int> betti_numbers(cpx_->dimension(), 0);
     
     for (auto pair : persistent_pairs_) {
-      if (is_infinity(cpx_->filtration(get<1>(pair)))) {
+      // Count never ended persistence intervals
+      if (cpx_->null_simplex() == get<1>(pair)) {
         // Increment corresponding betti number
         betti_numbers[cpx_->dimension(get<0>(pair))] += 1;
       }
@@ -616,16 +610,56 @@ class Persistent_cohomology {
     return betti_numbers;
   }
 
-  /** @brief Returns the Betti number passed by parameter.
+  /** @brief Returns the Betti number of the dimension passed by parameter.
    * @param[in] dimension The Betti number dimension to get.
-   * @return Betti number
+   * @return Betti number of the given dimension
    *
    */
   int betti_number(int dimension) const {
     int betti_number = 0;
     
     for (auto pair : persistent_pairs_) {
-      if (is_infinity(cpx_->filtration(get<1>(pair)))) {
+      // Count never ended persistence intervals
+      if (cpx_->null_simplex() == get<1>(pair)) {
+        if (cpx_->dimension(get<0>(pair)) == dimension) {
+          // Increment betti number found
+          ++betti_number;
+        }
+      }
+    }
+    return betti_number;
+  }
+  
+  /** @brief Returns the persistent Betti numbers.
+   * @param[in] from The persistence birth limit to be added in the number \f$(persistent birth \leq from)\f$.
+   * @param[in] to The persistence death limit to be added in the number  \f$(persistent death > from)\f$.
+   * @return A vector of persistent Betti numbers.
+   */
+  std::vector<int> persistent_betti_numbers(Filtration_value from, Filtration_value to) const {
+    // Init Betti numbers vector with zeros until Simplicial complex dimension
+    std::vector<int> betti_numbers(cpx_->dimension(), 0);
+    
+    for (auto pair : persistent_pairs_) {
+      // Count persistence intervals that covers the given interval
+      if (cpx_->filtration(get<0>(pair)) <= from && cpx_->filtration(get<1>(pair)) > to) {
+        // Increment corresponding betti number
+        betti_numbers[cpx_->dimension(get<0>(pair))] += 1;
+      }
+    }
+    return betti_numbers;
+  }
+
+  /** @brief Returns the persistentBetti number of the dimension passed by parameter.
+   * @param[in] dimension The Betti number dimension to get.
+   * @return Betti number of the given dimension
+   *
+   */
+  int persistent_betti_number(int dimension, Filtration_value from, Filtration_value to) const {
+    int betti_number = 0;
+    
+    for (auto pair : persistent_pairs_) {
+      // Count persistence intervals that covers the given interval
+      if (cpx_->filtration(get<0>(pair)) <= from && cpx_->filtration(get<1>(pair)) > to) {
         if (cpx_->dimension(get<0>(pair)) == dimension) {
           // Increment betti number found
           ++betti_number;
diff --git a/src/Persistent_cohomology/test/CMakeLists.txt b/src/Persistent_cohomology/test/CMakeLists.txt
index 459cc000..a034031a 100644
--- a/src/Persistent_cohomology/test/CMakeLists.txt
+++ b/src/Persistent_cohomology/test/CMakeLists.txt
@@ -12,8 +12,11 @@ endif()
 
 add_executable ( PersistentCohomologyUT persistent_cohomology_unit_test.cpp )
 target_link_libraries(PersistentCohomologyUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+add_executable ( BettiNumbersUT betti_numbers_unit_test.cpp )
+target_link_libraries(BettiNumbersUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
 if (TBB_FOUND)
   target_link_libraries(PersistentCohomologyUT ${TBB_RELEASE_LIBRARY})
+  target_link_libraries(BettiNumbersUT ${TBB_RELEASE_LIBRARY})
 endif()
 
 # Unitary tests
@@ -23,6 +26,11 @@ add_test(NAME PersistentCohomologyUT
     # XML format for Jenkins xUnit plugin 
     --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/PersistentCohomologyUT.xml --log_level=test_suite --report_level=no)
 
+add_test(NAME BettiNumbersUT 
+  COMMAND ${CMAKE_CURRENT_BINARY_DIR}/BettiNumbersUT
+    # XML format for Jenkins xUnit plugin 
+    --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/BettiNumbersUT.xml --log_level=test_suite --report_level=no)
+
 if(GMPXX_FOUND AND GMP_FOUND)
   add_executable ( PersistentCohomologyMultiFieldUT persistent_cohomology_unit_test_multi_field.cpp )
 target_link_libraries(PersistentCohomologyMultiFieldUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES})
diff --git a/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp
new file mode 100644
index 00000000..a4e00b45
--- /dev/null
+++ b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp
@@ -0,0 +1,234 @@
+#include <iostream>
+#include <string>
+#include <algorithm>
+#include <utility> // std::pair, std::make_pair
+#include <cmath> // float comparison
+#include <limits>
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "betti_numbers"
+#include <boost/test/unit_test.hpp>
+
+#include "gudhi/Simplex_tree.h"
+#include "gudhi/Persistent_cohomology.h"
+
+struct MyOptions : Gudhi::Simplex_tree_options_full_featured {
+  // Implicitly use 0 as filtration value for all simplices
+  static const bool store_filtration = false;
+  // The persistence algorithm needs this
+  static const bool store_key = true;
+  // I have few vertices
+  typedef short Vertex_handle;
+};
+
+using Mini_simplex_tree = Gudhi::Simplex_tree<MyOptions>;
+using Mini_st_persistence =
+    Gudhi::persistent_cohomology::Persistent_cohomology<Mini_simplex_tree, Gudhi::persistent_cohomology::Field_Zp>;
+
+/*
+ * Compare two intervals by dimension, then by length.
+ */
+template<class Simplicial_complex>
+struct cmp_intervals_by_dim_then_length {
+  explicit cmp_intervals_by_dim_then_length(Simplicial_complex * sc)
+      : sc_(sc) { }
+
+  template<typename Persistent_interval>
+  bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) {
+    if (sc_->dimension(get < 0 > (p1)) == sc_->dimension(get < 0 > (p2)))
+      return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1))
+              > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2)));
+    else
+      return (sc_->dimension(get < 0 > (p1)) > sc_->dimension(get < 0 > (p2)));
+  }
+  Simplicial_complex* sc_;
+};
+
+BOOST_AUTO_TEST_CASE( plain_homology_betti_numbers )
+{
+  Mini_simplex_tree st;
+
+  /* Complex to build. */
+  /*    1   4          */
+  /*    o---o          */
+  /*   /3\ /           */
+  /*  o---o   o        */
+  /*  2   0   5        */
+  const short tetra0123[] = {0, 1, 2, 3};
+  const short edge04[] = {0, 4};
+  const short edge14[] = {1, 4};
+  const short vertex5[] = {5};
+  st.insert_simplex_and_subfaces(tetra0123);
+  st.insert_simplex_and_subfaces(edge04);
+  st.insert_simplex(edge14);
+  st.insert_simplex(vertex5);
+  // FIXME: Remove this line
+  st.set_dimension(3);
+
+  // Sort the simplices in the order of the filtration
+  st.initialize_filtration();
+
+  // Class for homology computation
+  Mini_st_persistence pcoh(st);
+
+  // Initialize the coefficient field Z/3Z for homology
+  pcoh.init_coefficients(3);
+
+  // Compute the persistence diagram of the complex
+  pcoh.compute_persistent_cohomology();
+
+  // Print the result. The format is, on each line: 2 dim 0 inf
+  // where 2 represents the field, dim the dimension of the feature.
+  // 2  0 0 inf 
+  // 2  0 0 inf 
+  // 2  1 0 inf 
+  // means that in Z/2Z-homology, the Betti numbers are b0=2 and b1=1.
+  
+  BOOST_CHECK(pcoh.betti_number(0) == 2);
+  BOOST_CHECK(pcoh.betti_number(1) == 1);
+  BOOST_CHECK(pcoh.betti_number(2) == 0);
+
+  std::vector<int> bns = pcoh.betti_numbers();
+  BOOST_CHECK(bns.size() == 3);
+  BOOST_CHECK(bns[0] == 2);
+  BOOST_CHECK(bns[1] == 1);
+  BOOST_CHECK(bns[2] == 0);
+  
+  // Custom sort and output persistence
+  cmp_intervals_by_dim_then_length<Mini_simplex_tree> cmp(&st);
+  auto persistent_pairs = pcoh.get_persistent_pairs();
+  
+  std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp);
+  
+  BOOST_CHECK(persistent_pairs.size() == 3);
+  // persistent_pairs[0] = 2  1 0 inf
+  BOOST_CHECK(st.dimension(get<0>(persistent_pairs[0])) == 1);
+  BOOST_CHECK(st.filtration(get<0>(persistent_pairs[0])) == 0);
+  BOOST_CHECK(get<1>(persistent_pairs[0]) == st.null_simplex());
+
+  // persistent_pairs[1] = 2  0 0 inf
+  BOOST_CHECK(st.dimension(get<0>(persistent_pairs[1])) == 0);
+  BOOST_CHECK(st.filtration(get<0>(persistent_pairs[1])) == 0);
+  BOOST_CHECK(get<1>(persistent_pairs[1]) == st.null_simplex());
+
+  // persistent_pairs[2] = 2  0 0 inf
+  BOOST_CHECK(st.dimension(get<0>(persistent_pairs[2])) == 0);
+  BOOST_CHECK(st.filtration(get<0>(persistent_pairs[2])) == 0);
+  BOOST_CHECK(get<1>(persistent_pairs[2]) == st.null_simplex());
+}
+
+using Simplex_tree = Gudhi::Simplex_tree<>;
+using St_persistence =
+    Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp>;
+
+BOOST_AUTO_TEST_CASE( betti_numbers )
+{
+  Simplex_tree st;
+
+  /* Complex to build. */
+  /*    1   4          */
+  /*    o---o          */
+  /*   /3\ /           */
+  /*  o---o   o        */
+  /*  2   0   5        */
+  const short tetra0123[] = {0, 1, 2, 3};
+  const short edge04[] = {0, 4};
+  const short edge14[] = {1, 4};
+  const short vertex5[] = {5};
+  st.insert_simplex_and_subfaces(tetra0123, 4.0);
+  st.insert_simplex_and_subfaces(edge04, 2.0);
+  st.insert_simplex(edge14, 2.0);
+  st.insert_simplex(vertex5, 1.0);
+  // FIXME: Remove this line
+  st.set_dimension(3);
+
+  // Sort the simplices in the order of the filtration
+  st.initialize_filtration();
+
+  // Class for homology computation
+  St_persistence pcoh(st);
+
+  // Initialize the coefficient field Z/3Z for homology
+  pcoh.init_coefficients(3);
+
+  // Compute the persistence diagram of the complex
+  pcoh.compute_persistent_cohomology();
+
+  // Check the Betti numbers are b0=2, b1=1 and b2=0.
+  BOOST_CHECK(pcoh.betti_number(0) == 2);
+  BOOST_CHECK(pcoh.betti_number(1) == 1);
+  BOOST_CHECK(pcoh.betti_number(2) == 0);
+
+  // Check the Betti numbers are b0=2, b1=1 and b2=0.
+  std::vector<int> bns = pcoh.betti_numbers();
+  BOOST_CHECK(bns.size() == 3);
+  BOOST_CHECK(bns[0] == 2);
+  BOOST_CHECK(bns[1] == 1);
+  BOOST_CHECK(bns[2] == 0);
+  
+  // Check the persistent Betti numbers in [4., 10.] are b0=2, b1=1 and b2=0.
+  BOOST_CHECK(pcoh.persistent_betti_number(0, 4., 10.) == 2);
+  BOOST_CHECK(pcoh.persistent_betti_number(1, 4., 10.) == 1);
+  BOOST_CHECK(pcoh.persistent_betti_number(2, 4., 10.) == 0);
+
+  // Check the persistent Betti numbers in [2., 100.] are b0=2, b1=0 and b2=0.
+  BOOST_CHECK(pcoh.persistent_betti_number(0, 2., 100.) == 2);
+  BOOST_CHECK(pcoh.persistent_betti_number(1, 2., 100.) == 0);
+  BOOST_CHECK(pcoh.persistent_betti_number(2, 2., 100.) == 0);
+
+  // Check the persistent Betti numbers in [1., 1000.] are b0=1, b1=0 and b2=0.
+  BOOST_CHECK(pcoh.persistent_betti_number(0, 1., 1000.) == 1);
+  BOOST_CHECK(pcoh.persistent_betti_number(1, 1., 1000.) == 0);
+  BOOST_CHECK(pcoh.persistent_betti_number(2, 1., 1000.) == 0);
+
+  // Check the persistent Betti numbers in [.9, 1000.] are b0=0, b1=0 and b2=0.
+  BOOST_CHECK(pcoh.persistent_betti_number(0, .9, 1000.) == 0);
+  BOOST_CHECK(pcoh.persistent_betti_number(1, .9, 1000.) == 0);
+  BOOST_CHECK(pcoh.persistent_betti_number(2, .9, 1000.) == 0);
+
+  // Check the persistent Betti numbers in [4.1, 10000.] are b0=2, b1=1 and b2=0.
+  bns = pcoh.persistent_betti_numbers(4.1, 10000.);
+  BOOST_CHECK(bns[0] == 2);
+  BOOST_CHECK(bns[1] == 1);
+  BOOST_CHECK(bns[2] == 0);
+  
+  // Check the persistent Betti numbers in [2.1, 100000.] are b0=2, b1=0 and b2=0.
+  bns = pcoh.persistent_betti_numbers(2.1, 100000.);
+  BOOST_CHECK(bns[0] == 2);
+  BOOST_CHECK(bns[1] == 0);
+  BOOST_CHECK(bns[2] == 0);
+
+  // Check the persistent Betti numbers in [1.1, 1000000.] are b0=1, b1=0 and b2=0.
+  bns = pcoh.persistent_betti_numbers(1.1, 1000000.);
+  BOOST_CHECK(bns[0] == 1);
+  BOOST_CHECK(bns[1] == 0);
+  BOOST_CHECK(bns[2] == 0);
+  
+  // Check the persistent Betti numbers in [.1, 10000000.] are b0=0, b1=0 and b2=0.
+  bns = pcoh.persistent_betti_numbers(.1, 10000000.);
+  BOOST_CHECK(bns[0] == 0);
+  BOOST_CHECK(bns[1] == 0);
+  BOOST_CHECK(bns[2] == 0);
+  
+  // Custom sort and output persistence
+  cmp_intervals_by_dim_then_length<Simplex_tree> cmp(&st);
+  auto persistent_pairs = pcoh.get_persistent_pairs();
+  
+  std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp);
+  
+  BOOST_CHECK(persistent_pairs.size() == 3);
+  // persistent_pairs[0] = 2  1 4 inf
+  BOOST_CHECK(st.dimension(get<0>(persistent_pairs[0])) == 1);
+  BOOST_CHECK(st.filtration(get<0>(persistent_pairs[0])) == 4);
+  BOOST_CHECK(get<1>(persistent_pairs[0]) == st.null_simplex());
+
+  // persistent_pairs[1] = 2  0 2 inf
+  BOOST_CHECK(st.dimension(get<0>(persistent_pairs[1])) == 0);
+  BOOST_CHECK(st.filtration(get<0>(persistent_pairs[1])) == 2);
+  BOOST_CHECK(get<1>(persistent_pairs[1]) == st.null_simplex());
+
+  // persistent_pairs[2] = 2  0 1 inf
+  BOOST_CHECK(st.dimension(get<0>(persistent_pairs[2])) == 0);
+  BOOST_CHECK(st.filtration(get<0>(persistent_pairs[2])) == 1);
+  BOOST_CHECK(get<1>(persistent_pairs[2]) == st.null_simplex());
+}
diff --git a/src/common/doc/main_page.h b/src/common/doc/main_page.h
index 19dece4b..2b9f2f40 100644
--- a/src/common/doc/main_page.h
+++ b/src/common/doc/main_page.h
@@ -209,6 +209,8 @@
  * Persistent_cohomology/alpha_complex_persistence.cpp</a>
  * \li <a href="_persistent_cohomology_2periodic_alpha_complex_3d_persistence_8cpp-example.html">
  * Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp</a>
+ * \li <a href="_persistent_cohomology_2custom_persistence_sort_8cpp-example.html">
+ * Persistent_cohomology/custom_persistence_sort.cpp</a>
  * 
  * \subsection eigen3 Eigen3:
  * <a target="_blank" href="http://eigen.tuxfamily.org/">Eigen3</a> is a C++ template library for linear algebra:
@@ -224,6 +226,8 @@
  * Persistent_cohomology/alpha_complex_persistence.cpp</a>
  * \li <a href="_persistent_cohomology_2periodic_alpha_complex_3d_persistence_8cpp-example.html">
  * Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp</a>
+ * \li <a href="_persistent_cohomology_2custom_persistence_sort_8cpp-example.html">
+ * Persistent_cohomology/custom_persistence_sort.cpp</a>
  * 
  * \subsection tbb Threading Building Blocks:
  * <a target="_blank" href="https://www.threadingbuildingblocks.org/">Intel&reg; TBB</a> lets you easily write parallel
@@ -271,6 +275,8 @@
  * Persistent_cohomology/rips_persistence.cpp</a>
  * \li <a href="_persistent_cohomology_2periodic_alpha_complex_3d_persistence_8cpp-example.html">
  * Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp</a>
+ * \li <a href="_persistent_cohomology_2custom_persistence_sort_8cpp-example.html">
+ * Persistent_cohomology/custom_persistence_sort.cpp</a>
  * 
  * \subsection demos Demos and examples
  * To build the demos and examples, run the following commands in a terminal:
@@ -325,6 +331,7 @@ make \endverbatim
  * @example Persistent_cohomology/plain_homology.cpp
  * @example Persistent_cohomology/rips_multifield_persistence.cpp
  * @example Persistent_cohomology/rips_persistence.cpp
+ * @example Persistent_cohomology/custom_persistence_sort.cpp
  * @example Simplex_tree/mini_simplex_tree.cpp
  * @example Simplex_tree/simple_simplex_tree.cpp
  * @example Simplex_tree/simplex_tree_from_alpha_shapes_3.cpp