3 files changed, 132 insertions, 3 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..d591a0d0
--- /dev/null
+++ b/src/Persistent_cohomology/example/custom_persistence_sort.cpp
@@ -0,0 +1,116 @@
+/*    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
+#include <algorithm>  // for std::sort
+
+
+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 0a692717..2156ffd9 100644
--- a/src/Persistent_cohomology/example/plain_homology.cpp
+++ b/src/Persistent_cohomology/example/plain_homology.cpp
@@ -24,6 +24,7 @@
 #include <gudhi/Persistent_cohomology.h>
 
 #include <iostream>
+#include <vector>
 
 using namespace Gudhi;
 
@@ -76,9 +77,16 @@ int main() {
 
   // 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 
+  // 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.
   pcoh.output_diagram();
+
+  // Print the Betti numbers are b0=2 and b1=1.
+  std::cout << std::endl;
+  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;
 }