GUDHI 2.1.0 as released by upstream in a tarball.upstream/2.1.0

70 files changed, 6086 insertions, 19 deletions

diff --git a/utilities/Alpha_complex/CMakeLists.txt b/utilities/Alpha_complex/CMakeLists.txt
new file mode 100644
index 00000000..a2dfac20
--- /dev/null
+++ b/utilities/Alpha_complex/CMakeLists.txt
@@ -0,0 +1,65 @@
+cmake_minimum_required(VERSION 2.6)
+project(Alpha_complex_utilities)
+
+if(CGAL_FOUND)
+  add_executable(alpha_complex_3d_persistence alpha_complex_3d_persistence.cpp)
+  target_link_libraries(alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+  add_executable(exact_alpha_complex_3d_persistence exact_alpha_complex_3d_persistence.cpp)
+  target_link_libraries(exact_alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+  add_executable(weighted_alpha_complex_3d_persistence weighted_alpha_complex_3d_persistence.cpp)
+  target_link_libraries(weighted_alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+
+  if (TBB_FOUND)
+    target_link_libraries(alpha_complex_3d_persistence ${TBB_LIBRARIES})
+    target_link_libraries(exact_alpha_complex_3d_persistence ${TBB_LIBRARIES})
+    target_link_libraries(weighted_alpha_complex_3d_persistence ${TBB_LIBRARIES})
+  endif(TBB_FOUND)
+  
+  add_test(NAME Alpha_complex_utilities_alpha_complex_3d_persistence COMMAND $<TARGET_FILE:alpha_complex_3d_persistence>
+      "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45")
+  add_test(NAME Alpha_complex_utilities_exact_alpha_complex_3d COMMAND $<TARGET_FILE:exact_alpha_complex_3d_persistence>
+      "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off"  "-p" "2" "-m" "0.45")
+  add_test(NAME Alpha_complex_utilities_weighted_alpha_complex_3d COMMAND $<TARGET_FILE:weighted_alpha_complex_3d_persistence>
+      "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.weights"  "-p" "2" "-m" "0.45")
+
+  install(TARGETS alpha_complex_3d_persistence DESTINATION bin)
+  install(TARGETS exact_alpha_complex_3d_persistence DESTINATION bin)
+  install(TARGETS weighted_alpha_complex_3d_persistence DESTINATION bin)
+
+  if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0)
+    add_executable (alpha_complex_persistence alpha_complex_persistence.cpp)
+    target_link_libraries(alpha_complex_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+
+    add_executable(periodic_alpha_complex_3d_persistence periodic_alpha_complex_3d_persistence.cpp)
+    target_link_libraries(periodic_alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+
+    if (TBB_FOUND)
+      target_link_libraries(alpha_complex_persistence ${TBB_LIBRARIES})
+      target_link_libraries(periodic_alpha_complex_3d_persistence ${TBB_LIBRARIES})
+    endif(TBB_FOUND)
+    add_test(NAME Alpha_complex_utilities_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence>
+        "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45")
+    add_test(NAME Alpha_complex_utilities_periodic_alpha_complex_3d_persistence COMMAND $<TARGET_FILE:periodic_alpha_complex_3d_persistence>
+        "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt" "-p" "2" "-m" "0")
+
+    install(TARGETS alpha_complex_persistence DESTINATION bin)
+    install(TARGETS periodic_alpha_complex_3d_persistence DESTINATION bin)
+
+  endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0)
+
+  if (NOT CGAL_VERSION VERSION_LESS 4.11.0)
+    add_executable(weighted_periodic_alpha_complex_3d_persistence weighted_periodic_alpha_complex_3d_persistence.cpp)
+    target_link_libraries(weighted_periodic_alpha_complex_3d_persistence ${CGAL_LIBRARY})
+    if (TBB_FOUND)
+      target_link_libraries(weighted_periodic_alpha_complex_3d_persistence ${TBB_LIBRARIES})
+    endif(TBB_FOUND)
+
+    add_test(NAME Alpha_complex_utilities_weigted_periodic_alpha_complex_3d COMMAND $<TARGET_FILE:weighted_periodic_alpha_complex_3d_persistence>
+        "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.weights"
+        "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt" "3" "1.0")
+
+    install(TARGETS weighted_periodic_alpha_complex_3d_persistence DESTINATION bin)
+
+  endif (NOT CGAL_VERSION VERSION_LESS 4.11.0)
+
+endif(CGAL_FOUND)
diff --git a/utilities/Alpha_complex/alpha_complex_3d_helper.h b/utilities/Alpha_complex/alpha_complex_3d_helper.h
new file mode 100644
index 00000000..a59f0654
--- /dev/null
+++ b/utilities/Alpha_complex/alpha_complex_3d_helper.h
@@ -0,0 +1,74 @@
+/*    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/>.
+ */
+
+#ifndef ALPHA_COMPLEX_3D_HELPER_H_
+#define ALPHA_COMPLEX_3D_HELPER_H_
+
+template <class Vertex_list, class Cell_handle>
+Vertex_list from_cell(const Cell_handle& ch) {
+  Vertex_list the_list;
+  for (auto i = 0; i < 4; i++) {
+#ifdef DEBUG_TRACES
+    std::cout << "from cell[" << i << "]=" << ch->vertex(i)->point() << std::endl;
+#endif  // DEBUG_TRACES
+    the_list.push_back(ch->vertex(i));
+  }
+  return the_list;
+}
+
+template <class Vertex_list, class Facet>
+Vertex_list from_facet(const Facet& fct) {
+  Vertex_list the_list;
+  for (auto i = 0; i < 4; i++) {
+    if (fct.second != i) {
+#ifdef DEBUG_TRACES
+      std::cout << "from facet=[" << i << "]" << fct.first->vertex(i)->point() << std::endl;
+#endif  // DEBUG_TRACES
+      the_list.push_back(fct.first->vertex(i));
+    }
+  }
+  return the_list;
+}
+
+template <class Vertex_list, class Edge_3>
+Vertex_list from_edge(const Edge_3& edg) {
+  Vertex_list the_list;
+  for (auto i : {edg.second, edg.third}) {
+#ifdef DEBUG_TRACES
+    std::cout << "from edge[" << i << "]=" << edg.first->vertex(i)->point() << std::endl;
+#endif  // DEBUG_TRACES
+    the_list.push_back(edg.first->vertex(i));
+  }
+  return the_list;
+}
+
+template <class Vertex_list, class Vertex_handle>
+Vertex_list from_vertex(const Vertex_handle& vh) {
+  Vertex_list the_list;
+#ifdef DEBUG_TRACES
+  std::cout << "from vertex=" << vh->point() << std::endl;
+#endif  // DEBUG_TRACES
+  the_list.push_back(vh);
+  return the_list;
+}
+
+#endif  // ALPHA_COMPLEX_3D_HELPER_H_
diff --git a/utilities/Alpha_complex/alpha_complex_3d_persistence.cpp b/utilities/Alpha_complex/alpha_complex_3d_persistence.cpp
new file mode 100644
index 00000000..8ef5ffb2
--- /dev/null
+++ b/utilities/Alpha_complex/alpha_complex_3d_persistence.cpp
@@ -0,0 +1,269 @@
+/*    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
+ *
+ *    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 <boost/version.hpp>
+#include <boost/program_options.hpp>
+#include <boost/variant.hpp>
+
+#if BOOST_VERSION >= 105400
+#include <boost/container/static_vector.hpp>
+#endif
+
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Points_3D_off_io.h>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Delaunay_triangulation_3.h>
+#include <CGAL/Alpha_shape_3.h>
+#include <CGAL/iterator.h>
+
+#include <fstream>
+#include <cmath>
+#include <string>
+#include <tuple>
+#include <map>
+#include <utility>
+#include <vector>
+#include <cstdlib>
+
+#include "alpha_complex_3d_helper.h"
+
+// Alpha_shape_3 templates type definitions
+using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel>;
+using Fb = CGAL::Alpha_shape_cell_base_3<Kernel>;
+using Tds = CGAL::Triangulation_data_structure_3<Vb, Fb>;
+using Triangulation_3 = CGAL::Delaunay_triangulation_3<Kernel, Tds>;
+using Alpha_shape_3 = CGAL::Alpha_shape_3<Triangulation_3>;
+
+// From file type definition
+using Point_3 = Kernel::Point_3;
+
+// filtration with alpha values needed type definition
+using Alpha_value_type = Alpha_shape_3::FT;
+using Object = CGAL::Object;
+using Dispatch =
+    CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<Object, Alpha_value_type>,
+                                   CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<Object> >,
+                                                      std::back_insert_iterator<std::vector<Alpha_value_type> > > >;
+using Cell_handle = Alpha_shape_3::Cell_handle;
+using Facet = Alpha_shape_3::Facet;
+using Edge_3 = Alpha_shape_3::Edge;
+using Vertex_handle = Alpha_shape_3::Vertex_handle;
+
+#if BOOST_VERSION >= 105400
+using Vertex_list = boost::container::static_vector<Alpha_shape_3::Vertex_handle, 4>;
+#else
+using Vertex_list = std::vector<Alpha_shape_3::Vertex_handle>;
+#endif
+
+// gudhi type definition
+using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
+using Filtration_value = ST::Filtration_value;
+using Simplex_tree_vertex = ST::Vertex_handle;
+using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>;
+using Simplex_tree_vector_vertex = std::vector<Simplex_tree_vertex>;
+using Persistent_cohomology =
+    Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>;
+
+void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag,
+                     int &coeff_field_characteristic, Filtration_value &min_persistence);
+
+int main(int argc, char **argv) {
+  std::string off_file_points;
+  std::string output_file_diag;
+  int coeff_field_characteristic;
+  Filtration_value min_persistence;
+
+  program_options(argc, argv, off_file_points, output_file_diag, coeff_field_characteristic, min_persistence);
+
+  // Read the OFF file (input file name given as parameter) and triangulate points
+  Gudhi::Points_3D_off_reader<Point_3> off_reader(off_file_points);
+  // Check the read operation was correct
+  if (!off_reader.is_valid()) {
+    std::cerr << "Unable to read file " << off_file_points << std::endl;
+    exit(-1);
+  }
+
+  // Retrieve the points
+  std::vector<Point_3> lp = off_reader.get_point_cloud();
+
+  // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode.
+  Alpha_shape_3 as(lp.begin(), lp.end(), 0, Alpha_shape_3::GENERAL);
+#ifdef DEBUG_TRACES
+  std::cout << "Alpha shape computed in GENERAL mode" << std::endl;
+#endif  // DEBUG_TRACES
+
+  // filtration with alpha values from alpha shape
+  std::vector<Object> the_objects;
+  std::vector<Alpha_value_type> the_alpha_values;
+
+  Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
+                                                                  std::back_inserter(the_alpha_values));
+
+  as.filtration_with_alpha_values(disp);
+#ifdef DEBUG_TRACES
+  std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
+#endif  // DEBUG_TRACES
+
+  Alpha_shape_3::size_type count_vertices = 0;
+  Alpha_shape_3::size_type count_edges = 0;
+  Alpha_shape_3::size_type count_facets = 0;
+  Alpha_shape_3::size_type count_cells = 0;
+
+  // Loop on objects vector
+  Vertex_list vertex_list;
+  ST simplex_tree;
+  Alpha_shape_simplex_tree_map map_cgal_simplex_tree;
+  std::vector<Alpha_value_type>::iterator the_alpha_value_iterator = the_alpha_values.begin();
+  for (auto object_iterator : the_objects) {
+    // Retrieve Alpha shape vertex list from object
+    if (const Cell_handle *cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
+      vertex_list = from_cell<Vertex_list, Cell_handle>(*cell);
+      count_cells++;
+    } else if (const Facet *facet = CGAL::object_cast<Facet>(&object_iterator)) {
+      vertex_list = from_facet<Vertex_list, Facet>(*facet);
+      count_facets++;
+    } else if (const Edge_3 *edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
+      vertex_list = from_edge<Vertex_list, Edge_3>(*edge);
+      count_edges++;
+    } else if (const Vertex_handle *vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) {
+      count_vertices++;
+      vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex);
+    }
+    // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex
+    Simplex_tree_vector_vertex the_simplex;
+    for (auto the_alpha_shape_vertex : vertex_list) {
+      Alpha_shape_simplex_tree_map::iterator the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex);
+      if (the_map_iterator == map_cgal_simplex_tree.end()) {
+        // alpha shape not found
+        Simplex_tree_vertex vertex = map_cgal_simplex_tree.size();
+#ifdef DEBUG_TRACES
+        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl;
+#endif  // DEBUG_TRACES
+        the_simplex.push_back(vertex);
+        map_cgal_simplex_tree.emplace(the_alpha_shape_vertex, vertex);
+      } else {
+        // alpha shape found
+        Simplex_tree_vertex vertex = the_map_iterator->second;
+#ifdef DEBUG_TRACES
+        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl;
+#endif  // DEBUG_TRACES
+        the_simplex.push_back(vertex);
+      }
+    }
+    // Construction of the simplex_tree
+    Filtration_value filtr = /*std::sqrt*/ (*the_alpha_value_iterator);
+#ifdef DEBUG_TRACES
+    std::cout << "filtration = " << filtr << std::endl;
+#endif  // DEBUG_TRACES
+    simplex_tree.insert_simplex(the_simplex, filtr);
+    GUDHI_CHECK(the_alpha_value_iterator != the_alpha_values.end(), "CGAL provided more simplices than values");
+    ++the_alpha_value_iterator;
+  }
+
+#ifdef DEBUG_TRACES
+  std::cout << "vertices \t\t" << count_vertices << std::endl;
+  std::cout << "edges \t\t" << count_edges << std::endl;
+  std::cout << "facets \t\t" << count_facets << std::endl;
+  std::cout << "cells \t\t" << count_cells << std::endl;
+
+  std::cout << "Information of the Simplex Tree: " << std::endl;
+  std::cout << "  Number of vertices = " << simplex_tree.num_vertices() << " ";
+  std::cout << "  Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl;
+  std::cout << "  Dimension = " << simplex_tree.dimension() << " ";
+#endif  // DEBUG_TRACES
+
+#ifdef DEBUG_TRACES
+  std::cout << "Iterator on vertices: " << std::endl;
+  for (auto vertex : simplex_tree.complex_vertex_range()) {
+    std::cout << vertex << " ";
+  }
+#endif  // DEBUG_TRACES
+
+  // Sort the simplices in the order of the filtration
+  simplex_tree.initialize_filtration();
+
+  std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl;
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(simplex_tree, true);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(coeff_field_characteristic);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  // Output the diagram in filediag
+  if (output_file_diag.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::cout << "Result in file: " << output_file_diag << std::endl;
+    std::ofstream out(output_file_diag);
+    pcoh.output_diagram(out);
+    out.close();
+  }
+
+  return 0;
+}
+
+void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag,
+                     int &coeff_field_characteristic, Filtration_value &min_persistence) {
+  namespace po = boost::program_options;
+  po::options_description hidden("Hidden options");
+  hidden.add_options()("input-file", po::value<std::string>(&off_file_points),
+                       "Name of file containing a point set. Format is one point per line:   X1 ... Xd ");
+
+  po::options_description visible("Allowed options", 100);
+  visible.add_options()("help,h", "produce help message")(
+      "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()),
+      "Name of file in which the persistence diagram is written. Default print in std::cout")(
+      "field-charac,p", po::value<int>(&coeff_field_characteristic)->default_value(11),
+      "Characteristic p of the coefficient field Z/pZ for computing homology.")(
+      "min-persistence,m", po::value<Filtration_value>(&min_persistence),
+      "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length "
+      "intervals");
+
+  po::positional_options_description pos;
+  pos.add("input-file", 1);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+
+  po::variables_map vm;
+  po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  if (vm.count("help") || !vm.count("input-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::cout << "of a 3D Alpha complex defined on a set of input points.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
+  }
+}
diff --git a/utilities/Alpha_complex/alpha_complex_persistence.cpp b/utilities/Alpha_complex/alpha_complex_persistence.cpp
new file mode 100644
index 00000000..2105220a
--- /dev/null
+++ b/utilities/Alpha_complex/alpha_complex_persistence.cpp
@@ -0,0 +1,116 @@
+#include <boost/program_options.hpp>
+
+#include <CGAL/Epick_d.h>
+
+#include <gudhi/Alpha_complex.h>
+#include <gudhi/Persistent_cohomology.h>
+// to construct a simplex_tree from alpha complex
+#include <gudhi/Simplex_tree.h>
+
+#include <iostream>
+#include <string>
+#include <limits>  // for numeric_limits
+
+using Simplex_tree = Gudhi::Simplex_tree<>;
+using Filtration_value = Simplex_tree::Filtration_value;
+
+void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag,
+                     Filtration_value &alpha_square_max_value, int &coeff_field_characteristic,
+                     Filtration_value &min_persistence);
+
+int main(int argc, char **argv) {
+  std::string off_file_points;
+  std::string output_file_diag;
+  Filtration_value alpha_square_max_value;
+  int coeff_field_characteristic;
+  Filtration_value min_persistence;
+
+  program_options(argc, argv, off_file_points, output_file_diag, alpha_square_max_value, coeff_field_characteristic,
+                  min_persistence);
+
+  // ----------------------------------------------------------------------------
+  // Init of an alpha complex from an OFF file
+  // ----------------------------------------------------------------------------
+  using Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
+  Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_points);
+
+  Simplex_tree simplex;
+  if (alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) {
+    // ----------------------------------------------------------------------------
+    // Display information about the alpha complex
+    // ----------------------------------------------------------------------------
+    std::cout << "Simplicial complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices()
+              << " simplices - " << simplex.num_vertices() << " vertices." << std::endl;
+
+    // Sort the simplices in the order of the filtration
+    simplex.initialize_filtration();
+
+    std::cout << "Simplex_tree dim: " << simplex.dimension() << std::endl;
+    // Compute the persistence diagram of the complex
+    Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(
+        simplex);
+    // initializes the coefficient field for homology
+    pcoh.init_coefficients(coeff_field_characteristic);
+
+    pcoh.compute_persistent_cohomology(min_persistence);
+
+    // Output the diagram in filediag
+    if (output_file_diag.empty()) {
+      pcoh.output_diagram();
+    } else {
+      std::cout << "Result in file: " << output_file_diag << std::endl;
+      std::ofstream out(output_file_diag);
+      pcoh.output_diagram(out);
+      out.close();
+    }
+  }
+
+  return 0;
+}
+
+void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag,
+                     Filtration_value &alpha_square_max_value, int &coeff_field_characteristic,
+                     Filtration_value &min_persistence) {
+  namespace po = boost::program_options;
+  po::options_description hidden("Hidden options");
+  hidden.add_options()("input-file", po::value<std::string>(&off_file_points),
+                       "Name of file containing a point set. Format is one point per line:   X1 ... Xd ");
+
+  po::options_description visible("Allowed options", 100);
+  visible.add_options()("help,h", "produce help message")(
+      "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()),
+      "Name of file in which the persistence diagram is written. Default print in std::cout")(
+      "max-alpha-square-value,r", po::value<Filtration_value>(&alpha_square_max_value)
+                                      ->default_value(std::numeric_limits<Filtration_value>::infinity()),
+      "Maximal alpha square value for the Alpha complex construction.")(
+      "field-charac,p", po::value<int>(&coeff_field_characteristic)->default_value(11),
+      "Characteristic p of the coefficient field Z/pZ for computing homology.")(
+      "min-persistence,m", po::value<Filtration_value>(&min_persistence),
+      "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length "
+      "intervals");
+
+  po::positional_options_description pos;
+  pos.add("input-file", 1);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+
+  po::variables_map vm;
+  po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  if (vm.count("help") || !vm.count("input-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::cout << "of an Alpha complex defined on a set of input points.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
+  }
+}
diff --git a/utilities/Alpha_complex/alphacomplex.md b/utilities/Alpha_complex/alphacomplex.md
new file mode 100644
index 00000000..aace85d3
--- /dev/null
+++ b/utilities/Alpha_complex/alphacomplex.md
@@ -0,0 +1,158 @@
+

+

+# Alpha complex #

+

+

+## alpha_complex_persistence ##

+

+This program computes the persistent homology with coefficient field Z/pZ of the dD alpha complex built from a dD point cloud.

+The output diagram contains one bar per line, written with the convention:

+

+```

+   p dim birth death

+```

+

+where `dim` is the dimension of the homological feature, `birth` and `death` are respectively the birth and death of the feature,

+and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p` must be a prime number).

+

+**Usage**

+

+```

+   alpha_complex_persistence [options] <input OFF file>

+```

+

+where

+`<input OFF file>` is the path to the input point cloud in [nOFF ASCII format](http://www.geomview.org/docs/html/OFF.html).

+

+**Allowed options**

+

+* `-h [ --help ]` Produce help message

+* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output.

+* `-r [ --max-alpha-square-value ]` (default = inf) Maximal alpha square value for the Alpha complex construction.

+* `-p [ --field-charac ]` (default = 11)     Characteristic p of the coefficient field Z/pZ for computing homology.

+* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals.

+

+**Example**

+

+```

+   alpha_complex_persistence -r 32 -p 2 -m 0.45 ../../data/points/tore3D_300.off

+```

+

+N.B.:

+

+* Filtration values are alpha square values.

+

+

+## alpha_complex_3d_persistence ##

+This program computes the persistent homology with coefficient field Z/pZ of the 3D alpha complex built from a 3D point cloud. The output diagram contains one bar per line, written with the convention:

+

+```

+p dim birth death

+```

+

+where `dim` is the dimension of the homological feature, `birth` and `death` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p` must be a prime number).

+

+**Usage**

+

+```

+   alpha_complex_3d_persistence [options] <input OFF file>

+```

+

+where `<input OFF file>` is the path to the input point cloud in [nOFF ASCII format](http://www.geomview.org/docs/html/OFF.html).

+

+**Allowed options**

+

+* `-h [ --help ]` Produce help message

+* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output.

+* `-p [ --field-charac ]` (default=11) Characteristic p of the coefficient field Z/pZ for computing homology.

+* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals.

+

+**Example**

+

+```

+alpha_complex_3d_persistence ../../data/points/tore3D_300.off -p 2 -m 0.45

+```

+

+N.B.:

+

+* `alpha_complex_3d_persistence` only accepts OFF files in dimension 3.

+* Filtration values are alpha square values.

+

+

+## exact_alpha_complex_3d_persistence ##

+

+Same as `alpha_complex_3d_persistence`, but using exact computation.

+It is slower, but it is necessary when points are on a grid for instance.

+

+

+

+## weighted_alpha_complex_3d_persistence ##

+

+Same as `alpha_complex_3d_persistence`, but using weighted points.

+

+**Usage**

+

+```

+   weighted_alpha_complex_3d_persistence [options] <input OFF file> <weights input file>

+```

+

+where

+

+* `<input OFF file>` is the path to the input point cloud in [nOFF ASCII format](http://www.geomview.org/docs/html/OFF.html).

+* `<input weights file>` is the path to the file containing the weights of the points (one value per line).

+

+**Allowed options**

+

+* `-h [ --help ]` Produce help message

+* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output.

+* `-p [ --field-charac ]` (default=11) Characteristic p of the coefficient field Z/pZ for computing homology.

+* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals.

+

+**Example**

+

+```

+   weighted_alpha_complex_3d_persistence ../../data/points/tore3D_300.off ../../data/points/tore3D_300.weights -p 2 -m 0.45

+```

+

+

+N.B.:

+

+* Weights values are explained on CGAL [Alpha shape](https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title0)

+and [Regular triangulation](https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation) documentation.

+* Filtration values are alpha square values.

+

+

+## periodic_alpha_complex_3d_persistence ##

+Same as `alpha_complex_3d_persistence`, but using periodic alpha shape 3d.

+Refer to the [CGAL's 3D Periodic Triangulations User Manual](https://doc.cgal.org/latest/Periodic_3_triangulation_3/index.html) for more details.

+

+**Usage**

+

+```

+   periodic_alpha_complex_3d_persistence [options] <input OFF file> <cuboid file>

+```

+

+where

+

+* `<input OFF file>` is the path to the input point cloud in [nOFF ASCII format](http://www.geomview.org/docs/html/OFF.html).

+* `<cuboid file>` is the path to the file describing the periodic domain. It must be in the format described

+[here](/doc/latest/fileformats.html#FileFormatsIsoCuboid).

+

+**Allowed options**

+

+* `-h [ --help ]` Produce help message

+* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output.

+* `-p [ --field-charac ]` (default=11) Characteristic p of the coefficient field Z/pZ for computing homology.

+* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals

+

+

+**Example**

+

+```

+periodic_alpha_complex_3d_persistence ../../data/points/grid_10_10_10_in_0_1.off ../../data/points/iso_cuboid_3_in_0_1.txt -p 3 -m 1.0

+```

+

+N.B.:

+

+* Cuboid file must be in the format described [here](/doc/latest/fileformats.html#FileFormatsIsoCuboid).

+* Filtration values are alpha square values.

diff --git a/utilities/Alpha_complex/exact_alpha_complex_3d_persistence.cpp b/utilities/Alpha_complex/exact_alpha_complex_3d_persistence.cpp
new file mode 100644
index 00000000..cceac46e
--- /dev/null
+++ b/utilities/Alpha_complex/exact_alpha_complex_3d_persistence.cpp
@@ -0,0 +1,263 @@
+/*    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
+ *
+ *    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 <boost/program_options.hpp>
+#include <boost/variant.hpp>
+
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Points_3D_off_io.h>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Delaunay_triangulation_3.h>
+#include <CGAL/Alpha_shape_3.h>
+#include <CGAL/iterator.h>
+
+#include <fstream>
+#include <cmath>
+#include <string>
+#include <tuple>
+#include <map>
+#include <utility>
+#include <vector>
+#include <cstdlib>
+
+#include "alpha_complex_3d_helper.h"
+
+// Alpha_shape_3 templates type definitions
+using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Exact_tag = CGAL::Tag_true;
+using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel, CGAL::Default, Exact_tag>;
+using Fb = CGAL::Alpha_shape_cell_base_3<Kernel, CGAL::Default, Exact_tag>;
+using Tds = CGAL::Triangulation_data_structure_3<Vb, Fb>;
+using Triangulation_3 = CGAL::Delaunay_triangulation_3<Kernel, Tds>;
+using Alpha_shape_3 = CGAL::Alpha_shape_3<Triangulation_3, Exact_tag>;
+
+// From file type definition
+using Point_3 = Kernel::Point_3;
+
+// filtration with alpha values needed type definition
+using Alpha_value_type = Alpha_shape_3::FT;
+using Object = CGAL::Object;
+using Dispatch =
+    CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<Object, Alpha_value_type>,
+                                   CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<Object> >,
+                                                      std::back_insert_iterator<std::vector<Alpha_value_type> > > >;
+using Cell_handle = Alpha_shape_3::Cell_handle;
+using Facet = Alpha_shape_3::Facet;
+using Edge_3 = Alpha_shape_3::Edge;
+using Vertex_handle = Alpha_shape_3::Vertex_handle;
+using Vertex_list = std::vector<Vertex_handle>;
+
+// gudhi type definition
+using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
+using Filtration_value = ST::Filtration_value;
+using Simplex_tree_vertex = ST::Vertex_handle;
+using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>;
+using Simplex_tree_vector_vertex = std::vector<Simplex_tree_vertex>;
+using Persistent_cohomology =
+    Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>;
+
+void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag,
+                     int &coeff_field_characteristic, Filtration_value &min_persistence);
+
+int main(int argc, char **argv) {
+  std::string off_file_points;
+  std::string output_file_diag;
+  int coeff_field_characteristic;
+  Filtration_value min_persistence;
+
+  program_options(argc, argv, off_file_points, output_file_diag, coeff_field_characteristic, min_persistence);
+
+  // Read the OFF file (input file name given as parameter) and triangulate points
+  Gudhi::Points_3D_off_reader<Point_3> off_reader(off_file_points);
+  // Check the read operation was correct
+  if (!off_reader.is_valid()) {
+    std::cerr << "Unable to read file " << off_file_points << std::endl;
+    exit(-1);
+  }
+
+  // Retrieve the points
+  std::vector<Point_3> lp = off_reader.get_point_cloud();
+
+  // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode.
+  Alpha_shape_3 as(lp.begin(), lp.end(), 0, Alpha_shape_3::GENERAL);
+#ifdef DEBUG_TRACES
+  std::cout << "Alpha shape computed in GENERAL mode" << std::endl;
+#endif  // DEBUG_TRACES
+
+  // filtration with alpha values from alpha shape
+  std::vector<Object> the_objects;
+  std::vector<Alpha_value_type> the_alpha_values;
+
+  Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
+                                                                  std::back_inserter(the_alpha_values));
+
+  as.filtration_with_alpha_values(disp);
+#ifdef DEBUG_TRACES
+  std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
+#endif  // DEBUG_TRACES
+
+  Alpha_shape_3::size_type count_vertices = 0;
+  Alpha_shape_3::size_type count_edges = 0;
+  Alpha_shape_3::size_type count_facets = 0;
+  Alpha_shape_3::size_type count_cells = 0;
+
+  // Loop on objects vector
+  Vertex_list vertex_list;
+  ST simplex_tree;
+  Alpha_shape_simplex_tree_map map_cgal_simplex_tree;
+  std::vector<Alpha_value_type>::iterator the_alpha_value_iterator = the_alpha_values.begin();
+  for (auto object_iterator : the_objects) {
+    // Retrieve Alpha shape vertex list from object
+    if (const Cell_handle *cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
+      vertex_list = from_cell<Vertex_list, Cell_handle>(*cell);
+      count_cells++;
+    } else if (const Facet *facet = CGAL::object_cast<Facet>(&object_iterator)) {
+      vertex_list = from_facet<Vertex_list, Facet>(*facet);
+      count_facets++;
+    } else if (const Edge_3 *edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
+      vertex_list = from_edge<Vertex_list, Edge_3>(*edge);
+      count_edges++;
+    } else if (const Vertex_handle *vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) {
+      count_vertices++;
+      vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex);
+    }
+    // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex
+    Simplex_tree_vector_vertex the_simplex;
+    for (auto the_alpha_shape_vertex : vertex_list) {
+      Alpha_shape_simplex_tree_map::iterator the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex);
+      if (the_map_iterator == map_cgal_simplex_tree.end()) {
+        // alpha shape not found
+        Simplex_tree_vertex vertex = map_cgal_simplex_tree.size();
+#ifdef DEBUG_TRACES
+        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl;
+#endif  // DEBUG_TRACES
+        the_simplex.push_back(vertex);
+        map_cgal_simplex_tree.emplace(the_alpha_shape_vertex, vertex);
+      } else {
+        // alpha shape found
+        Simplex_tree_vertex vertex = the_map_iterator->second;
+#ifdef DEBUG_TRACES
+        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl;
+#endif  // DEBUG_TRACES
+        the_simplex.push_back(vertex);
+      }
+    }
+    // Construction of the simplex_tree
+    // you can also use the_alpha_value_iterator->exact()
+    Filtration_value filtr = /*std::sqrt*/ CGAL::to_double(the_alpha_value_iterator->exact());
+#ifdef DEBUG_TRACES
+    std::cout << "filtration = " << filtr << std::endl;
+#endif  // DEBUG_TRACES
+    simplex_tree.insert_simplex(the_simplex, filtr);
+    if (the_alpha_value_iterator != the_alpha_values.end())
+      ++the_alpha_value_iterator;
+    else
+      std::cout << "This shall not happen" << std::endl;
+  }
+
+#ifdef DEBUG_TRACES
+  std::cout << "vertices \t\t" << count_vertices << std::endl;
+  std::cout << "edges \t\t" << count_edges << std::endl;
+  std::cout << "facets \t\t" << count_facets << std::endl;
+  std::cout << "cells \t\t" << count_cells << std::endl;
+
+  std::cout << "Information of the Simplex Tree: " << std::endl;
+  std::cout << "  Number of vertices = " << simplex_tree.num_vertices() << " ";
+  std::cout << "  Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl;
+  std::cout << "  Dimension = " << simplex_tree.dimension() << " ";
+#endif  // DEBUG_TRACES
+
+#ifdef DEBUG_TRACES
+  std::cout << "Iterator on vertices: " << std::endl;
+  for (auto vertex : simplex_tree.complex_vertex_range()) {
+    std::cout << vertex << " ";
+  }
+#endif  // DEBUG_TRACES
+
+  // Sort the simplices in the order of the filtration
+  simplex_tree.initialize_filtration();
+
+  std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl;
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(simplex_tree, true);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(coeff_field_characteristic);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  // Output the diagram in filediag
+  if (output_file_diag.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::cout << "Result in file: " << output_file_diag << std::endl;
+    std::ofstream out(output_file_diag);
+    pcoh.output_diagram(out);
+    out.close();
+  }
+
+  return 0;
+}
+
+void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag,
+                     int &coeff_field_characteristic, Filtration_value &min_persistence) {
+  namespace po = boost::program_options;
+  po::options_description hidden("Hidden options");
+  hidden.add_options()("input-file", po::value<std::string>(&off_file_points),
+                       "Name of file containing a point set. Format is one point per line:   X1 ... Xd ");
+
+  po::options_description visible("Allowed options", 100);
+  visible.add_options()("help,h", "produce help message")(
+      "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()),
+      "Name of file in which the persistence diagram is written. Default print in std::cout")(
+      "field-charac,p", po::value<int>(&coeff_field_characteristic)->default_value(11),
+      "Characteristic p of the coefficient field Z/pZ for computing homology.")(
+      "min-persistence,m", po::value<Filtration_value>(&min_persistence),
+      "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length "
+      "intervals");
+
+  po::positional_options_description pos;
+  pos.add("input-file", 1);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+
+  po::variables_map vm;
+  po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  if (vm.count("help") || !vm.count("input-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::cout << "of a 3D Alpha complex defined on a set of input points.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
+  }
+}
diff --git a/utilities/Alpha_complex/periodic_alpha_complex_3d_persistence.cpp b/utilities/Alpha_complex/periodic_alpha_complex_3d_persistence.cpp
new file mode 100644
index 00000000..188cf604
--- /dev/null
+++ b/utilities/Alpha_complex/periodic_alpha_complex_3d_persistence.cpp
@@ -0,0 +1,300 @@
+/*    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
+ *                     Pawel Dlotko - 2017 - Swansea University, UK
+ *
+ *    Copyright (C) 2014  INRIA
+ *
+ *    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 <boost/program_options.hpp>
+#include <boost/variant.hpp>
+
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Points_3D_off_io.h>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
+#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
+#include <CGAL/Alpha_shape_3.h>
+#include <CGAL/iterator.h>
+
+#include <fstream>
+#include <cmath>
+#include <string>
+#include <tuple>
+#include <map>
+#include <utility>
+#include <vector>
+#include <cstdlib>
+
+#include "alpha_complex_3d_helper.h"
+
+// Traits
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using PK = CGAL::Periodic_3_Delaunay_triangulation_traits_3<K>;
+// Vertex type
+using DsVb = CGAL::Periodic_3_triangulation_ds_vertex_base_3<>;
+using Vb = CGAL::Triangulation_vertex_base_3<PK, DsVb>;
+using AsVb = CGAL::Alpha_shape_vertex_base_3<PK, Vb>;
+// Cell type
+using DsCb = CGAL::Periodic_3_triangulation_ds_cell_base_3<>;
+using Cb = CGAL::Triangulation_cell_base_3<PK, DsCb>;
+using AsCb = CGAL::Alpha_shape_cell_base_3<PK, Cb>;
+using Tds = CGAL::Triangulation_data_structure_3<AsVb, AsCb>;
+using P3DT3 = CGAL::Periodic_3_Delaunay_triangulation_3<PK, Tds>;
+using Alpha_shape_3 = CGAL::Alpha_shape_3<P3DT3>;
+using Point_3 = PK::Point_3;
+
+// filtration with alpha values needed type definition
+using Alpha_value_type = Alpha_shape_3::FT;
+using Object = CGAL::Object;
+using Dispatch =
+    CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<Object, Alpha_value_type>,
+                                   CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<Object> >,
+                                                      std::back_insert_iterator<std::vector<Alpha_value_type> > > >;
+using Cell_handle = Alpha_shape_3::Cell_handle;
+using Facet = Alpha_shape_3::Facet;
+using Edge_3 = Alpha_shape_3::Edge;
+using Vertex_handle = Alpha_shape_3::Vertex_handle;
+using Vertex_list = std::vector<Alpha_shape_3::Vertex_handle>;
+
+// gudhi type definition
+using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
+using Filtration_value = ST::Filtration_value;
+using Simplex_tree_vertex = ST::Vertex_handle;
+using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>;
+using Simplex_tree_vector_vertex = std::vector<Simplex_tree_vertex>;
+using Persistent_cohomology =
+    Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>;
+
+void program_options(int argc, char *argv[], std::string &off_file_points, std::string &cuboid_file,
+                     std::string &output_file_diag, int &coeff_field_characteristic, Filtration_value &min_persistence);
+
+int main(int argc, char **argv) {
+  std::string off_file_points;
+  std::string cuboid_file;
+  std::string output_file_diag;
+  int coeff_field_characteristic;
+  Filtration_value min_persistence;
+
+  program_options(argc, argv, off_file_points, cuboid_file, output_file_diag, coeff_field_characteristic,
+                  min_persistence);
+
+  // Read the OFF file (input file name given as parameter) and triangulate points
+  Gudhi::Points_3D_off_reader<Point_3> off_reader(off_file_points);
+  // Check the read operation was correct
+  if (!off_reader.is_valid()) {
+    std::cerr << "Unable to read OFF file " << off_file_points << std::endl;
+    exit(-1);
+  }
+
+  // Read iso_cuboid_3 information from file
+  std::ifstream iso_cuboid_str(cuboid_file);
+  double x_min, y_min, z_min, x_max, y_max, z_max;
+  if (iso_cuboid_str.good()) {
+    iso_cuboid_str >> x_min >> y_min >> z_min >> x_max >> y_max >> z_max;
+  } else {
+    std::cerr << "Unable to read file " << cuboid_file << std::endl;
+    exit(-1);
+  }
+  // Checking if the cuboid is the same in x,y and z direction. If not, CGAL will not process it.
+  if ((x_max - x_min != y_max - y_min) || (x_max - x_min != z_max - z_min) || (z_max - z_min != y_max - y_min)) {
+    std::cerr << "The size of the cuboid in every directions is not the same." << std::endl;
+    exit(-1);
+  }
+
+  // Retrieve the points
+  std::vector<Point_3> lp = off_reader.get_point_cloud();
+
+  // Define the periodic cube
+  P3DT3 pdt(PK::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max));
+  // Heuristic for inserting large point sets (if pts is reasonably large)
+  pdt.insert(lp.begin(), lp.end(), true);
+  // As pdt won't be modified anymore switch to 1-sheeted cover if possible
+  if (pdt.is_triangulation_in_1_sheet()) {
+    pdt.convert_to_1_sheeted_covering();
+  } else {
+    std::cerr << "ERROR: we were not able to construct a triangulation within a single periodic domain." << std::endl;
+    exit(-1);
+  }
+  std::cout << "Periodic Delaunay computed." << std::endl;
+
+  // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. This is the default mode
+  // Maybe need to set it to GENERAL mode
+  Alpha_shape_3 as(pdt, 0, Alpha_shape_3::GENERAL);
+
+  // filtration with alpha values from alpha shape
+  std::vector<Object> the_objects;
+  std::vector<Alpha_value_type> the_alpha_values;
+
+  Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
+                                                                  std::back_inserter(the_alpha_values));
+
+  as.filtration_with_alpha_values(disp);
+#ifdef DEBUG_TRACES
+  std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
+#endif  // DEBUG_TRACES
+
+  Alpha_shape_3::size_type count_vertices = 0;
+  Alpha_shape_3::size_type count_edges = 0;
+  Alpha_shape_3::size_type count_facets = 0;
+  Alpha_shape_3::size_type count_cells = 0;
+
+  // Loop on objects vector
+  Vertex_list vertex_list;
+  ST simplex_tree;
+  Alpha_shape_simplex_tree_map map_cgal_simplex_tree;
+  std::vector<Alpha_value_type>::iterator the_alpha_value_iterator = the_alpha_values.begin();
+  for (auto object_iterator : the_objects) {
+    // Retrieve Alpha shape vertex list from object
+    if (const Cell_handle *cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
+      vertex_list = from_cell<Vertex_list, Cell_handle>(*cell);
+      count_cells++;
+    } else if (const Facet *facet = CGAL::object_cast<Facet>(&object_iterator)) {
+      vertex_list = from_facet<Vertex_list, Facet>(*facet);
+      count_facets++;
+    } else if (const Edge_3 *edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
+      vertex_list = from_edge<Vertex_list, Edge_3>(*edge);
+      count_edges++;
+    } else if (const Vertex_handle *vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) {
+      count_vertices++;
+      vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex);
+    }
+    // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex
+    Simplex_tree_vector_vertex the_simplex;
+    for (auto the_alpha_shape_vertex : vertex_list) {
+      Alpha_shape_simplex_tree_map::iterator the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex);
+      if (the_map_iterator == map_cgal_simplex_tree.end()) {
+        // alpha shape not found
+        Simplex_tree_vertex vertex = map_cgal_simplex_tree.size();
+#ifdef DEBUG_TRACES
+        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl;
+#endif  // DEBUG_TRACES
+        the_simplex.push_back(vertex);
+        map_cgal_simplex_tree.emplace(the_alpha_shape_vertex, vertex);
+      } else {
+        // alpha shape found
+        Simplex_tree_vertex vertex = the_map_iterator->second;
+#ifdef DEBUG_TRACES
+        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl;
+#endif  // DEBUG_TRACES
+        the_simplex.push_back(vertex);
+      }
+    }
+    // Construction of the simplex_tree
+    Filtration_value filtr = /*std::sqrt*/ (*the_alpha_value_iterator);
+#ifdef DEBUG_TRACES
+    std::cout << "filtration = " << filtr << std::endl;
+#endif  // DEBUG_TRACES
+    simplex_tree.insert_simplex(the_simplex, filtr);
+    if (the_alpha_value_iterator != the_alpha_values.end())
+      ++the_alpha_value_iterator;
+    else
+      std::cout << "This shall not happen" << std::endl;
+  }
+
+#ifdef DEBUG_TRACES
+  std::cout << "vertices \t\t" << count_vertices << std::endl;
+  std::cout << "edges \t\t" << count_edges << std::endl;
+  std::cout << "facets \t\t" << count_facets << std::endl;
+  std::cout << "cells \t\t" << count_cells << std::endl;
+
+  std::cout << "Information of the Simplex Tree: " << std::endl;
+  std::cout << "  Number of vertices = " << simplex_tree.num_vertices() << " ";
+  std::cout << "  Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl;
+  std::cout << "  Dimension = " << simplex_tree.dimension() << " ";
+#endif  // DEBUG_TRACES
+
+#ifdef DEBUG_TRACES
+  std::cout << "Iterator on vertices: " << std::endl;
+  for (auto vertex : simplex_tree.complex_vertex_range()) {
+    std::cout << vertex << " ";
+  }
+#endif  // DEBUG_TRACES
+
+  // Sort the simplices in the order of the filtration
+  simplex_tree.initialize_filtration();
+
+  std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl;
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(simplex_tree, true);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(coeff_field_characteristic);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  // Output the diagram in filediag
+  if (output_file_diag.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::cout << "Result in file: " << output_file_diag << std::endl;
+    std::ofstream out(output_file_diag);
+    pcoh.output_diagram(out);
+    out.close();
+  }
+
+  return 0;
+}
+
+void program_options(int argc, char *argv[], std::string &off_file_points, std::string &cuboid_file,
+                     std::string &output_file_diag, int &coeff_field_characteristic,
+                     Filtration_value &min_persistence) {
+  namespace po = boost::program_options;
+  po::options_description hidden("Hidden options");
+  hidden.add_options()("input-file", po::value<std::string>(&off_file_points),
+                       "Name of file containing a point set. Format is one point per line:   X1 ... Xd ")(
+      "cuboid-file", po::value<std::string>(&cuboid_file),
+      "Name of file describing the periodic domain. Format is: min_hx min_hy min_hz\nmax_hx max_hy max_hz");
+
+  po::options_description visible("Allowed options", 100);
+  visible.add_options()("help,h", "produce help message")(
+      "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()),
+      "Name of file in which the persistence diagram is written. Default print in std::cout")(
+      "field-charac,p", po::value<int>(&coeff_field_characteristic)->default_value(11),
+      "Characteristic p of the coefficient field Z/pZ for computing homology.")(
+      "min-persistence,m", po::value<Filtration_value>(&min_persistence),
+      "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length "
+      "intervals");
+
+  po::positional_options_description pos;
+  pos.add("input-file", 1);
+  pos.add("cuboid-file", 2);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+
+  po::variables_map vm;
+  po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  if (vm.count("help") || !vm.count("input-file") || !vm.count("cuboid-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::cout << "of a periodic 3D Alpha complex defined on a set of input points.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file cuboid-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
+  }
+}
diff --git a/utilities/Alpha_complex/weighted_alpha_complex_3d_persistence.cpp b/utilities/Alpha_complex/weighted_alpha_complex_3d_persistence.cpp
new file mode 100644
index 00000000..93be8a05
--- /dev/null
+++ b/utilities/Alpha_complex/weighted_alpha_complex_3d_persistence.cpp
@@ -0,0 +1,314 @@
+/*    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
+ *
+ *    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 <boost/program_options.hpp>
+#include <boost/variant.hpp>
+
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Points_3D_off_io.h>
+
+#include <CGAL/config.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Regular_triangulation_3.h>
+#include <CGAL/Alpha_shape_3.h>
+#include <CGAL/iterator.h>
+
+// For CGAL < 4.11
+#if CGAL_VERSION_NR < 1041100000
+#include <CGAL/Regular_triangulation_euclidean_traits_3.h>
+#endif  // CGAL_VERSION_NR < 1041100000
+
+#include <fstream>
+#include <cmath>
+#include <string>
+#include <tuple>
+#include <map>
+#include <utility>
+#include <vector>
+#include <cstdlib>
+
+#include "alpha_complex_3d_helper.h"
+
+// Alpha_shape_3 templates type definitions
+using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
+
+// For CGAL < 4.11
+#if CGAL_VERSION_NR < 1041100000
+using Gt = CGAL::Regular_triangulation_euclidean_traits_3<Kernel>;
+using Vb = CGAL::Alpha_shape_vertex_base_3<Gt>;
+using Fb = CGAL::Alpha_shape_cell_base_3<Gt>;
+using Tds = CGAL::Triangulation_data_structure_3<Vb, Fb>;
+using Triangulation_3 = CGAL::Regular_triangulation_3<Gt, Tds>;
+
+// From file type definition
+using Point_3 = Gt::Bare_point;
+using Weighted_point_3 = Gt::Weighted_point;
+
+// For CGAL >= 4.11
+#else   // CGAL_VERSION_NR < 1041100000
+using Rvb = CGAL::Regular_triangulation_vertex_base_3<Kernel>;
+using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel, Rvb>;
+using Rcb = CGAL::Regular_triangulation_cell_base_3<Kernel>;
+using Cb = CGAL::Alpha_shape_cell_base_3<Kernel, Rcb>;
+using Tds = CGAL::Triangulation_data_structure_3<Vb, Cb>;
+using Triangulation_3 = CGAL::Regular_triangulation_3<Kernel, Tds>;
+
+// From file type definition
+using Point_3 = Triangulation_3::Bare_point;
+using Weighted_point_3 = Triangulation_3::Weighted_point;
+#endif  // CGAL_VERSION_NR < 1041100000
+
+using Alpha_shape_3 = CGAL::Alpha_shape_3<Triangulation_3>;
+
+// filtration with alpha values needed type definition
+using Alpha_value_type = Alpha_shape_3::FT;
+using Object = CGAL::Object;
+using Dispatch =
+    CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<Object, Alpha_value_type>,
+                                   CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<Object> >,
+                                                      std::back_insert_iterator<std::vector<Alpha_value_type> > > >;
+using Cell_handle = Alpha_shape_3::Cell_handle;
+using Facet = Alpha_shape_3::Facet;
+using Edge_3 = Alpha_shape_3::Edge;
+using Vertex_handle = Alpha_shape_3::Vertex_handle;
+using Vertex_list = std::vector<Alpha_shape_3::Vertex_handle>;
+
+// gudhi type definition
+using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
+using Filtration_value = ST::Filtration_value;
+using Simplex_tree_vertex = ST::Vertex_handle;
+using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>;
+using Simplex_tree_vector_vertex = std::vector<Simplex_tree_vertex>;
+using Persistent_cohomology =
+    Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>;
+
+void program_options(int argc, char *argv[], std::string &off_file_points, std::string &weight_file,
+                     std::string &output_file_diag, int &coeff_field_characteristic, Filtration_value &min_persistence);
+
+int main(int argc, char **argv) {
+  std::string off_file_points;
+  std::string weight_file;
+  std::string output_file_diag;
+  int coeff_field_characteristic;
+  Filtration_value min_persistence;
+
+  program_options(argc, argv, off_file_points, weight_file, output_file_diag, coeff_field_characteristic,
+                  min_persistence);
+
+  // Read the OFF file (input file name given as parameter) and triangulate points
+  Gudhi::Points_3D_off_reader<Point_3> off_reader(off_file_points);
+  // Check the read operation was correct
+  if (!off_reader.is_valid()) {
+    std::cerr << "Unable to read OFF file " << off_file_points << std::endl;
+    exit(-1);
+  }
+
+  // Retrieve the points
+  std::vector<Point_3> lp = off_reader.get_point_cloud();
+
+  // Read weights information from file
+  std::ifstream weights_ifstr(weight_file);
+  std::vector<Weighted_point_3> wp;
+  if (weights_ifstr.good()) {
+    double weight = 0.0;
+    std::size_t index = 0;
+    wp.reserve(lp.size());
+    // Attempt read the weight in a double format, return false if it fails
+    while ((weights_ifstr >> weight) && (index < lp.size())) {
+      wp.push_back(Weighted_point_3(lp[index], weight));
+      index++;
+    }
+    if (index != lp.size()) {
+      std::cerr << "Bad number of weights in file " << weight_file << std::endl;
+      exit(-1);
+    }
+  } else {
+    std::cerr << "Unable to read weights file " << weight_file << std::endl;
+    exit(-1);
+  }
+
+  // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode.
+  Alpha_shape_3 as(wp.begin(), wp.end(), 0, Alpha_shape_3::GENERAL);
+#ifdef DEBUG_TRACES
+  std::cout << "Alpha shape computed in GENERAL mode" << std::endl;
+#endif  // DEBUG_TRACES
+
+  // filtration with alpha values from alpha shape
+  std::vector<Object> the_objects;
+  std::vector<Alpha_value_type> the_alpha_values;
+
+  Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
+                                                                  std::back_inserter(the_alpha_values));
+
+  as.filtration_with_alpha_values(disp);
+#ifdef DEBUG_TRACES
+  std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
+#endif  // DEBUG_TRACES
+
+  Alpha_shape_3::size_type count_vertices = 0;
+  Alpha_shape_3::size_type count_edges = 0;
+  Alpha_shape_3::size_type count_facets = 0;
+  Alpha_shape_3::size_type count_cells = 0;
+
+  // Loop on objects vector
+  Vertex_list vertex_list;
+  ST simplex_tree;
+  Alpha_shape_simplex_tree_map map_cgal_simplex_tree;
+  std::vector<Alpha_value_type>::iterator the_alpha_value_iterator = the_alpha_values.begin();
+  for (auto object_iterator : the_objects) {
+    // Retrieve Alpha shape vertex list from object
+    if (const Cell_handle *cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
+      vertex_list = from_cell<Vertex_list, Cell_handle>(*cell);
+      count_cells++;
+    } else if (const Facet *facet = CGAL::object_cast<Facet>(&object_iterator)) {
+      vertex_list = from_facet<Vertex_list, Facet>(*facet);
+      count_facets++;
+    } else if (const Edge_3 *edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
+      vertex_list = from_edge<Vertex_list, Edge_3>(*edge);
+      count_edges++;
+    } else if (const Vertex_handle *vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) {
+      count_vertices++;
+      vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex);
+    }
+    // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex
+    Simplex_tree_vector_vertex the_simplex;
+    for (auto the_alpha_shape_vertex : vertex_list) {
+      Alpha_shape_simplex_tree_map::iterator the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex);
+      if (the_map_iterator == map_cgal_simplex_tree.end()) {
+        // alpha shape not found
+        Simplex_tree_vertex vertex = map_cgal_simplex_tree.size();
+#ifdef DEBUG_TRACES
+        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl;
+#endif  // DEBUG_TRACES
+        the_simplex.push_back(vertex);
+        map_cgal_simplex_tree.emplace(the_alpha_shape_vertex, vertex);
+      } else {
+        // alpha shape found
+        Simplex_tree_vertex vertex = the_map_iterator->second;
+#ifdef DEBUG_TRACES
+        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl;
+#endif  // DEBUG_TRACES
+        the_simplex.push_back(vertex);
+      }
+    }
+    // Construction of the simplex_tree
+    Filtration_value filtr = /*std::sqrt*/ (*the_alpha_value_iterator);
+#ifdef DEBUG_TRACES
+    std::cout << "filtration = " << filtr << std::endl;
+#endif  // DEBUG_TRACES
+    simplex_tree.insert_simplex(the_simplex, filtr);
+    if (the_alpha_value_iterator != the_alpha_values.end())
+      ++the_alpha_value_iterator;
+    else
+      std::cout << "This shall not happen" << std::endl;
+  }
+
+#ifdef DEBUG_TRACES
+  std::cout << "vertices \t\t" << count_vertices << std::endl;
+  std::cout << "edges \t\t" << count_edges << std::endl;
+  std::cout << "facets \t\t" << count_facets << std::endl;
+  std::cout << "cells \t\t" << count_cells << std::endl;
+
+  std::cout << "Information of the Simplex Tree: " << std::endl;
+  std::cout << "  Number of vertices = " << simplex_tree.num_vertices() << " ";
+  std::cout << "  Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl;
+  std::cout << "  Dimension = " << simplex_tree.dimension() << " ";
+#endif  // DEBUG_TRACES
+
+#ifdef DEBUG_TRACES
+  std::cout << "Iterator on vertices: " << std::endl;
+  for (auto vertex : simplex_tree.complex_vertex_range()) {
+    std::cout << vertex << " ";
+  }
+#endif  // DEBUG_TRACES
+
+  // Sort the simplices in the order of the filtration
+  simplex_tree.initialize_filtration();
+
+  std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl;
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(simplex_tree, true);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(coeff_field_characteristic);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  // Output the diagram in filediag
+  if (output_file_diag.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::cout << "Result in file: " << output_file_diag << std::endl;
+    std::ofstream out(output_file_diag);
+    pcoh.output_diagram(out);
+    out.close();
+  }
+
+  return 0;
+}
+
+void program_options(int argc, char *argv[], std::string &off_file_points, std::string &weight_file,
+                     std::string &output_file_diag, int &coeff_field_characteristic,
+                     Filtration_value &min_persistence) {
+  namespace po = boost::program_options;
+  po::options_description hidden("Hidden options");
+  hidden.add_options()("input-file", po::value<std::string>(&off_file_points),
+                       "Name of file containing a point set. Format is one point per line:   X1 ... Xd ")(
+      "weight-file", po::value<std::string>(&weight_file),
+      "Name of file containing a point weights. Format is one weigt per line:   W1\n...\nWn ");
+
+  po::options_description visible("Allowed options", 100);
+  visible.add_options()("help,h", "produce help message")(
+      "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()),
+      "Name of file in which the persistence diagram is written. Default print in std::cout")(
+      "field-charac,p", po::value<int>(&coeff_field_characteristic)->default_value(11),
+      "Characteristic p of the coefficient field Z/pZ for computing homology.")(
+      "min-persistence,m", po::value<Filtration_value>(&min_persistence),
+      "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length "
+      "intervals");
+
+  po::positional_options_description pos;
+  pos.add("input-file", 1);
+  pos.add("weight-file", 2);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+
+  po::variables_map vm;
+  po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  if (vm.count("help") || !vm.count("input-file") || !vm.count("weight-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::cout << "of a weighted 3D Alpha complex defined on a set of input points.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file weight-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
+  }
+}
diff --git a/utilities/Alpha_complex/weighted_periodic_alpha_complex_3d_persistence.cpp b/utilities/Alpha_complex/weighted_periodic_alpha_complex_3d_persistence.cpp
new file mode 100644
index 00000000..5321bb0a
--- /dev/null
+++ b/utilities/Alpha_complex/weighted_periodic_alpha_complex_3d_persistence.cpp
@@ -0,0 +1,286 @@
+/*    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
+ *                     Pawel Dlotko - 2017 - Swansea University, UK
+ *
+ *    Copyright (C) 2014  INRIA
+ *
+ *    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 <boost/variant.hpp>
+
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Points_3D_off_io.h>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Periodic_3_regular_triangulation_traits_3.h>
+#include <CGAL/Periodic_3_regular_triangulation_3.h>
+#include <CGAL/Alpha_shape_3.h>
+#include <CGAL/iterator.h>
+
+#include <fstream>
+#include <cmath>
+#include <string>
+#include <tuple>
+#include <map>
+#include <utility>
+#include <vector>
+#include <cstdlib>
+
+#include "alpha_complex_3d_helper.h"
+
+// Traits
+using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
+using PK = CGAL::Periodic_3_regular_triangulation_traits_3<Kernel>;
+
+// Vertex type
+using DsVb = CGAL::Periodic_3_triangulation_ds_vertex_base_3<>;
+using Vb = CGAL::Regular_triangulation_vertex_base_3<PK, DsVb>;
+using AsVb = CGAL::Alpha_shape_vertex_base_3<PK, Vb>;
+// Cell type
+using DsCb = CGAL::Periodic_3_triangulation_ds_cell_base_3<>;
+using Cb = CGAL::Regular_triangulation_cell_base_3<PK, DsCb>;
+using AsCb = CGAL::Alpha_shape_cell_base_3<PK, Cb>;
+using Tds = CGAL::Triangulation_data_structure_3<AsVb, AsCb>;
+using P3RT3 = CGAL::Periodic_3_regular_triangulation_3<PK, Tds>;
+using Alpha_shape_3 = CGAL::Alpha_shape_3<P3RT3>;
+
+using Point_3 = P3RT3::Bare_point;
+using Weighted_point_3 = P3RT3::Weighted_point;
+
+// filtration with alpha values needed type definition
+using Alpha_value_type = Alpha_shape_3::FT;
+using Object = CGAL::Object;
+using Dispatch =
+    CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<Object, Alpha_value_type>,
+                                   CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<Object> >,
+                                                      std::back_insert_iterator<std::vector<Alpha_value_type> > > >;
+using Cell_handle = Alpha_shape_3::Cell_handle;
+using Facet = Alpha_shape_3::Facet;
+using Edge_3 = Alpha_shape_3::Edge;
+using Vertex_handle = Alpha_shape_3::Vertex_handle;
+using Vertex_list = std::vector<Alpha_shape_3::Vertex_handle>;
+
+// gudhi type definition
+using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
+using Filtration_value = ST::Filtration_value;
+using Simplex_tree_vertex = ST::Vertex_handle;
+using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>;
+using Simplex_tree_vector_vertex = std::vector<Simplex_tree_vertex>;
+using Persistent_cohomology =
+    Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>;
+
+void usage(const std::string& progName) {
+  std::cerr << "Usage: " << progName << " path_to_the_OFF_file path_to_weight_file path_to_the_cuboid_file "
+                                        "coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n";
+  exit(-1);
+}
+
+int main(int argc, char* const argv[]) {
+  // program args management
+  if (argc != 6) {
+    std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n";
+    usage(argv[0]);
+  }
+
+  int coeff_field_characteristic = atoi(argv[4]);
+  Filtration_value min_persistence = strtof(argv[5], nullptr);
+
+  // Read points from file
+  std::string offInputFile(argv[1]);
+  // Read the OFF file (input file name given as parameter) and triangulate points
+  Gudhi::Points_3D_off_reader<Point_3> off_reader(offInputFile);
+  // Check the read operation was correct
+  if (!off_reader.is_valid()) {
+    std::cerr << "Unable to read file " << offInputFile << std::endl;
+    usage(argv[0]);
+  }
+
+  // Retrieve the points
+  std::vector<Point_3> lp = off_reader.get_point_cloud();
+
+  // Read iso_cuboid_3 information from file
+  std::ifstream iso_cuboid_str(argv[3]);
+  double x_min, y_min, z_min, x_max, y_max, z_max;
+  if (iso_cuboid_str.is_open()) {
+    if (!(iso_cuboid_str >> x_min >> y_min >> z_min >> x_max >> y_max >> z_max)) {
+      std::cerr << argv[3] << " - Bad file format." << std::endl;
+      usage(argv[0]);
+    }
+
+  } else {
+    std::cerr << "Unable to read file " << argv[3] << std::endl;
+    usage(argv[0]);
+  }
+  // Checking if the cuboid is the same in x,y and z direction. If not, CGAL will not process it.
+  if ((x_max - x_min != y_max - y_min) || (x_max - x_min != z_max - z_min) || (z_max - z_min != y_max - y_min)) {
+    std::cerr << "The size of the cuboid in every directions is not the same." << std::endl;
+    exit(-1);
+  }
+
+  double maximal_possible_weight = 0.015625 * (x_max - x_min) * (x_max - x_min);
+
+  // Read weights information from file
+  std::ifstream weights_ifstr(argv[2]);
+  std::vector<Weighted_point_3> wp;
+  if (weights_ifstr.is_open()) {
+    double weight = 0.0;
+    std::size_t index = 0;
+    wp.reserve(lp.size());
+    // Attempt read the weight in a double format, return false if it fails
+    while ((weights_ifstr >> weight) && (index < lp.size())) {
+      if ((weight >= maximal_possible_weight) || (weight < 0)) {
+        std::cerr << "At line " << (index + 1) << ", the weight (" << weight
+                  << ") is negative or more than or equal to maximal possible weight (" << maximal_possible_weight
+                  << ") = 1/64*cuboid length squared, which is not an acceptable input." << std::endl;
+        exit(-1);
+      }
+
+      wp.push_back(Weighted_point_3(lp[index], weight));
+      index++;
+    }
+    if (index != lp.size()) {
+      std::cerr << "Bad number of weights in file " << argv[2] << std::endl;
+      usage(argv[0]);
+    }
+  } else {
+    std::cerr << "Unable to read file " << argv[2] << std::endl;
+    usage(argv[0]);
+  }
+
+  // Define the periodic cube
+  P3RT3 prt(PK::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max));
+  // Heuristic for inserting large point sets (if pts is reasonably large)
+  prt.insert(wp.begin(), wp.end(), true);
+  // As prt won't be modified anymore switch to 1-sheeted cover if possible
+  if (prt.is_triangulation_in_1_sheet()) {
+    prt.convert_to_1_sheeted_covering();
+  } else {
+    std::cerr << "ERROR: we were not able to construct a triangulation within a single periodic domain." << std::endl;
+    exit(-1);
+  }
+  std::cout << "Weighted Periodic Delaunay computed." << std::endl;
+
+  // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. This is the default mode
+  // Maybe need to set it to GENERAL mode
+  Alpha_shape_3 as(prt, 0, Alpha_shape_3::GENERAL);
+
+  // filtration with alpha values from alpha shape
+  std::vector<Object> the_objects;
+  std::vector<Alpha_value_type> the_alpha_values;
+
+  Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
+                                                                  std::back_inserter(the_alpha_values));
+
+  as.filtration_with_alpha_values(disp);
+#ifdef DEBUG_TRACES
+  std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
+#endif  // DEBUG_TRACES
+
+  Alpha_shape_3::size_type count_vertices = 0;
+  Alpha_shape_3::size_type count_edges = 0;
+  Alpha_shape_3::size_type count_facets = 0;
+  Alpha_shape_3::size_type count_cells = 0;
+
+  // Loop on objects vector
+  Vertex_list vertex_list;
+  ST simplex_tree;
+  Alpha_shape_simplex_tree_map map_cgal_simplex_tree;
+  std::vector<Alpha_value_type>::iterator the_alpha_value_iterator = the_alpha_values.begin();
+  for (auto object_iterator : the_objects) {
+    // Retrieve Alpha shape vertex list from object
+    if (const Cell_handle* cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
+      vertex_list = from_cell<Vertex_list, Cell_handle>(*cell);
+      count_cells++;
+    } else if (const Facet* facet = CGAL::object_cast<Facet>(&object_iterator)) {
+      vertex_list = from_facet<Vertex_list, Facet>(*facet);
+      count_facets++;
+    } else if (const Edge_3* edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
+      vertex_list = from_edge<Vertex_list, Edge_3>(*edge);
+      count_edges++;
+    } else if (const Vertex_handle* vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) {
+      count_vertices++;
+      vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex);
+    }
+    // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex
+    Simplex_tree_vector_vertex the_simplex;
+    for (auto the_alpha_shape_vertex : vertex_list) {
+      Alpha_shape_simplex_tree_map::iterator the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex);
+      if (the_map_iterator == map_cgal_simplex_tree.end()) {
+        // alpha shape not found
+        Simplex_tree_vertex vertex = map_cgal_simplex_tree.size();
+#ifdef DEBUG_TRACES
+        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl;
+#endif  // DEBUG_TRACES
+        the_simplex.push_back(vertex);
+        map_cgal_simplex_tree.emplace(the_alpha_shape_vertex, vertex);
+      } else {
+        // alpha shape found
+        Simplex_tree_vertex vertex = the_map_iterator->second;
+#ifdef DEBUG_TRACES
+        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl;
+#endif  // DEBUG_TRACES
+        the_simplex.push_back(vertex);
+      }
+    }
+    // Construction of the simplex_tree
+    Filtration_value filtr = /*std::sqrt*/ (*the_alpha_value_iterator);
+#ifdef DEBUG_TRACES
+    std::cout << "filtration = " << filtr << std::endl;
+#endif  // DEBUG_TRACES
+    simplex_tree.insert_simplex(the_simplex, filtr);
+    if (the_alpha_value_iterator != the_alpha_values.end())
+      ++the_alpha_value_iterator;
+    else
+      std::cout << "This shall not happen" << std::endl;
+  }
+
+#ifdef DEBUG_TRACES
+  std::cout << "vertices \t\t" << count_vertices << std::endl;
+  std::cout << "edges \t\t" << count_edges << std::endl;
+  std::cout << "facets \t\t" << count_facets << std::endl;
+  std::cout << "cells \t\t" << count_cells << std::endl;
+
+  std::cout << "Information of the Simplex Tree: " << std::endl;
+  std::cout << "  Number of vertices = " << simplex_tree.num_vertices() << " ";
+  std::cout << "  Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl;
+  std::cout << "  Dimension = " << simplex_tree.dimension() << " ";
+#endif  // DEBUG_TRACES
+
+#ifdef DEBUG_TRACES
+  std::cout << "Iterator on vertices: " << std::endl;
+  for (auto vertex : simplex_tree.complex_vertex_range()) {
+    std::cout << vertex << " ";
+  }
+#endif  // DEBUG_TRACES
+
+  // Sort the simplices in the order of the filtration
+  simplex_tree.initialize_filtration();
+
+  std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl;
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(simplex_tree, true);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(coeff_field_characteristic);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  pcoh.output_diagram();
+
+  return 0;
+}
diff --git a/utilities/Bitmap_cubical_complex/CMakeLists.txt b/utilities/Bitmap_cubical_complex/CMakeLists.txt
new file mode 100644
index 00000000..676a730a
--- /dev/null
+++ b/utilities/Bitmap_cubical_complex/CMakeLists.txt
@@ -0,0 +1,29 @@
+cmake_minimum_required(VERSION 2.6)
+project(Bitmap_cubical_complex_utilities)
+
+add_executable ( cubical_complex_persistence cubical_complex_persistence.cpp )
+if (TBB_FOUND)
+  target_link_libraries(cubical_complex_persistence ${TBB_LIBRARIES})
+endif()
+
+add_test(NAME Bitmap_cubical_complex_utility_persistence_one_sphere COMMAND $<TARGET_FILE:cubical_complex_persistence>
+    "${CMAKE_SOURCE_DIR}/data/bitmap/CubicalOneSphere.txt")
+
+add_test(NAME Bitmap_cubical_complex_utility_persistence_two_sphere COMMAND $<TARGET_FILE:cubical_complex_persistence>
+    "${CMAKE_SOURCE_DIR}/data/bitmap/CubicalTwoSphere.txt")
+
+add_executable ( periodic_cubical_complex_persistence periodic_cubical_complex_persistence.cpp )
+if (TBB_FOUND)
+  target_link_libraries(periodic_cubical_complex_persistence ${TBB_LIBRARIES})
+endif()
+
+add_test(NAME Bitmap_cubical_complex_utility_periodic_boundary_conditions_2d_torus
+    COMMAND $<TARGET_FILE:periodic_cubical_complex_persistence>
+    "${CMAKE_SOURCE_DIR}/data/bitmap/2d_torus.txt")
+
+add_test(NAME Bitmap_cubical_complex_utility_periodic_boundary_conditions_3d_torus
+    COMMAND $<TARGET_FILE:periodic_cubical_complex_persistence>
+    "${CMAKE_SOURCE_DIR}/data/bitmap/3d_torus.txt")
+
+install(TARGETS cubical_complex_persistence DESTINATION bin)
+install(TARGETS periodic_cubical_complex_persistence DESTINATION bin)
diff --git a/utilities/Bitmap_cubical_complex/cubical_complex_persistence.cpp b/utilities/Bitmap_cubical_complex/cubical_complex_persistence.cpp
new file mode 100644
index 00000000..9d1bc08c
--- /dev/null
+++ b/utilities/Bitmap_cubical_complex/cubical_complex_persistence.cpp
@@ -0,0 +1,80 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2015  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 <gudhi/reader_utils.h>
+#include <gudhi/Bitmap_cubical_complex.h>
+#include <gudhi/Persistent_cohomology.h>
+
+// standard stuff
+#include <iostream>
+#include <string>
+#include <vector>
+#include <cstddef>
+
+int main(int argc, char** argv) {
+  std::cout
+      << "This program computes persistent homology, by using bitmap_cubical_complex class, of cubical "
+      << "complexes provided in text files in Perseus style (the only numbered in the first line is a dimension D of a"
+      << "bitmap. In the lines I between 2 and D+1 there are numbers of top dimensional cells in the direction I. Let "
+      << "N denote product of the numbers in the lines between 2 and D. In the lines D+2 to D+2+N there are "
+      << "filtrations of top dimensional cells. We assume that the cells are in the lexicographical order. See "
+      << "CubicalOneSphere.txt or CubicalTwoSphere.txt for example.\n"
+      << std::endl;
+
+  if (argc != 2) {
+    std::cerr << "Wrong number of parameters. Please provide the name of a file with a Perseus style bitmap at "
+              << "the input. The program will now terminate.\n";
+    return 1;
+  }
+
+  typedef Gudhi::cubical_complex::Bitmap_cubical_complex_base<double> Bitmap_cubical_complex_base;
+  typedef Gudhi::cubical_complex::Bitmap_cubical_complex<Bitmap_cubical_complex_base> Bitmap_cubical_complex;
+  typedef Gudhi::persistent_cohomology::Field_Zp Field_Zp;
+  typedef Gudhi::persistent_cohomology::Persistent_cohomology<Bitmap_cubical_complex, Field_Zp> Persistent_cohomology;
+
+  Bitmap_cubical_complex b(argv[1]);
+
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(b);
+  int p = 11;
+  double min_persistence = 0;
+
+  pcoh.init_coefficients(p);  // initializes the coefficient field for homology
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  std::string output_file_name(argv[1]);
+  output_file_name += "_persistence";
+
+  std::size_t last_in_path = output_file_name.find_last_of("/\\");
+
+  if (last_in_path != std::string::npos) {
+    output_file_name = output_file_name.substr(last_in_path + 1);
+  }
+
+  std::ofstream out(output_file_name.c_str());
+  pcoh.output_diagram(out);
+  out.close();
+
+  std::cout << "Result in file: " << output_file_name << "\n";
+
+  return 0;
+}
diff --git a/utilities/Bitmap_cubical_complex/cubicalcomplex.md b/utilities/Bitmap_cubical_complex/cubicalcomplex.md
new file mode 100644
index 00000000..6e1b2578
--- /dev/null
+++ b/utilities/Bitmap_cubical_complex/cubicalcomplex.md
@@ -0,0 +1,29 @@
+

+

+# Cubical complex#

+

+## cubical_complex_persistence ##

+This program computes persistent homology, by using the Bitmap_cubical_complex class, of cubical complexes provided in text files in Perseus style.

+See [here](/doc/latest/fileformats.html#FileFormatsPerseus) for a description of the file format.

+

+**Example**

+

+```

+   cubical_complex_persistence data/bitmap/CubicalTwoSphere.txt

+```

+

+* Creates a Cubical Complex from the Perseus style file `CubicalTwoSphere.txt`,

+computes Persistence cohomology from it and writes the results in a persistence file `CubicalTwoSphere.txt_persistence`.

+

+## periodic_cubical_complex_persistence ##

+

+Same as above, but with periodic boundary conditions.

+

+**Example**

+

+```

+   periodic_cubical_complex_persistence data/bitmap/3d_torus.txt

+```

+

+* Creates a Periodical Cubical Complex from the Perseus style file `3d_torus.txt`,

+computes Persistence cohomology from it and writes the results in a persistence file `3d_torus.txt_persistence`.

diff --git a/utilities/Bitmap_cubical_complex/periodic_cubical_complex_persistence.cpp b/utilities/Bitmap_cubical_complex/periodic_cubical_complex_persistence.cpp
new file mode 100644
index 00000000..c812cb3a
--- /dev/null
+++ b/utilities/Bitmap_cubical_complex/periodic_cubical_complex_persistence.cpp
@@ -0,0 +1,82 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2015  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 <gudhi/reader_utils.h>
+#include <gudhi/Bitmap_cubical_complex.h>
+#include <gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h>
+#include <gudhi/Persistent_cohomology.h>
+
+// standard stuff
+#include <iostream>
+#include <sstream>
+#include <vector>
+#include <string>
+
+int main(int argc, char** argv) {
+  std::cout
+      << "This program computes persistent homology, by using "
+      << "Bitmap_cubical_complex_periodic_boundary_conditions class, of cubical complexes provided in text files in "
+      << "Perseus style (the only numbered in the first line is a dimension D of a bitmap. In the lines I between 2 "
+      << "and D+1 there are numbers of top dimensional cells in the direction I. Let N denote product of the numbers "
+      << "in the lines between 2 and D. In the lines D+2 to D+2+N there are filtrations of top dimensional cells. We "
+      << "assume that the cells are in the lexicographical order. See CubicalOneSphere.txt or CubicalTwoSphere.txt for"
+      << " example.\n"
+      << std::endl;
+
+  if (argc != 2) {
+    std::cerr << "Wrong number of parameters. Please provide the name of a file with a Perseus style bitmap at "
+              << "the input. The program will now terminate.\n";
+    return 1;
+  }
+
+  typedef Gudhi::cubical_complex::Bitmap_cubical_complex_periodic_boundary_conditions_base<double> Bitmap_base;
+  typedef Gudhi::cubical_complex::Bitmap_cubical_complex<Bitmap_base> Bitmap_cubical_complex;
+
+  Bitmap_cubical_complex b(argv[1]);
+
+  typedef Gudhi::persistent_cohomology::Field_Zp Field_Zp;
+  typedef Gudhi::persistent_cohomology::Persistent_cohomology<Bitmap_cubical_complex, Field_Zp> Persistent_cohomology;
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(b, true);
+
+  int p = 11;
+  double min_persistence = 0;
+  pcoh.init_coefficients(p);  // initializes the coefficient field for homology
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  std::string output_file_name(argv[1]);
+  output_file_name += "_persistence";
+
+  std::size_t last_in_path = output_file_name.find_last_of("/\\");
+
+  if (last_in_path != std::string::npos) {
+    output_file_name = output_file_name.substr(last_in_path + 1);
+  }
+
+  std::ofstream out(output_file_name.c_str());
+  pcoh.output_diagram(out);
+  out.close();
+
+  std::cout << "Result in file: " << output_file_name << "\n";
+
+  return 0;
+}
diff --git a/utilities/Bottleneck_distance/CMakeLists.txt b/utilities/Bottleneck_distance/CMakeLists.txt
new file mode 100644
index 00000000..d19e3b1c
--- /dev/null
+++ b/utilities/Bottleneck_distance/CMakeLists.txt
@@ -0,0 +1,16 @@
+cmake_minimum_required(VERSION 2.6)
+project(Bottleneck_distance_utilities)
+
+if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1)
+  add_executable (bottleneck_distance bottleneck_distance.cpp)
+  if (TBB_FOUND)
+    target_link_libraries(bottleneck_distance ${TBB_LIBRARIES})
+  endif(TBB_FOUND)
+
+  add_test(NAME Bottleneck_distance_utilities_Bottleneck_read_file
+      COMMAND $<TARGET_FILE:bottleneck_distance>
+      "${CMAKE_SOURCE_DIR}/data/persistence_diagram/first.pers" "${CMAKE_SOURCE_DIR}/data/persistence_diagram/second.pers")
+
+  install(TARGETS bottleneck_distance DESTINATION bin)
+
+endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1)
diff --git a/utilities/Bottleneck_distance/bottleneck_distance.cpp b/utilities/Bottleneck_distance/bottleneck_distance.cpp
new file mode 100644
index 00000000..9dd52b31
--- /dev/null
+++ b/utilities/Bottleneck_distance/bottleneck_distance.cpp
@@ -0,0 +1,50 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Authors:       Francois Godi, small modifications by Pawel Dlotko
+ *
+ *    Copyright (C) 2015  INRIA
+ *
+ *    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 <gudhi/Bottleneck.h>
+#include <gudhi/reader_utils.h>
+#include <iostream>
+#include <vector>
+#include <utility>  // for pair
+#include <string>
+#include <limits>  // for numeric_limits
+
+int main(int argc, char** argv) {
+  if (argc < 3) {
+    std::cout << "To run this program please provide as an input two files with persistence diagrams. Each file" <<
+        " should contain a birth-death pair per line. Third, optional parameter is an error bound on the bottleneck" <<
+        " distance (set by default to the smallest positive double value). If you set the error bound to 0, be" <<
+        " aware this version is exact but expensive. The program will now terminate \n";
+    return -1;
+  }
+  std::vector<std::pair<double, double>> diag1 = Gudhi::read_persistence_intervals_in_dimension(argv[1]);
+  std::vector<std::pair<double, double>> diag2 = Gudhi::read_persistence_intervals_in_dimension(argv[2]);
+
+  double tolerance = std::numeric_limits<double>::min();
+  if (argc == 4) {
+    tolerance = atof(argv[3]);
+  }
+  double b = Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, tolerance);
+  std::cout << "The distance between the diagrams is : " << b << ". The tolerance is : " << tolerance << std::endl;
+
+  return 0;
+}
diff --git a/utilities/Bottleneck_distance/bottleneckdistance.md b/utilities/Bottleneck_distance/bottleneckdistance.md
new file mode 100644
index 00000000..526f5822
--- /dev/null
+++ b/utilities/Bottleneck_distance/bottleneckdistance.md
@@ -0,0 +1,18 @@
+

+

+# Bottleneck distance #

+

+## bottleneck_read_file_example ##

+

+This program computes the Bottleneck distance between two persistence diagram files.

+

+**Usage**

+

+```

+   bottleneck_read_file_example <file_1.pers> <file_2.pers> [<tolerance>]

+```

+

+where

+

+* `<file_1.pers>` and `<file_2.pers>` must be in the format described [here](/doc/latest/fileformats.html#FileFormatsPers).

+* `<tolerance>` is an error bound on the bottleneck distance (set by default to the smallest positive double value).

diff --git a/utilities/Nerve_GIC/CMakeLists.txt b/utilities/Nerve_GIC/CMakeLists.txt
new file mode 100644
index 00000000..7762c8a0
--- /dev/null
+++ b/utilities/Nerve_GIC/CMakeLists.txt
@@ -0,0 +1,24 @@
+cmake_minimum_required(VERSION 2.6)
+project(Nerve_GIC_examples)
+
+if (NOT CGAL_VERSION VERSION_LESS 4.8.1)
+
+  add_executable ( Nerve Nerve.cpp )
+  add_executable ( VoronoiGIC VoronoiGIC.cpp )
+
+  if (TBB_FOUND)
+    target_link_libraries(Nerve ${TBB_LIBRARIES})
+    target_link_libraries(VoronoiGIC ${TBB_LIBRARIES})
+  endif()
+
+  file(COPY KeplerMapperVisuFromTxtFile.py km.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
+  # Copy files for not to pollute sources when testing
+  file(COPY "${CMAKE_SOURCE_DIR}/data/points/human.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
+
+  add_test(NAME Nerve_GIC_utilities_nerve COMMAND $<TARGET_FILE:Nerve>
+      "human.off" "2" "10" "0.3")
+
+  add_test(NAME Nerve_GIC_utilities_VoronoiGIC COMMAND $<TARGET_FILE:VoronoiGIC>
+      "human.off" "100")
+
+endif (NOT CGAL_VERSION VERSION_LESS 4.8.1)
diff --git a/utilities/Nerve_GIC/KeplerMapperVisuFromTxtFile.py b/utilities/Nerve_GIC/KeplerMapperVisuFromTxtFile.py
new file mode 100755
index 00000000..c811f610
--- /dev/null
+++ b/utilities/Nerve_GIC/KeplerMapperVisuFromTxtFile.py
@@ -0,0 +1,89 @@
+#!/usr/bin/env python
+
+import km
+import numpy as np
+from collections import defaultdict
+import argparse
+
+"""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):       Mathieu Carriere
+
+   Copyright (C) 2017 INRIA
+
+   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/>.
+"""
+
+__author__ = "Mathieu Carriere"
+__copyright__ = "Copyright (C) 2017 INRIA"
+__license__ = "GPL v3"
+
+parser = argparse.ArgumentParser(description='Creates an html Keppler Mapper '
+                                 'file to visualize a SC.txt file.',
+                                 epilog='Example: '
+                                 './KeplerMapperVisuFromTxtFile.py '
+                                 '-f ../../data/points/human.off_sc.txt'
+                                 '- Constructs an human.off_sc.html file.')
+parser.add_argument("-f", "--file", type=str, required=True)
+
+args = parser.parse_args()
+
+with open(args.file, 'r') as f:
+    network = {}
+    mapper = km.KeplerMapper(verbose=0)
+    data = np.zeros((3,3))
+    projected_data = mapper.fit_transform( data, projection="sum", scaler=None )
+
+    nodes = defaultdict(list)
+    links = defaultdict(list)
+    custom = defaultdict(list)
+
+    dat = f.readline()
+    lens = f.readline()
+    color = f.readline();
+    param = [float(i) for i in f.readline().split(" ")]
+
+    nums = [int(i) for i in f.readline().split(" ")]
+    num_nodes = nums[0]
+    num_edges = nums[1]
+
+    for i in range(0,num_nodes):
+        point = [float(j) for j in f.readline().split(" ")]
+        nodes[  str(int(point[0]))  ] = [  int(point[0]), point[1], int(point[2])  ]
+        links[  str(int(point[0]))  ] = []
+        custom[  int(point[0])  ] = point[1]
+
+    m = min([custom[i] for i in range(0,num_nodes)])
+    M = max([custom[i] for i in range(0,num_nodes)])
+
+    for i in range(0,num_edges):
+        edge = [int(j) for j in f.readline().split(" ")]
+        links[  str(edge[0])  ].append(  str(edge[1])  )
+        links[  str(edge[1])  ].append(  str(edge[0])  )
+
+    network["nodes"] = nodes
+    network["links"] = links
+    network["meta"] = lens
+
+    html_output_filename = args.file.rsplit('.', 1)[0] + '.html'
+    mapper.visualize(network, color_function = color, path_html=html_output_filename, title=dat,
+    graph_link_distance=30, graph_gravity=0.1, graph_charge=-120, custom_tooltips=custom, width_html=0,
+    height_html=0, show_tooltips=True, show_title=True, show_meta=True, res=param[0],gain=param[1], minimum=m,maximum=M)
+    message = repr(html_output_filename) + " is generated. You can now use your favorite web browser to visualize it."
+    print(message)
+
+
+    f.close()
diff --git a/utilities/Nerve_GIC/Nerve.cpp b/utilities/Nerve_GIC/Nerve.cpp
new file mode 100644
index 00000000..aefc3874
--- /dev/null
+++ b/utilities/Nerve_GIC/Nerve.cpp
@@ -0,0 +1,96 @@
+/*    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):       Mathieu Carrière
+ *
+ *    Copyright (C) 2017  INRIA
+ *
+ *    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 <gudhi/GIC.h>
+
+#include <string>
+#include <vector>
+
+void usage(int nbArgs, char *const progName) {
+  std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n";
+  std::cerr << "Usage: " << progName << " filename.off coordinate resolution gain [-v] \n";
+  std::cerr << "       i.e.: " << progName << " ../../data/points/human.off 2 10 0.3 -v \n";
+  exit(-1);  // ----- >>
+}
+
+int main(int argc, char **argv) {
+  if ((argc != 5) && (argc != 6)) usage(argc, argv[0]);
+
+  using Point = std::vector<float>;
+
+  std::string off_file_name(argv[1]);
+  int coord = atoi(argv[2]);
+  int resolution = atoi(argv[3]);
+  double gain = atof(argv[4]);
+  bool verb = 0;
+  if (argc == 6) verb = 1;
+
+  // --------------------------------
+  // Init of a Nerve from an OFF file
+  // --------------------------------
+
+  Gudhi::cover_complex::Cover_complex<Point> SC;
+  SC.set_verbose(verb);
+
+  bool check = SC.read_point_cloud(off_file_name);
+
+  if (!check) {
+    std::cout << "Incorrect OFF file." << std::endl;
+  } else {
+    SC.set_type("Nerve");
+
+    SC.set_color_from_coordinate(coord);
+    SC.set_function_from_coordinate(coord);
+
+    SC.set_graph_from_OFF();
+    SC.set_resolution_with_interval_number(resolution);
+    SC.set_gain(gain);
+    SC.set_cover_from_function();
+
+    SC.find_simplices();
+
+    SC.write_info();
+
+    Gudhi::Simplex_tree<> stree;
+    SC.create_complex(stree);
+    SC.compute_PD();
+
+    // ----------------------------------------------------------------------------
+    // Display information about the graph induced complex
+    // ----------------------------------------------------------------------------
+
+    if (verb) {
+      std::cout << "Nerve is of dimension " << stree.dimension() << " - " << stree.num_simplices() << " simplices - "
+                << stree.num_vertices() << " vertices." << std::endl;
+
+      std::cout << "Iterator on Nerve simplices" << std::endl;
+      for (auto f_simplex : stree.filtration_simplex_range()) {
+        for (auto vertex : stree.simplex_vertex_range(f_simplex)) {
+          std::cout << vertex << " ";
+        }
+        std::cout << std::endl;
+      }
+    }
+  }
+
+  return 0;
+}
diff --git a/utilities/Nerve_GIC/Nerve.txt b/utilities/Nerve_GIC/Nerve.txt
new file mode 100644
index 00000000..839ff45e
--- /dev/null
+++ b/utilities/Nerve_GIC/Nerve.txt
@@ -0,0 +1,63 @@
+Min function value = -0.979672 and Max function value = 0.816414
+Interval 0 = [-0.979672, -0.761576]
+Interval 1 = [-0.838551, -0.581967]
+Interval 2 = [-0.658942, -0.402359]
+Interval 3 = [-0.479334, -0.22275]
+Interval 4 = [-0.299725, -0.0431415]
+Interval 5 = [-0.120117, 0.136467]
+Interval 6 = [0.059492, 0.316076]
+Interval 7 = [0.239101, 0.495684]
+Interval 8 = [0.418709, 0.675293]
+Interval 9 = [0.598318, 0.816414]
+Computing preimages...
+Computing connected components...
+.txt generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox.
+5 interval(s) in dimension 0:
+  [-0.909111, 0.00817529]
+  [-0.171433, 0.367392]
+  [-0.171433, 0.367392]
+  [-0.909111, 0.745853]
+0 interval(s) in dimension 1:
+Nerve is of dimension 1 - 41 simplices - 21 vertices.
+Iterator on Nerve simplices
+1 
+0 
+4 
+4 0 
+2 
+2 1 
+8 
+8 2 
+5 
+5 4 
+9 
+9 8 
+13 
+13 5 
+14 
+14 9 
+19 
+19 13 
+25 
+32 
+20 
+32 20 
+33 
+33 25 
+26 
+26 14 
+26 19 
+42 
+42 26 
+34 
+34 33 
+27 
+27 20 
+35 
+35 27 
+35 34 
+42 35 
+44 
+44 35 
+54 
+54 44 
\ No newline at end of file
diff --git a/utilities/Nerve_GIC/VoronoiGIC.cpp b/utilities/Nerve_GIC/VoronoiGIC.cpp
new file mode 100644
index 00000000..54bb871e
--- /dev/null
+++ b/utilities/Nerve_GIC/VoronoiGIC.cpp
@@ -0,0 +1,90 @@
+/*    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):       Mathieu Carrière
+ *
+ *    Copyright (C) 2017  INRIA
+ *
+ *    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 <gudhi/GIC.h>
+
+#include <string>
+#include <vector>
+
+void usage(int nbArgs, char *const progName) {
+  std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n";
+  std::cerr << "Usage: " << progName << " filename.off N [-v] \n";
+  std::cerr << "       i.e.: " << progName << " ../../data/points/human.off 100 -v \n";
+  exit(-1);  // ----- >>
+}
+
+int main(int argc, char **argv) {
+  if ((argc != 3) && (argc != 4)) usage(argc, argv[0]);
+
+  using Point = std::vector<float>;
+
+  std::string off_file_name(argv[1]);
+  int m = atoi(argv[2]);
+  bool verb = 0;
+  if (argc == 4) verb = 1;
+
+  // ----------------------------------------------------------------------------
+  // Init of a graph induced complex from an OFF file
+  // ----------------------------------------------------------------------------
+
+  Gudhi::cover_complex::Cover_complex<Point> GIC;
+  GIC.set_verbose(verb);
+
+  bool check = GIC.read_point_cloud(off_file_name);
+
+  if (!check) {
+    std::cout << "Incorrect OFF file." << std::endl;
+  } else {
+    GIC.set_type("GIC");
+
+    GIC.set_color_from_coordinate();
+
+    GIC.set_graph_from_OFF();
+    GIC.set_cover_from_Voronoi(Gudhi::Euclidean_distance(), m);
+
+    GIC.find_simplices();
+
+    GIC.plot_OFF();
+
+    Gudhi::Simplex_tree<> stree;
+    GIC.create_complex(stree);
+
+    // ----------------------------------------------------------------------------
+    // Display information about the graph induced complex
+    // ----------------------------------------------------------------------------
+
+    if (verb) {
+      std::cout << "Graph induced complex is of dimension " << stree.dimension() << " - " << stree.num_simplices()
+                << " simplices - " << stree.num_vertices() << " vertices." << std::endl;
+
+      std::cout << "Iterator on graph induced complex simplices" << std::endl;
+      for (auto f_simplex : stree.filtration_simplex_range()) {
+        for (auto vertex : stree.simplex_vertex_range(f_simplex)) {
+          std::cout << vertex << " ";
+        }
+        std::cout << std::endl;
+      }
+    }
+  }
+
+  return 0;
+}
diff --git a/utilities/Nerve_GIC/covercomplex.md b/utilities/Nerve_GIC/covercomplex.md
new file mode 100644
index 00000000..f33cb2e0
--- /dev/null
+++ b/utilities/Nerve_GIC/covercomplex.md
@@ -0,0 +1,62 @@
+
+
+# Cover complex #
+
+
+## Nerve ##
+This program builds the Nerve of a point cloud sampled on an OFF file.
+The cover C comes from the preimages of intervals covering a coordinate function,
+which are then refined into their connected components using the triangulation of the .OFF file.
+
+The program also writes a file SC.txt.
+The first three lines in this file are the location of the input point cloud and the function used to compute the cover.
+The fourth line contains the number of vertices nv and edges ne of the Nerve. The next nv lines represent the vertices.
+Each line contains the vertex ID, the number of data points it contains, and their average color function value.
+Finally, the next ne lines represent the edges, characterized by the ID of their vertices.
+
+**Usage**
+
+`Nerve <OFF input file> coordinate resolution gain [-v]`
+
+where
+
+* `coordinate` is the coordinate function to cover
+* `resolution` is the number of the intervals
+* `gain` is the gain for each interval
+* `-v` is optional, it activates verbose mode.
+
+**Example**
+
+`Nerve ../../data/points/human.off 2 10 0.3`
+
+* Builds the Nerve of a point cloud sampled on a 3D human shape (human.off).
+The cover C comes from the preimages of intervals (10 intervals with gain 0.3) covering the height function (coordinate 2).
+
+`python KeplerMapperVisuFromTxtFile.py -f ../../data/points/human.off_sc.txt`
+
+* Constructs `human.off_sc.html` file. You can now use your favorite web browser to visualize it.
+
+## VoronoiGIC ##
+
+This util builds the Graph Induced Complex (GIC) of a point cloud.
+It subsamples *N* points in the point cloud, which act as seeds of a geodesic Voronoï diagram.
+Each cell of the diagram is then an element of C.
+
+The program also writes a file `*_sc.off`, that is an OFF file that can be visualized with GeomView.
+
+**Usage**
+
+`VoroniGIC <OFF input file> samples_number [-v]`
+
+where
+
+* `samples_number` is the number of samples to take from the point cloud
+* `-v` is optional, it activates verbose mode.
+
+**Example**
+
+`VoroniGIC ../../data/points/human.off 700`
+
+* Builds the Voronoi Graph Induced Complex with 700 subsamples from `human.off` file.
+`../../data/points/human_sc.off` can be visualized with GeomView.
+
diff --git a/utilities/Nerve_GIC/km.py b/utilities/Nerve_GIC/km.py
new file mode 100755
index 00000000..53024aab
--- /dev/null
+++ b/utilities/Nerve_GIC/km.py
@@ -0,0 +1,390 @@
+from __future__ import division
+import numpy as np
+from collections import defaultdict
+import json
+import itertools
+from sklearn import cluster, preprocessing, manifold
+from datetime import datetime
+import sys
+
+class KeplerMapper(object):
+  # With this class you can build topological networks from (high-dimensional) data.
+  #
+  # 1)   	Fit a projection/lens/function to a dataset and transform it. 
+  #     	For instance "mean_of_row(x) for x in X"
+  # 2)   	Map this projection with overlapping intervals/hypercubes. 
+  #    		Cluster the points inside the interval 
+  #    		(Note: we cluster on the inverse image/original data to lessen projection loss).
+  #    		If two clusters/nodes have the same members (due to the overlap), then: 
+  #    		connect these with an edge.
+  # 3)  	Visualize the network using HTML and D3.js.
+  # 
+  # functions
+  # ---------
+  # fit_transform:   Create a projection (lens) from a dataset
+  # map:         	Apply Mapper algorithm on this projection and build a simplicial complex
+  # visualize:    	Turns the complex dictionary into a HTML/D3.js visualization
+  
+  def __init__(self, verbose=2):
+    self.verbose = verbose
+    
+    self.chunk_dist = []
+    self.overlap_dist = []
+    self.d = []
+    self.nr_cubes = 0
+    self.overlap_perc = 0
+    self.clusterer = False
+
+  def fit_transform(self, X, projection="sum", scaler=preprocessing.MinMaxScaler()):
+    # Creates the projection/lens from X. 
+    #
+    # Input:      X. Input features as a numpy array.
+    # Output:     projected_X. original data transformed to a projection (lens).
+    # 
+    # parameters
+    # ----------
+    # projection:   Projection parameter is either a string, 
+    #               a scikit class with fit_transform, like manifold.TSNE(), 
+    #               or a list of dimension indices.
+    # scaler:       if None, do no scaling, else apply scaling to the projection
+    #               Default: Min-Max scaling
+    
+    self.scaler = scaler
+    self.projection = str(projection)
+    
+    # Detect if projection is a class (for scikit-learn)
+    #if str(type(projection))[1:6] == "class": #TODO: de-ugly-fy
+    #  reducer = projection
+    #  if self.verbose > 0:
+    #    try:    
+    #      projection.set_params(**{"verbose":self.verbose})
+    #    except:
+    #      pass
+    #    print("\n..Projecting data using: \n\t%s\n"%str(projection))
+    #  X = reducer.fit_transform(X)
+    
+    # Detect if projection is a string (for standard functions)
+    if isinstance(projection, str):
+      if self.verbose > 0:
+        print("\n..Projecting data using: %s"%(projection))
+      # Stats lenses
+      if projection == "sum": # sum of row
+        X = np.sum(X, axis=1).reshape((X.shape[0],1))
+      if projection == "mean": # mean of row
+        X = np.mean(X, axis=1).reshape((X.shape[0],1))
+      if projection == "median": # mean of row
+        X = np.median(X, axis=1).reshape((X.shape[0],1))
+      if projection == "max": # max of row
+        X = np.max(X, axis=1).reshape((X.shape[0],1))
+      if projection == "min": # min of row
+        X = np.min(X, axis=1).reshape((X.shape[0],1))
+      if projection == "std": # std of row
+        X = np.std(X, axis=1).reshape((X.shape[0],1))
+        
+      if projection == "dist_mean": # Distance of x to mean of X
+        X_mean = np.mean(X, axis=0) 
+        X = np.sum(np.sqrt((X - X_mean)**2), axis=1).reshape((X.shape[0],1))
+
+    # Detect if projection is a list (with dimension indices)
+    if isinstance(projection, list):
+      if self.verbose > 0:
+        print("\n..Projecting data using: %s"%(str(projection)))
+      X = X[:,np.array(projection)]
+      
+    # Scaling
+    if scaler is not None:
+      if self.verbose > 0:
+        print("\n..Scaling with: %s\n"%str(scaler))
+      X = scaler.fit_transform(X)
+    
+    return X
+
+  def map(self, projected_X, inverse_X=None, clusterer=cluster.DBSCAN(eps=0.5,min_samples=3), nr_cubes=10, overlap_perc=0.1):
+    # This maps the data to a simplicial complex. Returns a dictionary with nodes and links.
+    # 
+    # Input:    projected_X. A Numpy array with the projection/lens. 
+    # Output:    complex. A dictionary with "nodes", "links" and "meta information"
+    #
+    # parameters
+    # ----------
+    # projected_X  	projected_X. A Numpy array with the projection/lens. Required.
+    # inverse_X    	Numpy array or None. If None then the projection itself is used for clustering.
+    # clusterer    	Scikit-learn API compatible clustering algorithm. Default: DBSCAN
+    # nr_cubes    	Int. The number of intervals/hypercubes to create.
+    # overlap_perc  Float. The percentage of overlap "between" the intervals/hypercubes.
+    
+    start = datetime.now()
+    
+    # Helper function
+    def cube_coordinates_all(nr_cubes, nr_dimensions):
+      # Helper function to get origin coordinates for our intervals/hypercubes
+      # Useful for looping no matter the number of cubes or dimensions
+      # Example:   	if there are 4 cubes per dimension and 3 dimensions 
+      #       		return the bottom left (origin) coordinates of 64 hypercubes, 
+      #       		as a sorted list of Numpy arrays
+      # TODO: elegance-ify...
+      l = []
+      for x in range(nr_cubes):
+        l += [x] * nr_dimensions
+      return [np.array(list(f)) for f in sorted(set(itertools.permutations(l,nr_dimensions)))]
+    
+    nodes = defaultdict(list)
+    links = defaultdict(list)
+    complex = {}
+    self.nr_cubes = nr_cubes
+    self.clusterer = clusterer
+    self.overlap_perc = overlap_perc
+    
+    if self.verbose > 0:
+      print("Mapping on data shaped %s using dimensions\n"%(str(projected_X.shape)))
+    
+    # If inverse image is not provided, we use the projection as the inverse image (suffer projection loss)
+    if inverse_X is None:
+      inverse_X = projected_X
+      
+    # We chop up the min-max column ranges into 'nr_cubes' parts
+    self.chunk_dist = (np.max(projected_X, axis=0) - np.min(projected_X, axis=0))/nr_cubes
+
+    # We calculate the overlapping windows distance 
+    self.overlap_dist = self.overlap_perc * self.chunk_dist
+
+    # We find our starting point
+    self.d = np.min(projected_X, axis=0)
+    
+    # Use a dimension index array on the projected X
+    # (For now this uses the entire dimensionality, but we keep for experimentation)
+    di = np.array([x for x in range(projected_X.shape[1])])
+    
+    # Prefix'ing the data with ID's
+    ids = np.array([x for x in range(projected_X.shape[0])])
+    projected_X = np.c_[ids,projected_X]
+    inverse_X = np.c_[ids,inverse_X]
+
+    # Subdivide the projected data X in intervals/hypercubes with overlap
+    if self.verbose > 0:
+      total_cubes = len(cube_coordinates_all(nr_cubes,projected_X.shape[1]))
+      print("Creating %s hypercubes."%total_cubes)
+
+    for i, coor in enumerate(cube_coordinates_all(nr_cubes,di.shape[0])):
+      # Slice the hypercube
+      hypercube = projected_X[ np.invert(np.any((projected_X[:,di+1] >= self.d[di] + (coor * self.chunk_dist[di])) & 
+          (projected_X[:,di+1] < self.d[di] + (coor * self.chunk_dist[di]) + self.chunk_dist[di] + self.overlap_dist[di]) == False, axis=1 )) ]
+      
+      if self.verbose > 1:
+        print("There are %s points in cube_%s / %s with starting range %s"%
+              (hypercube.shape[0],i,total_cubes,self.d[di] + (coor * self.chunk_dist[di])))
+      
+      # If at least one sample inside the hypercube
+      if hypercube.shape[0] > 0:
+        # Cluster the data point(s) in the cube, skipping the id-column
+        # Note that we apply clustering on the inverse image (original data samples) that fall inside the cube.
+        inverse_x = inverse_X[[int(nn) for nn in hypercube[:,0]]]
+        
+        clusterer.fit(inverse_x[:,1:])
+        
+        if self.verbose > 1:
+          print("Found %s clusters in cube_%s\n"%(np.unique(clusterer.labels_[clusterer.labels_ > -1]).shape[0],i))
+        
+        #Now for every (sample id in cube, predicted cluster label)
+        for a in np.c_[hypercube[:,0],clusterer.labels_]:
+          if a[1] != -1: #if not predicted as noise
+            cluster_id = str(coor[0])+"_"+str(i)+"_"+str(a[1])+"_"+str(coor)+"_"+str(self.d[di] + (coor * self.chunk_dist[di])) # TODO: de-rudimentary-ify
+            nodes[cluster_id].append( int(a[0]) ) # Append the member id's as integers
+      else:
+        if self.verbose > 1:
+          print("Cube_%s is empty.\n"%(i))
+
+    # Create links when clusters from different hypercubes have members with the same sample id.
+    candidates = itertools.combinations(nodes.keys(),2)
+    for candidate in candidates:
+      # if there are non-unique members in the union
+      if len(nodes[candidate[0]]+nodes[candidate[1]]) != len(set(nodes[candidate[0]]+nodes[candidate[1]])):
+        links[candidate[0]].append( candidate[1] )
+
+    # Reporting
+    if self.verbose > 0:
+      nr_links = 0
+      for k in links:
+        nr_links += len(links[k])
+      print("\ncreated %s edges and %s nodes in %s."%(nr_links,len(nodes),str(datetime.now()-start)))
+    
+    complex["nodes"] = nodes
+    complex["links"] = links
+    complex["meta"] = self.projection
+
+    return complex
+
+  def visualize(self, complex, color_function="", path_html="mapper_visualization_output.html", title="My Data", 
+          graph_link_distance=30, graph_gravity=0.1, graph_charge=-120, custom_tooltips=None, width_html=0, 
+          height_html=0, show_tooltips=True, show_title=True, show_meta=True, res=0,gain=0,minimum=0,maximum=0):
+    # Turns the dictionary 'complex' in a html file with d3.js
+    #
+    # Input:      complex. Dictionary (output from calling .map())
+    # Output:      a HTML page saved as a file in 'path_html'.
+    # 
+    # parameters
+    # ----------
+    # color_function    	string. Not fully implemented. Default: "" (distance to origin)
+    # path_html        		file path as string. Where to save the HTML page.
+    # title          		string. HTML page document title and first heading.
+    # graph_link_distance  	int. Edge length.
+    # graph_gravity     	float. "Gravity" to center of layout.
+    # graph_charge      	int. charge between nodes.
+    # custom_tooltips   	None or Numpy Array. You could use "y"-label array for this.
+    # width_html        	int. Width of canvas. Default: 0 (full width)
+    # height_html       	int. Height of canvas. Default: 0 (full height)
+    # show_tooltips     	bool. default:True
+    # show_title      		bool. default:True
+    # show_meta        		bool. default:True
+
+    # Format JSON for D3 graph
+    json_s = {}
+    json_s["nodes"] = []
+    json_s["links"] = []
+    k2e = {} # a key to incremental int dict, used for id's when linking
+
+    for e, k in enumerate(complex["nodes"]):
+      # Tooltip and node color formatting, TODO: de-mess-ify
+      if custom_tooltips is not None:
+        tooltip_s = "<h2>Cluster %s</h2>"%k + " ".join(str(custom_tooltips[complex["nodes"][k][0]]).split(" "))
+        if maximum == minimum:
+          tooltip_i = 0
+        else:
+          tooltip_i = int(30*(custom_tooltips[complex["nodes"][k][0]]-minimum)/(maximum-minimum))
+        json_s["nodes"].append({"name": str(k), "tooltip": tooltip_s, "group": 2 * int(np.log(complex["nodes"][k][2])), "color": tooltip_i})
+      else:
+        tooltip_s = "<h2>Cluster %s</h2>Contains %s members."%(k,len(complex["nodes"][k]))
+        json_s["nodes"].append({"name": str(k), "tooltip": tooltip_s, "group": 2 * int(np.log(len(complex["nodes"][k]))), "color": str(k.split("_")[0])})
+      k2e[k] = e
+    for k in complex["links"]:
+      for link in complex["links"][k]:
+        json_s["links"].append({"source": k2e[k], "target":k2e[link],"value":1})
+
+    # Width and height of graph in HTML output
+    if width_html == 0:
+      width_css = "100%"
+      width_js = 'document.getElementById("holder").offsetWidth-20'
+    else:
+      width_css = "%spx" % width_html
+      width_js = "%s" % width_html
+    if height_html == 0:
+      height_css = "100%"
+      height_js = 'document.getElementById("holder").offsetHeight-20'
+    else:
+      height_css = "%spx" % height_html
+      height_js = "%s" % height_html
+    
+    # Whether to show certain UI elements or not
+    if show_tooltips == False:
+      tooltips_display = "display: none;"
+    else:
+      tooltips_display = ""
+      
+    if show_meta == False:
+      meta_display = "display: none;"
+    else:
+      meta_display = ""
+      
+    if show_title == False:
+      title_display = "display: none;"
+    else:
+      title_display = ""  
+    
+    with open(path_html,"wb") as outfile:
+      html = """<!DOCTYPE html>
+    <meta charset="utf-8">
+    <meta name="generator" content="KeplerMapper">
+    <title>%s | KeplerMapper</title>
+    <link href='https://fonts.googleapis.com/css?family=Roboto:700,300' rel='stylesheet' type='text/css'>
+    <style>
+    * {margin: 0; padding: 0;}
+    html { height: 100%%;}
+    body {background: #111; height: 100%%; font: 100 16px Roboto, Sans-serif;}
+    .link { stroke: #999; stroke-opacity: .333;  }
+    .divs div { border-radius: 50%%; background: red; position: absolute; }
+    .divs { position: absolute; top: 0; left: 0; }
+    #holder { position: relative; width: %s; height: %s; background: #111; display: block;}
+    h1 { %s padding: 20px; color: #fafafa; text-shadow: 0px 1px #000,0px -1px #000; position: absolute; font: 300 30px Roboto, Sans-serif;}
+    h2 { text-shadow: 0px 1px #000,0px -1px #000; font: 700 16px Roboto, Sans-serif;}
+    .meta {  position: absolute; opacity: 0.9; width: 220px; top: 80px; left: 20px; display: block; %s background: #000; line-height: 25px; color: #fafafa; border: 20px solid #000; font: 100 16px Roboto, Sans-serif;}
+    div.tooltip { position: absolute; width: 380px; display: block; %s padding: 20px; background: #000; border: 0px; border-radius: 3px; pointer-events: none; z-index: 999; color: #FAFAFA;}
+    }
+    </style>
+    <body>
+    <div id="holder">
+      <h1>%s</h1>
+      <p class="meta">
+      <b>Lens</b><br>%s<br><br>
+      <b>Length of intervals</b><br>%s<br><br>
+      <b>Overlap percentage</b><br>%s%%<br><br>
+      <b>Color Function</b><br>%s
+      </p>
+    </div>
+    <script src="https://cdnjs.cloudflare.com/ajax/libs/d3/3.5.5/d3.min.js"></script>
+    <script>
+    var width = %s,
+      height = %s;
+    var color = d3.scale.ordinal()
+      .domain(["0","1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13","14","15","16","17","18","19","20","21","22","23","24","25","26","27","28","29","30"])
+      .range(["#FF0000","#FF1400","#FF2800","#FF3c00","#FF5000","#FF6400","#FF7800","#FF8c00","#FFa000","#FFb400","#FFc800","#FFdc00","#FFf000","#fdff00","#b0ff00","#65ff00","#17ff00","#00ff36","#00ff83","#00ffd0","#00e4ff","#00c4ff","#00a4ff","#00a4ff","#0084ff","#0064ff","#0044ff","#0022ff","#0002ff","#0100ff","#0300ff","#0500ff"]);
+    var force = d3.layout.force()
+      .charge(%s)
+      .linkDistance(%s)
+      .gravity(%s)
+      .size([width, height]);
+    var svg = d3.select("#holder").append("svg")
+      .attr("width", width)
+      .attr("height", height);
+    
+    var div = d3.select("#holder").append("div")   
+      .attr("class", "tooltip")               
+      .style("opacity", 0.0);
+    
+    var divs = d3.select('#holder').append('div')
+      .attr('class', 'divs')
+      .attr('style', function(d) { return 'overflow: hidden; width: ' + width + 'px; height: ' + height + 'px;'; });  
+    
+      graph = %s;
+      force
+        .nodes(graph.nodes)
+        .links(graph.links)
+        .start();
+      var link = svg.selectAll(".link")
+        .data(graph.links)
+        .enter().append("line")
+        .attr("class", "link")
+        .style("stroke-width", function(d) { return Math.sqrt(d.value); });
+      var node = divs.selectAll('div')
+      .data(graph.nodes)
+        .enter().append('div')
+        .on("mouseover", function(d) {      
+          div.transition()        
+            .duration(200)      
+            .style("opacity", .9);
+          div .html(d.tooltip + "<br/>")  
+            .style("left", (d3.event.pageX + 100) + "px")     
+            .style("top", (d3.event.pageY - 28) + "px");    
+          })                  
+        .on("mouseout", function(d) {       
+          div.transition()        
+            .duration(500)      
+            .style("opacity", 0);   
+        })
+        .call(force.drag);
+      
+      node.append("title")
+        .text(function(d) { return d.name; });
+      force.on("tick", function() {
+      link.attr("x1", function(d) { return d.source.x; })
+        .attr("y1", function(d) { return d.source.y; })
+        .attr("x2", function(d) { return d.target.x; })
+        .attr("y2", function(d) { return d.target.y; });
+      node.attr("cx", function(d) { return d.x; })
+        .attr("cy", function(d) { return d.y; })
+        .attr('style', function(d) { return 'width: ' + (d.group * 2) + 'px; height: ' + (d.group * 2) + 'px; ' + 'left: '+(d.x-(d.group))+'px; ' + 'top: '+(d.y-(d.group))+'px; background: '+color(d.color)+'; box-shadow: 0px 0px 3px #111; box-shadow: 0px 0px 33px '+color(d.color)+', inset 0px 0px 5px rgba(0, 0, 0, 0.2);'})
+        ;
+      });
+    </script>"""%(title,width_css, height_css, title_display, meta_display, tooltips_display, title,complex["meta"],res,gain*100,color_function,width_js,height_js,graph_charge,graph_link_distance,graph_gravity,json.dumps(json_s))
+      outfile.write(html.encode("utf-8"))
+    if self.verbose > 0:
+      print("\nWrote d3.js graph to '%s'"%path_html)
diff --git a/utilities/Nerve_GIC/km.py.COPYRIGHT b/utilities/Nerve_GIC/km.py.COPYRIGHT
new file mode 100644
index 00000000..bef7b121
--- /dev/null
+++ b/utilities/Nerve_GIC/km.py.COPYRIGHT
@@ -0,0 +1,26 @@
+km.py is a fork of https://github.com/MLWave/kepler-mapper.
+Only the visualization part has been kept (Mapper part has been removed).
+
+This file has te following Copyright :
+
+The MIT License (MIT)
+
+Copyright (c) 2015 Triskelion - HJ van Veen - info@mlwave.com
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
diff --git a/utilities/Persistence_representations/CMakeLists.txt b/utilities/Persistence_representations/CMakeLists.txt
new file mode 100644
index 00000000..137eb0c1
--- /dev/null
+++ b/utilities/Persistence_representations/CMakeLists.txt
@@ -0,0 +1,53 @@
+# Copy files, otherwise result files are created in sources
+file(COPY "${CMAKE_SOURCE_DIR}/data/persistence_diagram/first.pers" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/")
+file(COPY "${CMAKE_SOURCE_DIR}/data/persistence_diagram/second.pers" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/")
+
+function(add_persistence_representation_creation_utility creation_utility)
+    add_executable ( ${creation_utility} ${creation_utility}.cpp )
+
+    # as the function is called in a subdirectory level, need to '../' to find persistence files
+    # ARGN will add all the other arguments (except creation_utility) sent to the CMake functions
+    add_test(NAME Persistence_representation_utilities_${creation_utility} COMMAND $<TARGET_FILE:${creation_utility}>
+        ${ARGN} "${CMAKE_CURRENT_BINARY_DIR}/../first.pers"
+        "${CMAKE_CURRENT_BINARY_DIR}/../second.pers")
+endfunction(add_persistence_representation_creation_utility)
+
+function(add_persistence_representation_plot_utility plot_utility tool_extension)
+    add_executable ( ${plot_utility} ${plot_utility}.cpp )
+
+    # as the function is called in a subdirectory level, need to '../' to find persistence heat maps files
+    add_test(NAME Persistence_representation_utilities_${plot_utility}_first COMMAND $<TARGET_FILE:${plot_utility}>
+        "${CMAKE_CURRENT_BINARY_DIR}/../first.pers${tool_extension}")
+    #add_test(NAME Persistence_representation_utilities_${plot_utility}_second COMMAND $<TARGET_FILE:${plot_utility}>
+    #    "${CMAKE_CURRENT_BINARY_DIR}/../second.pers${tool_extension}")
+    if(GNUPLOT_PATH)
+        add_test(NAME Persistence_representation_utilities_${plot_utility}_first_gnuplot COMMAND ${GNUPLOT_PATH}
+            "-e" "load '${CMAKE_CURRENT_BINARY_DIR}/../first.pers${tool_extension}_GnuplotScript'")
+        #add_test(NAME Persistence_representation_utilities_${plot_utility}_second_gnuplot COMMAND ${GNUPLOT_PATH}
+        #    "-e" "load '${CMAKE_CURRENT_BINARY_DIR}/../second.pers${tool_extension}_GnuplotScript'")
+    endif()
+endfunction(add_persistence_representation_plot_utility)
+
+function(add_persistence_representation_function_utility function_utility tool_extension)
+    add_executable ( ${function_utility} ${function_utility}.cpp )
+
+    # ARGV2 is an optional argument
+    if (${ARGV2})
+        # as the function is called in a subdirectory level, need to '../' to find persistence heat maps files
+        add_test(NAME Persistence_representation_utilities_${function_utility} COMMAND $<TARGET_FILE:${function_utility}>
+            "${ARGV2}"
+            "${CMAKE_CURRENT_BINARY_DIR}/../first.pers${tool_extension}"
+            "${CMAKE_CURRENT_BINARY_DIR}/../second.pers${tool_extension}")
+    else()
+        # as the function is called in a subdirectory level, need to '../' to find persistence heat maps files
+        add_test(NAME Persistence_representation_utilities_${function_utility} COMMAND $<TARGET_FILE:${function_utility}>
+            "${CMAKE_CURRENT_BINARY_DIR}/../first.pers${tool_extension}"
+            "${CMAKE_CURRENT_BINARY_DIR}/../second.pers${tool_extension}")
+    endif()
+endfunction(add_persistence_representation_function_utility)
+
+add_subdirectory(persistence_heat_maps)
+add_subdirectory(persistence_intervals)
+add_subdirectory(persistence_landscapes)
+add_subdirectory(persistence_landscapes_on_grid)
+add_subdirectory(persistence_vectors)
diff --git a/utilities/Persistence_representations/persistence_heat_maps/CMakeLists.txt b/utilities/Persistence_representations/persistence_heat_maps/CMakeLists.txt
new file mode 100644
index 00000000..386e9fa5
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_heat_maps/CMakeLists.txt
@@ -0,0 +1,15 @@
+cmake_minimum_required(VERSION 2.6)
+project(Persistence_representations_heat_maps_utilities)
+
+add_persistence_representation_creation_utility(create_pssk  "10" "-1" "-1" "4" "-1")
+add_persistence_representation_creation_utility(create_p_h_m_weighted_by_arctan_of_their_persistence "10" "-1" "-1" "4" "-1")
+add_persistence_representation_creation_utility(create_p_h_m_weighted_by_distance_from_diagonal "10" "-1" "-1" "4" "-1")
+add_persistence_representation_creation_utility(create_p_h_m_weighted_by_squared_diag_distance "10" "-1" "-1" "4" "-1")
+# Need to set grid min and max for further average, distance and scalar_product
+add_persistence_representation_creation_utility(create_persistence_heat_maps "10" "0" "35" "10" "-1")
+
+add_persistence_representation_plot_utility(plot_persistence_heat_map ".mps")
+
+add_persistence_representation_function_utility(average_persistence_heat_maps ".mps")
+add_persistence_representation_function_utility(compute_distance_of_persistence_heat_maps ".mps" "1")
+add_persistence_representation_function_utility(compute_scalar_product_of_persistence_heat_maps ".mps")
diff --git a/utilities/Persistence_representations/persistence_heat_maps/average_persistence_heat_maps.cpp b/utilities/Persistence_representations/persistence_heat_maps/average_persistence_heat_maps.cpp
new file mode 100644
index 00000000..6739e0b6
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_heat_maps/average_persistence_heat_maps.cpp
@@ -0,0 +1,63 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_heat_maps.h>
+
+#include <iostream>
+#include <vector>
+
+using constant_scaling_function = Gudhi::Persistence_representations::constant_scaling_function;
+using Persistence_heat_maps = Gudhi::Persistence_representations::Persistence_heat_maps<constant_scaling_function>;
+
+int main(int argc, char** argv) {
+  std::cout << "This program computes average of persistence heat maps stored in files (the files needs to be "
+            << "created beforehand).\n"
+            << "The parameters of this programs are names of files with persistence heat maps.\n";
+
+  if (argc < 3) {
+    std::cout << "Wrong number of parameters, the program will now terminate \n";
+    return 1;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 1; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+
+  std::vector<Persistence_heat_maps*> maps;
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    Persistence_heat_maps* l = new Persistence_heat_maps;
+    l->load_from_file(filenames[i]);
+    maps.push_back(l);
+  }
+
+  Persistence_heat_maps av;
+  av.compute_average(maps);
+  av.print_to_file("average.mps");
+
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    delete maps[i];
+  }
+
+  std::cout << "Average can be found in 'average.mps' file\n";
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_heat_maps/compute_distance_of_persistence_heat_maps.cpp b/utilities/Persistence_representations/persistence_heat_maps/compute_distance_of_persistence_heat_maps.cpp
new file mode 100644
index 00000000..ed8278a2
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_heat_maps/compute_distance_of_persistence_heat_maps.cpp
@@ -0,0 +1,94 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_heat_maps.h>
+
+#include <iostream>
+#include <sstream>
+#include <limits>
+#include <vector>
+
+using constant_scaling_function = Gudhi::Persistence_representations::constant_scaling_function;
+using Persistence_heat_maps = Gudhi::Persistence_representations::Persistence_heat_maps<constant_scaling_function>;
+
+int main(int argc, char** argv) {
+  std::cout << "This program computes distance of persistence heat maps stored in files (the files needs to be "
+            << "created beforehand).\n"
+            << "The first parameter of a program is an integer p. The program compute L^p distance of the two heat "
+            << "maps. For L^infty distance choose p = -1. \n"
+            << "The remaining parameters of this program are names of files with persistence heat maps.\n";
+
+  if (argc < 3) {
+    std::cout << "Wrong number of parameters, the program will now terminate \n";
+    return 1;
+  }
+
+  int pp = atoi(argv[1]);
+  double p = std::numeric_limits<double>::max();
+  if (pp != -1) {
+    p = pp;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 2; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+  std::vector<Persistence_heat_maps> maps;
+  maps.reserve(filenames.size());
+  for (size_t file_no = 0; file_no != filenames.size(); ++file_no) {
+    Persistence_heat_maps l;
+    l.load_from_file(filenames[file_no]);
+    maps.push_back(l);
+  }
+
+  // and now we will compute the scalar product of landscapes.
+
+  // first we prepare an array:
+  std::vector<std::vector<double> > distance(filenames.size());
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    std::vector<double> v(filenames.size(), 0);
+    distance[i] = v;
+  }
+
+  // and now we can compute the distances:
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    for (size_t j = i; j != filenames.size(); ++j) {
+      distance[i][j] = distance[j][i] = maps[i].distance(maps[j], p);
+    }
+  }
+
+  // and now output the result to the screen and a file:
+  std::ofstream out;
+  out.open("distance.mps");
+  for (size_t i = 0; i != distance.size(); ++i) {
+    for (size_t j = 0; j != distance.size(); ++j) {
+      std::cout << distance[i][j] << " ";
+      out << distance[i][j] << " ";
+    }
+    std::cout << std::endl;
+    out << std::endl;
+  }
+  out.close();
+
+  std::cout << "Distance can be found in 'distance.mps' file\n";
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_heat_maps/compute_scalar_product_of_persistence_heat_maps.cpp b/utilities/Persistence_representations/persistence_heat_maps/compute_scalar_product_of_persistence_heat_maps.cpp
new file mode 100644
index 00000000..63626853
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_heat_maps/compute_scalar_product_of_persistence_heat_maps.cpp
@@ -0,0 +1,85 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_heat_maps.h>
+
+#include <iostream>
+#include <sstream>
+#include <vector>
+
+using constant_scaling_function = Gudhi::Persistence_representations::constant_scaling_function;
+using Persistence_heat_maps = Gudhi::Persistence_representations::Persistence_heat_maps<constant_scaling_function>;
+
+int main(int argc, char** argv) {
+  std::cout << "This program computes scalar product of persistence heat maps stored in a file (the file needs to be "
+            << "created beforehand). \n"
+            << "The parameters of this programs are names of files with persistence heat maps.\n";
+
+  if (argc < 3) {
+    std::cout << "Wrong number of parameters, the program will now terminate \n";
+    return 1;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 1; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+  std::vector<Persistence_heat_maps> maps;
+  maps.reserve(filenames.size());
+  for (size_t file_no = 0; file_no != filenames.size(); ++file_no) {
+    Persistence_heat_maps l;
+    l.load_from_file(filenames[file_no]);
+    maps.push_back(l);
+  }
+
+  // and now we will compute the scalar product of landscapes.
+
+  // first we prepare an array:
+  std::vector<std::vector<double> > scalar_product(filenames.size());
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    std::vector<double> v(filenames.size(), 0);
+    scalar_product[i] = v;
+  }
+
+  // and now we can compute the scalar product:
+  for (size_t i = 0; i != maps.size(); ++i) {
+    for (size_t j = i; j != maps.size(); ++j) {
+      scalar_product[i][j] = scalar_product[j][i] = maps[i].compute_scalar_product(maps[j]);
+    }
+  }
+
+  // and now output the result to the screen and a file:
+  std::ofstream out;
+  out.open("scalar_product.mps");
+  for (size_t i = 0; i != scalar_product.size(); ++i) {
+    for (size_t j = 0; j != scalar_product.size(); ++j) {
+      std::cout << scalar_product[i][j] << " ";
+      out << scalar_product[i][j] << " ";
+    }
+    std::cout << std::endl;
+    out << std::endl;
+  }
+  out.close();
+
+  std::cout << "Distance can be found in 'scalar_product.mps' file\n";
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_arctan_of_their_persistence.cpp b/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_arctan_of_their_persistence.cpp
new file mode 100644
index 00000000..b4a1daa5
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_arctan_of_their_persistence.cpp
@@ -0,0 +1,81 @@
+/*    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):       Pawel Dlotko
+*
+*    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_heat_maps.h>
+
+#include <iostream>
+#include <sstream>
+#include <limits>
+#include <vector>
+
+using arc_tan_of_persistence_of_point = Gudhi::Persistence_representations::arc_tan_of_persistence_of_point;
+using Persistence_heat_maps =
+    Gudhi::Persistence_representations::Persistence_heat_maps<arc_tan_of_persistence_of_point>;
+
+int main(int argc, char** argv) {
+  std::cout << "This program creates persistence heat map files (*.mps) of persistence diagrams files (*.pers) "
+            << "provided as an input.The Gaussian kernels are weighted by the arc tangential of their persistence.\n"
+            << "The first parameter of a program is an integer, a size of a grid.\n"
+            << "The second and third parameters are min and max of the grid. If you want those numbers to be computed "
+            << "based on the data, set them both to -1 \n"
+            << "The fourth parameter is an integer, the standard deviation of a Gaussian kernel expressed in a number "
+            << "of pixels.\n"
+            << "The fifth parameter of this program is a dimension of persistence that will be used in creation of "
+            << "the persistence heat maps."
+            << "If your input files contains persistence pairs of various dimension, as a fifth parameter of the "
+            << "procedure please provide the dimension of persistence you want to use."
+            << "If in your files there are only birth-death pairs of the same dimension, set the fifth parameter to "
+            << "-1.\n"
+            << "The remaining parameters are the names of files with persistence diagrams. \n";
+
+  if (argc < 7) {
+    std::cout << "Wrong parameter list, the program will now terminate \n";
+    return 1;
+  }
+
+  size_t size_of_grid = (size_t)atoi(argv[1]);
+  double min_ = atof(argv[2]);
+  double max_ = atof(argv[3]);
+  size_t stdiv = atof(argv[4]);
+
+  unsigned dimension = std::numeric_limits<unsigned>::max();
+  int dim = atoi(argv[5]);
+  if (dim >= 0) {
+    dimension = (unsigned)dim;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 6; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+
+  std::vector<std::vector<double> > filter = Gudhi::Persistence_representations::create_Gaussian_filter(stdiv, 1);
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    std::cout << "Creating a heat map based on a file : " << filenames[i] << std::endl;
+    Persistence_heat_maps l(filenames[i], filter, false, size_of_grid, min_, max_, dimension);
+
+    std::stringstream ss;
+    ss << filenames[i] << ".mps";
+    l.print_to_file(ss.str().c_str());
+  }
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_distance_from_diagonal.cpp b/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_distance_from_diagonal.cpp
new file mode 100644
index 00000000..c50f9ddb
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_distance_from_diagonal.cpp
@@ -0,0 +1,81 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_heat_maps.h>
+
+#include <iostream>
+#include <sstream>
+#include <limits>
+#include <vector>
+
+using distance_from_diagonal_scaling = Gudhi::Persistence_representations::distance_from_diagonal_scaling;
+using Persistence_heat_maps = Gudhi::Persistence_representations::Persistence_heat_maps<distance_from_diagonal_scaling>;
+
+int main(int argc, char** argv) {
+  std::cout << "This program creates persistence heat map files (*.mps) of persistence diagrams files (*.pers) "
+            << "provided as an input.The Gaussian kernels are weighted by the distance of a center from the "
+            << "diagonal.\n"
+            << "The first parameter of a program is an integer, a size of a grid.\n"
+            << "The second and third parameters are min and max of the grid. If you want those numbers to be computed "
+            << "based on the data, set them both to -1 \n"
+            << "The fourth parameter is an integer, the standard deviation of a Gaussian kernel expressed in a number "
+            << "of pixels.\n"
+            << "The fifth parameter of this program is a dimension of persistence that will be used in creation of "
+            << "the persistence heat maps."
+            << "If your input files contains persistence pairs of various dimension, as a fifth parameter of the "
+            << "procedure please provide the dimension of persistence you want to use."
+            << "If in your files there are only birth-death pairs of the same dimension, set the fifth parameter to "
+            << "-1.\n"
+            << "The remaining parameters are the names of files with persistence diagrams. \n";
+
+  if (argc < 7) {
+    std::cout << "Wrong parameter list, the program will now terminate \n";
+    return 1;
+  }
+
+  size_t size_of_grid = (size_t)atoi(argv[1]);
+  double min_ = atof(argv[2]);
+  double max_ = atof(argv[3]);
+  size_t stdiv = atof(argv[4]);
+
+  unsigned dimension = std::numeric_limits<unsigned>::max();
+  int dim = atoi(argv[5]);
+  if (dim >= 0) {
+    dimension = (unsigned)dim;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 6; i != argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+
+  std::vector<std::vector<double> > filter = Gudhi::Persistence_representations::create_Gaussian_filter(stdiv, 1);
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    std::cout << "Creating a heat map based on a file : " << filenames[i] << std::endl;
+    Persistence_heat_maps l(filenames[i], filter, false, size_of_grid, min_, max_, dimension);
+
+    std::stringstream ss;
+    ss << filenames[i] << ".mps";
+    l.print_to_file(ss.str().c_str());
+  }
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_squared_diag_distance.cpp b/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_squared_diag_distance.cpp
new file mode 100644
index 00000000..59ff3c24
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_squared_diag_distance.cpp
@@ -0,0 +1,83 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_heat_maps.h>
+
+#include <iostream>
+#include <sstream>
+#include <limits>
+#include <vector>
+
+using squared_distance_from_diagonal_scaling =
+    Gudhi::Persistence_representations::squared_distance_from_diagonal_scaling;
+using Persistence_heat_maps =
+    Gudhi::Persistence_representations::Persistence_heat_maps<squared_distance_from_diagonal_scaling>;
+
+int main(int argc, char** argv) {
+  std::cout << "This program creates persistence heat map files (*.mps) of persistence diagrams files (*.pers) "
+            << "provided as an input.The Gaussian kernels are weighted by the square of distance of a center from the "
+            << "diagonal.\n"
+            << "The first parameter of a program is an integer, a size of a grid.\n"
+            << "The second and third parameters are min and max of the grid. If you want those numbers to be computed "
+            << "based on the data, set them both to -1 \n"
+            << "The fourth parameter is an integer, the standard deviation of a Gaussian kernel expressed in a number "
+            << "of pixels.\n"
+            << "The fifth parameter of this program is a dimension of persistence that will be used in creation of "
+            << "the persistence heat maps."
+            << "If your input files contains persistence pairs of various dimension, as a fifth parameter of the "
+            << "procedure please provide the dimension of persistence you want to use."
+            << "If in your files there are only birth-death pairs of the same dimension, set the fifth parameter to "
+            << "-1.\n"
+            << "The remaining parameters are the names of files with persistence diagrams. \n";
+
+  if (argc < 7) {
+    std::cout << "Wrong parameter list, the program will now terminate \n";
+    return 1;
+  }
+
+  size_t size_of_grid = (size_t)atoi(argv[1]);
+  double min_ = atof(argv[2]);
+  double max_ = atof(argv[3]);
+  size_t stdiv = atof(argv[4]);
+
+  unsigned dimension = std::numeric_limits<unsigned>::max();
+  int dim = atoi(argv[5]);
+  if (dim >= 0) {
+    dimension = (unsigned)dim;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 6; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+
+  std::vector<std::vector<double> > filter = Gudhi::Persistence_representations::create_Gaussian_filter(stdiv, 1);
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    std::cout << "Creating a heat map based on a file : " << filenames[i] << std::endl;
+    Persistence_heat_maps l(filenames[i], filter, false, size_of_grid, min_, max_, dimension);
+
+    std::stringstream ss;
+    ss << filenames[i] << ".mps";
+    l.print_to_file(ss.str().c_str());
+  }
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_heat_maps/create_persistence_heat_maps.cpp b/utilities/Persistence_representations/persistence_heat_maps/create_persistence_heat_maps.cpp
new file mode 100644
index 00000000..25cd1067
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_heat_maps/create_persistence_heat_maps.cpp
@@ -0,0 +1,78 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_heat_maps.h>
+
+#include <iostream>
+#include <sstream>
+#include <limits>
+#include <vector>
+
+using constant_scaling_function = Gudhi::Persistence_representations::constant_scaling_function;
+using Persistence_heat_maps = Gudhi::Persistence_representations::Persistence_heat_maps<constant_scaling_function>;
+
+int main(int argc, char** argv) {
+  std::cout << "This program creates persistence heat map files (*.mps) of persistence diagrams files (*.pers) "
+            << "provided as an input.\n"
+            << "The first parameter of a program is an integer, a size of a grid.\n"
+            << "The second and third parameters are min and max of the grid. If you want those numbers to be computed "
+            << "based on the data, set them both to -1 \n"
+            << "The fourth parameter is an integer, the standard deviation of a Gaussian kernel expressed in a number "
+            << "of pixels.\n"
+            << "The fifth parameter of this program is a dimension of persistence that will be used in creation of "
+            << "the persistence heat maps."
+            << "If your input files contains persistence pairs of various dimension, as a fifth parameter of the "
+            << "procedure please provide the dimension of persistence you want to use."
+            << "If in your files there are only birth-death pairs of the same dimension, set the fifth parameter to "
+            << "-1.\n"
+            << "The remaining parameters are the names of files with persistence diagrams. \n";
+
+  if (argc < 7) {
+    std::cout << "Wrong parameter list, the program will now terminate \n";
+    return 1;
+  }
+  size_t size_of_grid = (size_t)atoi(argv[1]);
+  double min_ = atof(argv[2]);
+  double max_ = atof(argv[3]);
+  size_t stdiv = atof(argv[4]);
+
+  unsigned dimension = std::numeric_limits<unsigned>::max();
+  int dim = atoi(argv[5]);
+  if (dim >= 0) {
+    dimension = (unsigned)dim;
+  }
+  std::vector<const char*> filenames;
+  for (int i = 6; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+
+  std::vector<std::vector<double> > filter = Gudhi::Persistence_representations::create_Gaussian_filter(stdiv, 1);
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    std::cout << "Creating a heat map based on file : " << filenames[i] << std::endl;
+    Persistence_heat_maps l(filenames[i], filter, false, size_of_grid, min_, max_, dimension);
+
+    std::stringstream ss;
+    ss << filenames[i] << ".mps";
+    l.print_to_file(ss.str().c_str());
+  }
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_heat_maps/create_pssk.cpp b/utilities/Persistence_representations/persistence_heat_maps/create_pssk.cpp
new file mode 100644
index 00000000..97ddb8f0
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_heat_maps/create_pssk.cpp
@@ -0,0 +1,79 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/PSSK.h>
+
+#include <iostream>
+#include <sstream>
+#include <limits>
+#include <vector>
+
+using PSSK = Gudhi::Persistence_representations::PSSK;
+
+int main(int argc, char** argv) {
+  std::cout << "This program creates PSSK files (*.pssk) of persistence diagrams files (*.pers) "
+            << "provided as an input.\n"
+            << "The first parameter of a program is an integer, a size of a grid.\n"
+            << "The second and third parameters are min and max of the grid. If you want those numbers to be computed "
+            << "based on the data, set them both to -1 \n"
+            << "The fourth parameter is an integer, the standard deviation of a Gaussian kernel expressed in a number "
+            << "of pixels.\n"
+            << "The fifth parameter of this program is a dimension of persistence that will be used in creation of "
+            << "the PSSK."
+            << "If your input files contains persistence pairs of various dimension, as a fifth parameter of the "
+            << "procedure please provide the dimension of persistence you want to use."
+            << "If in your files there are only birth-death pairs of the same dimension, set the first parameter to "
+            << "-1.\n"
+            << "The remaining parameters are the names of files with persistence diagrams. \n";
+
+  if (argc < 7) {
+    std::cout << "Wrong parameter list, the program will now terminate \n";
+    return 1;
+  }
+
+  size_t size_of_grid = (size_t)atoi(argv[1]);
+  double min_ = atof(argv[2]);
+  double max_ = atof(argv[3]);
+  size_t stdiv = atof(argv[4]);
+
+  unsigned dimension = std::numeric_limits<unsigned>::max();
+  int dim = atoi(argv[5]);
+  if (dim >= 0) {
+    dimension = (unsigned)dim;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 6; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+
+  std::vector<std::vector<double> > filter = Gudhi::Persistence_representations::create_Gaussian_filter(stdiv, 1);
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    std::cout << "Creating a PSSK based on a file : " << filenames[i] << std::endl;
+    PSSK l(filenames[i], filter, size_of_grid, min_, max_, dimension);
+
+    std::stringstream ss;
+    ss << filenames[i] << ".pssk";
+    l.print_to_file(ss.str().c_str());
+  }
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_heat_maps/plot_persistence_heat_map.cpp b/utilities/Persistence_representations/persistence_heat_maps/plot_persistence_heat_map.cpp
new file mode 100644
index 00000000..63711d83
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_heat_maps/plot_persistence_heat_map.cpp
@@ -0,0 +1,42 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_heat_maps.h>
+
+#include <iostream>
+#include <sstream>
+
+using constant_scaling_function = Gudhi::Persistence_representations::constant_scaling_function;
+using Persistence_heat_maps = Gudhi::Persistence_representations::Persistence_heat_maps<constant_scaling_function>;
+
+int main(int argc, char** argv) {
+  std::cout << "This program creates a gnuplot script from a persistence heat maps stored in a file (the file needs "
+            << "to be created beforehand). Please call the code with the name of a single heat maps file \n";
+  if (argc != 2) {
+    std::cout << "Wrong parameter list, the program will now terminate \n";
+    return 1;
+  }
+  Persistence_heat_maps l;
+  l.load_from_file(argv[1]);
+  l.plot(argv[1]);
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_intervals/CMakeLists.txt b/utilities/Persistence_representations/persistence_intervals/CMakeLists.txt
new file mode 100644
index 00000000..897e12a3
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_intervals/CMakeLists.txt
@@ -0,0 +1,32 @@
+cmake_minimum_required(VERSION 2.6)
+project(Persistence_representations_intervals_utilities)
+
+
+add_executable ( plot_histogram_of_intervals_lengths plot_histogram_of_intervals_lengths.cpp )
+
+add_test(NAME plot_histogram_of_intervals_lengths COMMAND $<TARGET_FILE:plot_histogram_of_intervals_lengths>
+    "${CMAKE_CURRENT_BINARY_DIR}/../first.pers" "-1")
+
+add_persistence_representation_plot_utility(plot_persistence_intervals "")
+add_persistence_representation_plot_utility(plot_persistence_Betti_numbers "")
+
+add_persistence_representation_creation_utility(compute_birth_death_range_in_persistence_diagram "-1")
+
+
+add_executable ( compute_number_of_dominant_intervals compute_number_of_dominant_intervals.cpp )
+add_test(NAME Persistence_representation_utilities_compute_number_of_dominant_intervals
+    COMMAND $<TARGET_FILE:compute_number_of_dominant_intervals>
+    "${CMAKE_CURRENT_BINARY_DIR}/../first.pers" "-1" "2")
+
+
+if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1)
+  add_executable ( compute_bottleneck_distance compute_bottleneck_distance.cpp )
+  if (TBB_FOUND)
+  target_link_libraries(compute_bottleneck_distance ${TBB_LIBRARIES})
+  endif(TBB_FOUND)
+  add_test(NAME Persistence_representation_utilities_compute_bottleneck_distance
+      COMMAND $<TARGET_FILE:compute_bottleneck_distance>
+      "-1"
+      "${CMAKE_CURRENT_BINARY_DIR}/../first.pers"
+      "${CMAKE_CURRENT_BINARY_DIR}/../second.pers")
+endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1)
diff --git a/utilities/Persistence_representations/persistence_intervals/compute_birth_death_range_in_persistence_diagram.cpp b/utilities/Persistence_representations/persistence_intervals/compute_birth_death_range_in_persistence_diagram.cpp
new file mode 100644
index 00000000..9102da79
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_intervals/compute_birth_death_range_in_persistence_diagram.cpp
@@ -0,0 +1,68 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_intervals.h>
+
+#include <iostream>
+#include <vector>
+#include <limits>
+#include <utility>
+
+using Persistence_intervals = Gudhi::Persistence_representations::Persistence_intervals;
+
+int main(int argc, char** argv) {
+  std::cout << "This program computes the range of birth and death times of persistence pairs in diagrams provided as "
+            << "an input.\n"
+            << "The first parameter is the dimension of persistence to be used to create persistence intervals. "
+            << "If your file contains the information about dimension of persistence pairs, please provide here the "
+            << "dimension of persistence pairs you want to use. "
+            << "If your input files consist only of birth-death pairs, please set this first parameter to -1.\n"
+            << "The remaining parameters of the program are the names of files with persistence diagrams.\n";
+
+  if (argc < 3) {
+    std::cout << "Wrong parameter list, the program will now terminate \n";
+    return 1;
+  }
+
+  unsigned dimension = std::numeric_limits<unsigned>::max();
+  int dim = atoi(argv[1]);
+  if (dim >= 0) {
+    dimension = (unsigned)dim;
+  }
+  std::vector<const char*> filenames;
+  for (int i = 2; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+
+  double min_ = std::numeric_limits<double>::max();
+  double max_ = -std::numeric_limits<double>::max();
+
+  for (size_t file_no = 0; file_no != filenames.size(); ++file_no) {
+    std::cout << "Creating diagram based on a file : " << filenames[file_no] << std::endl;
+    Persistence_intervals p(filenames[file_no], dimension);
+    std::pair<double, double> min_max_ = p.get_x_range();
+    if (min_max_.first < min_) min_ = min_max_.first;
+    if (min_max_.second > max_) max_ = min_max_.second;
+  }
+  std::cout << "Birth-death range : min: " << min_ << ", max: " << max_ << std::endl;
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_intervals/compute_bottleneck_distance.cpp b/utilities/Persistence_representations/persistence_intervals/compute_bottleneck_distance.cpp
new file mode 100644
index 00000000..c8290845
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_intervals/compute_bottleneck_distance.cpp
@@ -0,0 +1,95 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_intervals_with_distances.h>
+
+#include <iostream>
+#include <sstream>
+#include <limits>
+#include <vector>
+
+using Persistence_intervals_with_distances = Gudhi::Persistence_representations::Persistence_intervals_with_distances;
+
+int main(int argc, char** argv) {
+  std::cout << "This program computes the bottleneck distance of persistence pairs in diagrams provided as "
+            << "an input.\n"
+            << "The first parameter is the dimension of persistence to be used to create persistence intervals. "
+            << "If your file contains the information about dimension of persistence pairs, please provide here the "
+            << "dimension of persistence pairs you want to use. "
+            << "If your input files consist only of birth-death pairs, please set this first parameter to -1.\n"
+            << "The remaining parameters of the program are the names of files with persistence diagrams.\n";
+
+  if (argc < 3) {
+    std::cout << "Wrong number of parameters, the program will now terminate \n";
+    return 1;
+  }
+
+  unsigned dimension = std::numeric_limits<unsigned>::max();
+  int dim = atoi(argv[1]);
+  if (dim >= 0) {
+    dimension = (unsigned)dim;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 2; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+
+  // reading the persistence intervals:
+  std::vector<Persistence_intervals_with_distances> persistence_intervals;
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    Persistence_intervals_with_distances pers(filenames[i], dimension);
+    persistence_intervals.push_back(pers);
+  }
+
+  // and now we will compute the scalar product of landscapes.
+
+  // first we prepare an array:
+  std::vector<std::vector<double> > distance(filenames.size());
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    std::vector<double> v(filenames.size(), 0);
+    distance[i] = v;
+  }
+
+  // and now we can compute the distances:
+  for (size_t i = 0; i != persistence_intervals.size(); ++i) {
+    for (size_t j = i + 1; j != persistence_intervals.size(); ++j) {
+      distance[i][j] = distance[j][i] = persistence_intervals[i].distance(persistence_intervals[j]);
+    }
+  }
+
+  // and now output the result to the screen and a file:
+  std::ofstream out;
+  out.open("distance.itv");
+  for (size_t i = 0; i != distance.size(); ++i) {
+    for (size_t j = 0; j != distance.size(); ++j) {
+      std::cout << distance[i][j] << " ";
+      out << distance[i][j] << " ";
+    }
+    std::cout << std::endl;
+    out << std::endl;
+  }
+  out.close();
+
+  std::cout << "Distance can be found in 'distance.itv' file\n";
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_intervals/compute_number_of_dominant_intervals.cpp b/utilities/Persistence_representations/persistence_intervals/compute_number_of_dominant_intervals.cpp
new file mode 100644
index 00000000..b3d126f0
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_intervals/compute_number_of_dominant_intervals.cpp
@@ -0,0 +1,54 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_intervals.h>
+
+#include <iostream>
+#include <limits>
+#include <vector>
+#include <utility>
+
+using Persistence_intervals = Gudhi::Persistence_representations::Persistence_intervals;
+
+int main(int argc, char** argv) {
+  std::cout << "This program compute the dominant intervals. A number of intervals to be displayed is a parameter of "
+               "this program. \n";
+  if (argc != 4) {
+    std::cout << "To run this program, please provide the name of a file with persistence diagram, dimension of "
+                 "intervals that should be taken into account (if your file contains only persistence pairs in a "
+                 "single dimension, set it up to -1) and number of dominant intervals you would like to get \n";
+    return 1;
+  }
+  int dim = atoi(argv[2]);
+  unsigned dimension = std::numeric_limits<unsigned>::max();
+  if (dim >= 0) {
+    dimension = (unsigned)dim;
+  }
+  Persistence_intervals p(argv[1], dimension);
+  std::vector<std::pair<double, double> > dominant_intervals = p.dominant_intervals(atoi(argv[3]));
+  std::cout << "Here are the dominant intervals : " << std::endl;
+  for (size_t i = 0; i != dominant_intervals.size(); ++i) {
+    std::cout << " " << dominant_intervals[i].first << "," << dominant_intervals[i].second << " " << std::endl;
+  }
+
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_intervals/plot_histogram_of_intervals_lengths.cpp b/utilities/Persistence_representations/persistence_intervals/plot_histogram_of_intervals_lengths.cpp
new file mode 100644
index 00000000..ccb5b645
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_intervals/plot_histogram_of_intervals_lengths.cpp
@@ -0,0 +1,77 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_intervals.h>
+
+#include <iostream>
+#include <vector>
+#include <limits>
+#include <utility>
+
+using Persistence_intervals = Gudhi::Persistence_representations::Persistence_intervals;
+
+int main(int argc, char** argv) {
+  std::cout << "This program computes a histogram of barcode's length. A number of bins in the histogram is a "
+            << "parameter of this program. \n";
+  if ((argc != 3) && (argc != 4)) {
+    std::cout << "To run this program, please provide the name of a file with persistence diagram and number of "
+              << "dominant intervals you would like to get. Set a negative number dominant intervals value "
+              << "If your file contains only birth-death pairs.\n"
+              << "The third parameter is the dimension of the persistence that is to be used. If your "
+              << "file contains only birth-death pairs, you can skip this parameter\n";
+    return 1;
+  }
+
+  unsigned dominant_interval_number = std::numeric_limits<unsigned>::max();
+  int nbr = atoi(argv[2]);
+  if (nbr >= 0) {
+    dominant_interval_number = static_cast<unsigned>(nbr);
+  }
+
+  int persistence_dimension = -1;
+  if (argc == 4) {
+    persistence_dimension = atoi(argv[3]);
+  }
+
+  Persistence_intervals p(argv[1], dominant_interval_number);
+  std::vector<std::pair<double, double> > dominant_intervals = p.dominant_intervals(persistence_dimension);
+  std::vector<size_t> histogram = p.histogram_of_lengths(10);
+
+  std::stringstream gnuplot_script;
+  gnuplot_script << argv[1] << "_GnuplotScript";
+  std::ofstream out;
+  out.open(gnuplot_script.str().c_str());
+
+  out << "set style data histogram" << std::endl;
+  out << "set style histogram cluster gap 1" << std::endl;
+  out << "set style fill solid border -1" << std::endl;
+  out << "plot '-' notitle" << std::endl;
+  for (size_t i = 0; i != histogram.size(); ++i) {
+    out << histogram[i] << std::endl;
+  }
+  out << std::endl;
+  out.close();
+
+  std::cout << "To visualize, install gnuplot and type the command: gnuplot -persist -e \"load \'"
+            << gnuplot_script.str().c_str() << "\'\"" << std::endl;
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_intervals/plot_persistence_Betti_numbers.cpp b/utilities/Persistence_representations/persistence_intervals/plot_persistence_Betti_numbers.cpp
new file mode 100644
index 00000000..b433c2b3
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_intervals/plot_persistence_Betti_numbers.cpp
@@ -0,0 +1,87 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_intervals.h>
+
+#include <iostream>
+#include <vector>
+#include <limits>
+#include <utility>
+
+using Persistence_intervals = Gudhi::Persistence_representations::Persistence_intervals;
+
+int main(int argc, char** argv) {
+  if ((argc != 3) && (argc != 2)) {
+    std::cout << "This program creates a gnuplot script of Betti numbers from a single persistence diagram file"
+              << "(*.pers).\n"
+              << "To run this program, please provide the name of a file with persistence diagram.\n"
+              << "The second optional parameter of a program is the dimension of the persistence that is to be used. "
+              << "If your file contains only birth-death pairs, you can skip this parameter.\n";
+    return 1;
+  }
+
+  unsigned dimension = std::numeric_limits<unsigned>::max();
+  int dim = -1;
+  if (argc == 3) {
+    dim = atoi(argv[2]);
+  }
+  if (dim >= 0) {
+    dimension = (unsigned)dim;
+  }
+
+  Persistence_intervals p(argv[1], dimension);
+  std::vector<std::pair<double, size_t> > pbns = p.compute_persistent_betti_numbers();
+
+  // set up the ranges so that we see the image well.
+  double xRangeBegin = pbns[0].first;
+  double xRangeEnd = pbns[pbns.size() - 1].first;
+  double yRangeBegin = 0;
+  double yRangeEnd = 0;
+  for (size_t i = 0; i != pbns.size(); ++i) {
+    if (pbns[i].second > yRangeEnd) yRangeEnd = pbns[i].second;
+  }
+  xRangeBegin -= (xRangeEnd - xRangeBegin) / 100.0;
+  xRangeEnd += (xRangeEnd - xRangeBegin) / 100.0;
+  yRangeEnd += yRangeEnd / 100;
+
+  std::stringstream gnuplot_script;
+  gnuplot_script << argv[1] << "_GnuplotScript";
+  std::ofstream out;
+  out.open(gnuplot_script.str().c_str());
+
+  out << "set xrange [" << xRangeBegin << " : " << xRangeEnd << "]" << std::endl;
+  out << "set yrange [" << yRangeBegin << " : " << yRangeEnd << "]" << std::endl;
+  out << "plot '-' using 1:2 notitle with lp " << std::endl;
+  double previous_y = 0;
+  for (size_t i = 0; i != pbns.size(); ++i) {
+    out << pbns[i].first << " " << previous_y << std::endl;
+    out << pbns[i].first << " " << pbns[i].second << std::endl;
+    previous_y = pbns[i].second;
+  }
+  out << std::endl;
+  out.close();
+
+  std::cout << "To visualize, install gnuplot and type the command: gnuplot -persist -e \"load \'"
+            << gnuplot_script.str().c_str() << "\'\"" << std::endl;
+
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_intervals/plot_persistence_intervals.cpp b/utilities/Persistence_representations/persistence_intervals/plot_persistence_intervals.cpp
new file mode 100644
index 00000000..33387802
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_intervals/plot_persistence_intervals.cpp
@@ -0,0 +1,53 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_intervals.h>
+
+#include <iostream>
+#include <limits>
+#include <vector>
+#include <utility>
+
+using Persistence_intervals = Gudhi::Persistence_representations::Persistence_intervals;
+
+int main(int argc, char** argv) {
+  if ((argc != 3) && (argc != 2)) {
+    std::cout << "This program creates a gnuplot script from a single persistence diagram file (*.pers).\n"
+              << "To run this program, please provide the name of a file with persistence diagram.\n"
+              << "The second optional parameter of a program is the dimension of the persistence that is to be used. "
+              << "If your file contains only birth-death pairs, you can skip this parameter.\n";
+    return 1;
+  }
+  unsigned dimension = std::numeric_limits<unsigned>::max();
+  int dim = -1;
+  if (argc == 3) {
+    dim = atoi(argv[2]);
+  }
+  if (dim >= 0) {
+    dimension = (unsigned)dim;
+  }
+  std::vector<std::pair<double, double> > intervals =
+      Gudhi::Persistence_representations::read_persistence_intervals_in_one_dimension_from_file(argv[1], dimension);
+  Persistence_intervals b(intervals);
+  b.plot(argv[1]);
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_landscapes/CMakeLists.txt b/utilities/Persistence_representations/persistence_landscapes/CMakeLists.txt
new file mode 100644
index 00000000..d7087ed8
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_landscapes/CMakeLists.txt
@@ -0,0 +1,10 @@
+cmake_minimum_required(VERSION 2.6)
+project(Persistence_representations_landscapes_utilities)
+
+add_persistence_representation_creation_utility(create_landscapes "-1")
+
+add_persistence_representation_plot_utility(plot_landscapes ".land")
+
+add_persistence_representation_function_utility(average_landscapes ".land")
+add_persistence_representation_function_utility(compute_distance_of_landscapes ".land" "1")
+add_persistence_representation_function_utility(compute_scalar_product_of_landscapes ".land")
diff --git a/utilities/Persistence_representations/persistence_landscapes/average_landscapes.cpp b/utilities/Persistence_representations/persistence_landscapes/average_landscapes.cpp
new file mode 100644
index 00000000..1a59be8c
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_landscapes/average_landscapes.cpp
@@ -0,0 +1,63 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_landscape.h>
+
+#include <iostream>
+#include <vector>
+
+using Persistence_landscape = Gudhi::Persistence_representations::Persistence_landscape;
+
+int main(int argc, char** argv) {
+  std::cout << "This program computes average of persistence landscapes stored in files (the files needs to be "
+            << "created beforehand).\n"
+            << "The parameters of this programs are names of files with persistence landscapes.\n";
+  std::vector<const char*> filenames;
+
+  if (argc < 3) {
+    std::cout << "Wrong number of parameters, the program will now terminate \n";
+    return 1;
+  }
+
+  for (int i = 1; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+
+  std::vector<Persistence_landscape*> lands;
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    Persistence_landscape* l = new Persistence_landscape;
+    l->load_landscape_from_file(filenames[i]);
+    lands.push_back(l);
+  }
+
+  Persistence_landscape av;
+  av.compute_average(lands);
+
+  av.print_to_file("average.land");
+
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    delete lands[i];
+  }
+
+  std::cout << "Average can be found in 'average.land' file\n";
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_landscapes/compute_distance_of_landscapes.cpp b/utilities/Persistence_representations/persistence_landscapes/compute_distance_of_landscapes.cpp
new file mode 100644
index 00000000..5062f521
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_landscapes/compute_distance_of_landscapes.cpp
@@ -0,0 +1,93 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_landscape.h>
+
+#include <iostream>
+#include <sstream>
+#include <limits>
+#include <vector>
+
+using Persistence_landscape = Gudhi::Persistence_representations::Persistence_landscape;
+
+int main(int argc, char** argv) {
+  std::cout << "This program computes distance of persistence landscapes stored in files (the files needs to be "
+            << "created beforehand).\n"
+            << "The first parameter of a program is an integer p. The program compute L^p distance of the two heat "
+            << "maps. For L^infty distance choose p = -1. \n"
+            << "The remaining parameters of this program are names of files with persistence landscapes.\n";
+
+  if (argc < 3) {
+    std::cout << "Wrong number of parameters, the program will now terminate \n";
+    return 1;
+  }
+
+  int pp = atoi(argv[1]);
+  double p = std::numeric_limits<double>::max();
+  if (pp != -1) {
+    p = pp;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 2; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+  std::vector<Persistence_landscape> landscaspes;
+  landscaspes.reserve(filenames.size());
+  for (size_t file_no = 0; file_no != filenames.size(); ++file_no) {
+    Persistence_landscape l;
+    l.load_landscape_from_file(filenames[file_no]);
+    landscaspes.push_back(l);
+  }
+
+  // and now we will compute the scalar product of landscapes.
+
+  // first we prepare an array:
+  std::vector<std::vector<double> > distance(filenames.size());
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    std::vector<double> v(filenames.size(), 0);
+    distance[i] = v;
+  }
+
+  // and now we can compute the distances:
+  for (size_t i = 0; i != landscaspes.size(); ++i) {
+    for (size_t j = i; j != landscaspes.size(); ++j) {
+      distance[i][j] = distance[j][i] = compute_distance_of_landscapes(landscaspes[i], landscaspes[j], p);
+    }
+  }
+
+  // and now output the result to the screen and a file:
+  std::ofstream out;
+  out.open("distance.land");
+  for (size_t i = 0; i != distance.size(); ++i) {
+    for (size_t j = 0; j != distance.size(); ++j) {
+      std::cout << distance[i][j] << " ";
+      out << distance[i][j] << " ";
+    }
+    std::cout << std::endl;
+    out << std::endl;
+  }
+  out.close();
+
+  std::cout << "Distance can be found in 'distance.land' file\n";
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_landscapes/compute_scalar_product_of_landscapes.cpp b/utilities/Persistence_representations/persistence_landscapes/compute_scalar_product_of_landscapes.cpp
new file mode 100644
index 00000000..5b5e9fa3
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_landscapes/compute_scalar_product_of_landscapes.cpp
@@ -0,0 +1,84 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_landscape.h>
+
+#include <iostream>
+#include <sstream>
+#include <vector>
+
+using Persistence_landscape = Gudhi::Persistence_representations::Persistence_landscape;
+
+int main(int argc, char** argv) {
+  std::cout << "This program computes scalar product of persistence landscapes stored in a file (the file needs to be "
+            << "created beforehand). \n"
+            << "The parameters of this programs are names of files with persistence landscapes.\n";
+
+  if (argc < 3) {
+    std::cout << "Wrong number of parameters, the program will now terminate \n";
+    return 1;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 1; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+  std::vector<Persistence_landscape> landscaspes;
+  landscaspes.reserve(filenames.size());
+  for (size_t file_no = 0; file_no != filenames.size(); ++file_no) {
+    Persistence_landscape l;
+    l.load_landscape_from_file(filenames[file_no]);
+    landscaspes.push_back(l);
+  }
+
+  // and now we will compute the scalar product of landscapes.
+
+  // first we prepare an array:
+  std::vector<std::vector<double> > scalar_product(filenames.size());
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    std::vector<double> v(filenames.size(), 0);
+    scalar_product[i] = v;
+  }
+
+  // and now we can compute the scalar product:
+  for (size_t i = 0; i != landscaspes.size(); ++i) {
+    for (size_t j = i; j != landscaspes.size(); ++j) {
+      scalar_product[i][j] = scalar_product[j][i] = compute_inner_product(landscaspes[i], landscaspes[j]);
+    }
+  }
+
+  // and now output the result to the screen and a file:
+  std::ofstream out;
+  out.open("scalar_product.land");
+  for (size_t i = 0; i != scalar_product.size(); ++i) {
+    for (size_t j = 0; j != scalar_product.size(); ++j) {
+      std::cout << scalar_product[i][j] << " ";
+      out << scalar_product[i][j] << " ";
+    }
+    std::cout << std::endl;
+    out << std::endl;
+  }
+  out.close();
+
+  std::cout << "Distance can be found in 'scalar_product.land' file\n";
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_landscapes/create_landscapes.cpp b/utilities/Persistence_representations/persistence_landscapes/create_landscapes.cpp
new file mode 100644
index 00000000..6030e994
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_landscapes/create_landscapes.cpp
@@ -0,0 +1,65 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_landscape.h>
+
+#include <iostream>
+#include <sstream>
+#include <vector>
+#include <limits>
+
+using Persistence_landscape = Gudhi::Persistence_representations::Persistence_landscape;
+
+int main(int argc, char** argv) {
+  std::cout << "This program creates persistence landscapes files (*.land) of persistence diagrams files (*.pers) "
+            << "provided as an input.\n"
+            << "The first parameter of this program is a dimension of persistence that will be used in creation of "
+            << "the persistence heat maps."
+            << "If your input files contains persistence pairs of various dimension, as a first parameter of the "
+            << "procedure please provide the dimension of persistence you want to use."
+            << "If in your files there are only birth-death pairs of the same dimension, set the first parameter to "
+            << "-1.\n"
+            << "The remaining parameters are the names of files with persistence diagrams. \n";
+
+  if (argc < 3) {
+    std::cout << "Wrong parameter list, the program will now terminate \n";
+    return 1;
+  }
+  std::vector<const char*> filenames;
+  int dim = atoi(argv[1]);
+  unsigned dimension = std::numeric_limits<unsigned>::max();
+  if (dim >= 0) {
+    dimension = (unsigned)dim;
+  }
+  for (int i = 2; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    std::cout << "Creating a landscape based on file : " << filenames[i] << std::endl;
+    Persistence_landscape l(filenames[i], dimension);
+    std::stringstream ss;
+    ss << filenames[i] << ".land";
+    l.print_to_file(ss.str().c_str());
+  }
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_landscapes/plot_landscapes.cpp b/utilities/Persistence_representations/persistence_landscapes/plot_landscapes.cpp
new file mode 100644
index 00000000..c797a7a8
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_landscapes/plot_landscapes.cpp
@@ -0,0 +1,43 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_landscape.h>
+
+#include <iostream>
+#include <sstream>
+
+using Persistence_landscape = Gudhi::Persistence_representations::Persistence_landscape;
+
+int main(int argc, char** argv) {
+  std::cout << "This program creates a gnuplot script from a persistence landscape stored in a file (the file needs "
+            << "to be created beforehand). Please call the code with the name of a single landscape file.\n";
+  if (argc != 2) {
+    std::cout << "Wrong parameter list, the program will now terminate \n";
+    return 1;
+  }
+
+  Persistence_landscape l;
+  l.load_landscape_from_file(argv[1]);
+  l.plot(argv[1]);
+
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_landscapes_on_grid/CMakeLists.txt b/utilities/Persistence_representations/persistence_landscapes_on_grid/CMakeLists.txt
new file mode 100644
index 00000000..c5ea4bbf
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_landscapes_on_grid/CMakeLists.txt
@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.6)
+project(Persistence_representations_lanscapes_on_grid_utilities)
+
+# Need to set grid min and max for further average, distance and scalar_product
+add_persistence_representation_creation_utility(create_landscapes_on_grid "100" "0" "35" "-1")
+
+add_persistence_representation_plot_utility(plot_landscapes_on_grid ".g_land")
+
+add_persistence_representation_function_utility(average_landscapes_on_grid ".g_land")
+add_persistence_representation_function_utility(compute_distance_of_landscapes_on_grid ".g_land" "1")
+add_persistence_representation_function_utility(compute_scalar_product_of_landscapes_on_grid ".g_land")
diff --git a/utilities/Persistence_representations/persistence_landscapes_on_grid/average_landscapes_on_grid.cpp b/utilities/Persistence_representations/persistence_landscapes_on_grid/average_landscapes_on_grid.cpp
new file mode 100644
index 00000000..0b098d1a
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_landscapes_on_grid/average_landscapes_on_grid.cpp
@@ -0,0 +1,63 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_landscape_on_grid.h>
+
+#include <iostream>
+#include <vector>
+
+using Persistence_landscape_on_grid = Gudhi::Persistence_representations::Persistence_landscape_on_grid;
+
+int main(int argc, char** argv) {
+  std::cout << "This program computes average of persistence landscapes on grid stored in files (the files needs to "
+            << "be created beforehand).\n"
+            << "The parameters of this programs are names of files with persistence landscapes on grid.\n";
+
+  if (argc < 3) {
+    std::cout << "Wrong number of parameters, the program will now terminate \n";
+    return 1;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 1; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+
+  std::vector<Persistence_landscape_on_grid*> lands;
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    Persistence_landscape_on_grid* l = new Persistence_landscape_on_grid;
+    l->load_landscape_from_file(filenames[i]);
+    lands.push_back(l);
+  }
+
+  Persistence_landscape_on_grid av;
+  av.compute_average(lands);
+
+  av.print_to_file("average.g_land");
+
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    delete lands[i];
+  }
+
+  std::cout << "Average can be found in 'average.g_land' file\n";
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_landscapes_on_grid/compute_distance_of_landscapes_on_grid.cpp b/utilities/Persistence_representations/persistence_landscapes_on_grid/compute_distance_of_landscapes_on_grid.cpp
new file mode 100644
index 00000000..fd0fcd15
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_landscapes_on_grid/compute_distance_of_landscapes_on_grid.cpp
@@ -0,0 +1,93 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_landscape_on_grid.h>
+
+#include <iostream>
+#include <sstream>
+#include <limits>
+#include <vector>
+
+using Persistence_landscape_on_grid = Gudhi::Persistence_representations::Persistence_landscape_on_grid;
+
+int main(int argc, char** argv) {
+  std::cout << "This program computes distance of persistence landscapes on grid stored in files (the files needs to "
+            << "be created beforehand).\n"
+            << "The first parameter of a program is an integer p. The program compute L^p distance of the two heat "
+            << "maps. For L^infty distance choose p = -1. \n"
+            << "The remaining parameters of this program are names of files with persistence landscapes on grid.\n";
+
+  if (argc < 3) {
+    std::cout << "Wrong number of parameters, the program will now terminate \n";
+    return 1;
+  }
+
+  int pp = atoi(argv[1]);
+  double p = std::numeric_limits<double>::max();
+  if (pp != -1) {
+    p = pp;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 2; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+  std::vector<Persistence_landscape_on_grid> landscaspes;
+  landscaspes.reserve(filenames.size());
+  for (size_t file_no = 0; file_no != filenames.size(); ++file_no) {
+    Persistence_landscape_on_grid l;
+    l.load_landscape_from_file(filenames[file_no]);
+    landscaspes.push_back(l);
+  }
+
+  // and now we will compute the scalar product of landscapes.
+
+  // first we prepare an array:
+  std::vector<std::vector<double> > distance(filenames.size());
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    std::vector<double> v(filenames.size(), 0);
+    distance[i] = v;
+  }
+
+  // and now we can compute the scalar product:
+  for (size_t i = 0; i != landscaspes.size(); ++i) {
+    for (size_t j = i; j != landscaspes.size(); ++j) {
+      distance[i][j] = distance[j][i] = compute_distance_of_landscapes_on_grid(landscaspes[i], landscaspes[j], p);
+    }
+  }
+
+  // and now output the result to the screen and a file:
+  std::ofstream out;
+  out.open("distance.g_land");
+  for (size_t i = 0; i != distance.size(); ++i) {
+    for (size_t j = 0; j != distance.size(); ++j) {
+      std::cout << distance[i][j] << " ";
+      out << distance[i][j] << " ";
+    }
+    std::cout << std::endl;
+    out << std::endl;
+  }
+  out.close();
+
+  std::cout << "Distance can be found in 'distance.g_land' file\n";
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_landscapes_on_grid/compute_scalar_product_of_landscapes_on_grid.cpp b/utilities/Persistence_representations/persistence_landscapes_on_grid/compute_scalar_product_of_landscapes_on_grid.cpp
new file mode 100644
index 00000000..01de3dee
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_landscapes_on_grid/compute_scalar_product_of_landscapes_on_grid.cpp
@@ -0,0 +1,85 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_landscape_on_grid.h>
+
+#include <iostream>
+#include <sstream>
+#include <vector>
+
+using Persistence_landscape_on_grid = Gudhi::Persistence_representations::Persistence_landscape_on_grid;
+
+int main(int argc, char** argv) {
+  std::cout
+      << "This program computes scalar product of persistence landscapes on grid stored in a file (the file needs to "
+      << "be created beforehand). \n"
+      << "The parameters of this programs are names of files with persistence landscapes on grid.\n";
+
+  if (argc < 3) {
+    std::cout << "Wrong number of parameters, the program will now terminate \n";
+    return 1;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 1; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+  std::vector<Persistence_landscape_on_grid> landscaspes;
+  landscaspes.reserve(filenames.size());
+  for (size_t file_no = 0; file_no != filenames.size(); ++file_no) {
+    Persistence_landscape_on_grid l;
+    l.load_landscape_from_file(filenames[file_no]);
+    landscaspes.push_back(l);
+  }
+
+  // and now we will compute the scalar product of landscapes.
+
+  // first we prepare an array:
+  std::vector<std::vector<double> > scalar_product(filenames.size());
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    std::vector<double> v(filenames.size(), 0);
+    scalar_product[i] = v;
+  }
+
+  // and now we can compute the scalar product:
+  for (size_t i = 0; i != landscaspes.size(); ++i) {
+    for (size_t j = i; j != landscaspes.size(); ++j) {
+      scalar_product[i][j] = scalar_product[j][i] = compute_inner_product(landscaspes[i], landscaspes[j]);
+    }
+  }
+
+  // and now output the result to the screen and a file:
+  std::ofstream out;
+  out.open("scalar_product.g_land");
+  for (size_t i = 0; i != scalar_product.size(); ++i) {
+    for (size_t j = 0; j != scalar_product.size(); ++j) {
+      std::cout << scalar_product[i][j] << " ";
+      out << scalar_product[i][j] << " ";
+    }
+    std::cout << std::endl;
+    out << std::endl;
+  }
+  out.close();
+
+  std::cout << "Distance can be found in 'scalar_product.g_land' file\n";
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_landscapes_on_grid/create_landscapes_on_grid.cpp b/utilities/Persistence_representations/persistence_landscapes_on_grid/create_landscapes_on_grid.cpp
new file mode 100644
index 00000000..78e8ef57
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_landscapes_on_grid/create_landscapes_on_grid.cpp
@@ -0,0 +1,79 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_landscape_on_grid.h>
+
+#include <iostream>
+#include <sstream>
+#include <limits>
+#include <vector>
+
+using Persistence_landscape_on_grid = Gudhi::Persistence_representations::Persistence_landscape_on_grid;
+
+int main(int argc, char** argv) {
+  std::cout << "This program creates persistence landscapes on grid files (*.g_land) of persistence diagrams files "
+            << "(*.pers) provided as an input.\n"
+            << "The first parameter of a program is an integer, a size of a grid.\n"
+            << "The second and third parameters are min and max of the grid. If you want those numbers to be computed "
+            << "based on the data, set them both to -1 \n"
+            << "The fourth parameter of this program is a dimension of persistence that will be used in creation of "
+            << "the persistence heat maps."
+            << "If your input files contains persistence pairs of various dimension, as a fourth parameter of the "
+            << "procedure please provide the dimension of persistence you want to use."
+            << "If in your files there are only birth-death pairs of the same dimension, set the fourth parameter to "
+            << "-1.\n"
+            << "The remaining parameters are the names of files with persistence diagrams. \n";
+
+  if (argc < 6) {
+    std::cout << "Wrong parameter list, the program will now terminate \n";
+    return 1;
+  }
+
+  size_t size_of_grid = (size_t)atoi(argv[1]);
+  double min_ = atof(argv[2]);
+  double max_ = atof(argv[3]);
+  int dim = atoi(argv[4]);
+  unsigned dimension = std::numeric_limits<unsigned>::max();
+  if (dim >= 0) {
+    dimension = (unsigned)dim;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 5; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    std::cout << "Creating persistence landscape on a grid based on a file : " << filenames[i] << std::endl;
+    Persistence_landscape_on_grid l;
+    if ((min_ != -1) || (max_ != -1)) {
+      l = Persistence_landscape_on_grid(filenames[i], min_, max_, size_of_grid, dimension);
+    } else {
+      // (min_ == -1) && (max_ == -1), in this case the program will find min_ and max_ based on the data.
+      l = Persistence_landscape_on_grid(filenames[i], size_of_grid, dimension);
+    }
+    std::stringstream ss;
+    ss << filenames[i] << ".g_land";
+    l.print_to_file(ss.str().c_str());
+  }
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_landscapes_on_grid/plot_landscapes_on_grid.cpp b/utilities/Persistence_representations/persistence_landscapes_on_grid/plot_landscapes_on_grid.cpp
new file mode 100644
index 00000000..dddb3615
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_landscapes_on_grid/plot_landscapes_on_grid.cpp
@@ -0,0 +1,44 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_landscape_on_grid.h>
+
+#include <iostream>
+#include <sstream>
+
+using Persistence_landscape_on_grid = Gudhi::Persistence_representations::Persistence_landscape_on_grid;
+
+int main(int argc, char** argv) {
+  std::cout << "This program creates a gnuplot script from a persistence landscape on grid stored in a file (the file "
+            << "needs to be created beforehand). Please call the code with the name of a single landscape on grid file"
+            << ".\n";
+  if (argc != 2) {
+    std::cout << "Wrong parameter list, the program will now terminate \n";
+    return 1;
+  }
+
+  Persistence_landscape_on_grid l;
+  l.load_landscape_from_file(argv[1]);
+  l.plot(argv[1]);
+
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_vectors/CMakeLists.txt b/utilities/Persistence_representations/persistence_vectors/CMakeLists.txt
new file mode 100644
index 00000000..a401c955
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_vectors/CMakeLists.txt
@@ -0,0 +1,10 @@
+cmake_minimum_required(VERSION 2.6)
+project(Persistence_vectors_utilities)
+
+add_persistence_representation_creation_utility(create_persistence_vectors "-1")
+
+add_persistence_representation_plot_utility(plot_persistence_vectors ".vect")
+
+add_persistence_representation_function_utility(average_persistence_vectors ".vect")
+add_persistence_representation_function_utility(compute_distance_of_persistence_vectors ".vect" "1")
+add_persistence_representation_function_utility(compute_scalar_product_of_persistence_vectors ".vect")
diff --git a/utilities/Persistence_representations/persistence_vectors/average_persistence_vectors.cpp b/utilities/Persistence_representations/persistence_vectors/average_persistence_vectors.cpp
new file mode 100644
index 00000000..0144e76f
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_vectors/average_persistence_vectors.cpp
@@ -0,0 +1,65 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_vectors.h>
+
+#include <iostream>
+#include <vector>
+
+using Euclidean_distance = Gudhi::Euclidean_distance;
+using Vector_distances_in_diagram = Gudhi::Persistence_representations::Vector_distances_in_diagram<Euclidean_distance>;
+
+int main(int argc, char** argv) {
+  std::cout << "This program computes average of persistence vectors stored in files (the files needs to "
+            << "be created beforehand).\n"
+            << "The parameters of this programs are names of files with persistence vectors.\n";
+
+  if (argc < 3) {
+    std::cout << "Wrong number of parameters, the program will now terminate \n";
+    return 1;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 1; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+
+  std::vector<Vector_distances_in_diagram*> lands;
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    Vector_distances_in_diagram* l = new Vector_distances_in_diagram;
+    l->load_from_file(filenames[i]);
+    lands.push_back(l);
+  }
+
+  Vector_distances_in_diagram av;
+  av.compute_average(lands);
+
+  av.print_to_file("average.vect");
+
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    delete lands[i];
+  }
+
+  std::cout << "Done \n";
+
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_vectors/compute_distance_of_persistence_vectors.cpp b/utilities/Persistence_representations/persistence_vectors/compute_distance_of_persistence_vectors.cpp
new file mode 100644
index 00000000..7e66d25e
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_vectors/compute_distance_of_persistence_vectors.cpp
@@ -0,0 +1,94 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_vectors.h>
+
+#include <iostream>
+#include <sstream>
+#include <limits>
+#include <vector>
+
+using Euclidean_distance = Gudhi::Euclidean_distance;
+using Vector_distances_in_diagram = Gudhi::Persistence_representations::Vector_distances_in_diagram<Euclidean_distance>;
+
+int main(int argc, char** argv) {
+  std::cout << "This program compute distance of persistence vectors stored in a file (the file needs to be created "
+               "beforehand). \n";
+  std::cout << "The first parameter of a program is an integer p. The program compute l^p distance of the vectors. For "
+               "l^infty distance choose p = -1. \n";
+  std::cout << "The remaining parameters of this programs are names of files with persistence vectors.\n";
+
+  if (argc < 3) {
+    std::cout << "Wrong number of parameters, the program will now terminate \n";
+    return 1;
+  }
+
+  int pp = atoi(argv[1]);
+  double p = std::numeric_limits<double>::max();
+  if (pp != -1) {
+    p = pp;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 2; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+  std::vector<Vector_distances_in_diagram> vectors;
+  vectors.reserve(filenames.size());
+  for (size_t file_no = 0; file_no != filenames.size(); ++file_no) {
+    Vector_distances_in_diagram l;
+    l.load_from_file(filenames[file_no]);
+    vectors.push_back(l);
+  }
+
+  // and now we will compute the scalar product of landscapes.
+
+  // first we prepare an array:
+  std::vector<std::vector<double> > distance(filenames.size());
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    std::vector<double> v(filenames.size(), 0);
+    distance[i] = v;
+  }
+
+  // and now we can compute the distances:
+  for (size_t i = 0; i != vectors.size(); ++i) {
+    for (size_t j = i + 1; j != vectors.size(); ++j) {
+      distance[i][j] = distance[j][i] = vectors[i].distance(vectors[j], p);
+    }
+  }
+
+  // and now output the result to the screen and a file:
+  std::ofstream out;
+  out.open("distance.vect");
+  for (size_t i = 0; i != distance.size(); ++i) {
+    for (size_t j = 0; j != distance.size(); ++j) {
+      std::cout << distance[i][j] << " ";
+      out << distance[i][j] << " ";
+    }
+    std::cout << std::endl;
+    out << std::endl;
+  }
+  out.close();
+
+  std::cout << "Distance can be found in 'distance.vect' file\n";
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_vectors/compute_scalar_product_of_persistence_vectors.cpp b/utilities/Persistence_representations/persistence_vectors/compute_scalar_product_of_persistence_vectors.cpp
new file mode 100644
index 00000000..303c6e3e
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_vectors/compute_scalar_product_of_persistence_vectors.cpp
@@ -0,0 +1,86 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_vectors.h>
+
+#include <iostream>
+#include <sstream>
+#include <limits>
+#include <vector>
+
+using Euclidean_distance = Gudhi::Euclidean_distance;
+using Vector_distances_in_diagram = Gudhi::Persistence_representations::Vector_distances_in_diagram<Euclidean_distance>;
+
+int main(int argc, char** argv) {
+  std::cout << "This program computes scalar product of persistence vectors stored in a file (the file needs to "
+            << "be created beforehand). \n"
+            << "The parameters of this programs are names of files with persistence vectors.\n";
+
+  if (argc < 3) {
+    std::cout << "Wrong number of parameters, the program will now terminate \n";
+    return 1;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 1; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+  std::vector<Vector_distances_in_diagram> vectors;
+  vectors.reserve(filenames.size());
+  for (size_t file_no = 0; file_no != filenames.size(); ++file_no) {
+    Vector_distances_in_diagram l;
+    l.load_from_file(filenames[file_no]);
+    vectors.push_back(l);
+  }
+
+  // and now we will compute the scalar product of landscapes.
+
+  // first we prepare an array:
+  std::vector<std::vector<double> > scalar_product(filenames.size());
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    std::vector<double> v(filenames.size(), 0);
+    scalar_product[i] = v;
+  }
+
+  // and now we can compute the scalar product:
+  for (size_t i = 0; i != vectors.size(); ++i) {
+    for (size_t j = i; j != vectors.size(); ++j) {
+      scalar_product[i][j] = scalar_product[j][i] = vectors[i].compute_scalar_product(vectors[j]);
+    }
+  }
+
+  // and now output the result to the screen and a file:
+  std::ofstream out;
+  out.open("scalar_product.vect");
+  for (size_t i = 0; i != scalar_product.size(); ++i) {
+    for (size_t j = 0; j != scalar_product.size(); ++j) {
+      std::cout << scalar_product[i][j] << " ";
+      out << scalar_product[i][j] << " ";
+    }
+    std::cout << std::endl;
+    out << std::endl;
+  }
+  out.close();
+
+  std::cout << "Distance can be found in 'scalar_product.vect' file\n";
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_vectors/create_persistence_vectors.cpp b/utilities/Persistence_representations/persistence_vectors/create_persistence_vectors.cpp
new file mode 100644
index 00000000..cc5e5393
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_vectors/create_persistence_vectors.cpp
@@ -0,0 +1,69 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_vectors.h>
+
+#include <iostream>
+#include <sstream>
+#include <limits>
+#include <vector>
+
+using Euclidean_distance = Gudhi::Euclidean_distance;
+using Vector_distances_in_diagram = Gudhi::Persistence_representations::Vector_distances_in_diagram<Euclidean_distance>;
+
+int main(int argc, char** argv) {
+  std::cout << "This program creates persistence vectors files (*.vect) of persistence diagrams files (*.pers) "
+            << "provided as an input.\n"
+            << "The first parameter of this program is a dimension of persistence that will be used in creation of "
+            << "the persistence heat maps."
+            << "If your input files contains persistence pairs of various dimension, as a first parameter of the "
+            << "procedure please provide the dimension of persistence you want to use."
+            << "If in your files there are only birth-death pairs of the same dimension, set the first parameter to "
+            << "-1.\n"
+            << "The remaining parameters are the names of files with persistence diagrams. \n";
+
+  if (argc < 3) {
+    std::cout << "Wrong parameter list, the program will now terminate \n";
+    return 1;
+  }
+
+  std::cout << "The remaining parameters are the names of files with persistence diagrams. \n";
+  int dim = atoi(argv[1]);
+  unsigned dimension = std::numeric_limits<unsigned>::max();
+  if (dim >= 0) {
+    dimension = (unsigned)dim;
+  }
+
+  std::vector<const char*> filenames;
+  for (int i = 2; i < argc; ++i) {
+    filenames.push_back(argv[i]);
+  }
+
+  for (size_t i = 0; i != filenames.size(); ++i) {
+    std::cerr << "Creating persistence vectors based on a file : " << filenames[i] << std::endl;
+    Vector_distances_in_diagram l(filenames[i], dimension);
+    std::stringstream ss;
+    ss << filenames[i] << ".vect";
+    l.print_to_file(ss.str().c_str());
+  }
+  return 0;
+}
diff --git a/utilities/Persistence_representations/persistence_vectors/plot_persistence_vectors.cpp b/utilities/Persistence_representations/persistence_vectors/plot_persistence_vectors.cpp
new file mode 100644
index 00000000..aa33107d
--- /dev/null
+++ b/utilities/Persistence_representations/persistence_vectors/plot_persistence_vectors.cpp
@@ -0,0 +1,43 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Persistence_vectors.h>
+
+#include <iostream>
+#include <sstream>
+
+using Euclidean_distance = Gudhi::Euclidean_distance;
+using Vector_distances_in_diagram = Gudhi::Persistence_representations::Vector_distances_in_diagram<Euclidean_distance>;
+
+int main(int argc, char** argv) {
+  std::cout << "This program create a Gnuplot script to plot persistence vector. Please call this program with the "
+               "name of file with persistence vector. \n";
+  if (argc != 2) {
+    std::cout << "Wrong number of parameters, the program will now terminate. \n";
+    return 1;
+  }
+  Vector_distances_in_diagram l;
+  l.load_from_file(argv[1]);
+  l.plot(argv[1]);
+
+  return 0;
+}
diff --git a/utilities/Rips_complex/CMakeLists.txt b/utilities/Rips_complex/CMakeLists.txt
new file mode 100644
index 00000000..baa571fa
--- /dev/null
+++ b/utilities/Rips_complex/CMakeLists.txt
@@ -0,0 +1,21 @@
+cmake_minimum_required(VERSION 2.6)
+project(Rips_complex_utilities)
+
+add_executable(rips_distance_matrix_persistence rips_distance_matrix_persistence.cpp)
+target_link_libraries(rips_distance_matrix_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY})
+
+add_executable(rips_persistence rips_persistence.cpp)
+target_link_libraries(rips_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY})
+
+if (TBB_FOUND)
+  target_link_libraries(rips_distance_matrix_persistence ${TBB_LIBRARIES})
+  target_link_libraries(rips_persistence ${TBB_LIBRARIES})
+endif()
+
+add_test(NAME Rips_complex_utility_from_rips_distance_matrix COMMAND $<TARGET_FILE:rips_distance_matrix_persistence>
+    "${CMAKE_SOURCE_DIR}/data/distance_matrix/full_square_distance_matrix.csv" "-r" "1.0" "-d" "3" "-p" "3" "-m" "0")
+add_test(NAME Rips_complex_utility_from_rips_on_tore_3D COMMAND $<TARGET_FILE:rips_persistence>
+    "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3")
+
+install(TARGETS rips_distance_matrix_persistence DESTINATION bin)
+install(TARGETS rips_persistence DESTINATION bin)
diff --git a/utilities/Rips_complex/rips_distance_matrix_persistence.cpp b/utilities/Rips_complex/rips_distance_matrix_persistence.cpp
new file mode 100644
index 00000000..ca3c0327
--- /dev/null
+++ b/utilities/Rips_complex/rips_distance_matrix_persistence.cpp
@@ -0,0 +1,133 @@
+/*    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):       Pawel Dlotko, Vincent Rouvreau
+ *
+ *    Copyright (C) 2016  INRIA
+ *
+ *    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 <gudhi/Rips_complex.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/reader_utils.h>
+
+#include <boost/program_options.hpp>
+
+#include <string>
+#include <vector>
+#include <limits>  // infinity
+
+// Types definition
+using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
+using Filtration_value = Simplex_tree::Filtration_value;
+using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>;
+using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
+using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>;
+using Distance_matrix = std::vector<std::vector<Filtration_value>>;
+
+void program_options(int argc, char* argv[], std::string& csv_matrix_file, std::string& filediag,
+                     Filtration_value& threshold, int& dim_max, int& p, Filtration_value& min_persistence);
+
+int main(int argc, char* argv[]) {
+  std::string csv_matrix_file;
+  std::string filediag;
+  Filtration_value threshold;
+  int dim_max;
+  int p;
+  Filtration_value min_persistence;
+
+  program_options(argc, argv, csv_matrix_file, filediag, threshold, dim_max, p, min_persistence);
+
+  Distance_matrix distances = Gudhi::read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_matrix_file);
+  Rips_complex rips_complex_from_file(distances, threshold);
+
+  // Construct the Rips complex in a Simplex Tree
+  Simplex_tree simplex_tree;
+
+  rips_complex_from_file.create_complex(simplex_tree, dim_max);
+  std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n";
+  std::cout << "   and has dimension " << simplex_tree.dimension() << " \n";
+
+  // Sort the simplices in the order of the filtration
+  simplex_tree.initialize_filtration();
+
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(simplex_tree);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(p);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  // Output the diagram in filediag
+  if (filediag.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::ofstream out(filediag);
+    pcoh.output_diagram(out);
+    out.close();
+  }
+  return 0;
+}
+
+void program_options(int argc, char* argv[], std::string& csv_matrix_file, std::string& filediag,
+                     Filtration_value& threshold, int& dim_max, int& p, Filtration_value& min_persistence) {
+  namespace po = boost::program_options;
+  po::options_description hidden("Hidden options");
+  hidden.add_options()(
+      "input-file", po::value<std::string>(&csv_matrix_file),
+      "Name of file containing a distance matrix. Can be square or lower triangular matrix. Separator is ';'.");
+
+  po::options_description visible("Allowed options", 100);
+  visible.add_options()("help,h", "produce help message")(
+      "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()),
+      "Name of file in which the persistence diagram is written. Default print in std::cout")(
+      "max-edge-length,r",
+      po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()),
+      "Maximal length of an edge for the Rips complex construction.")(
+      "cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),
+      "Maximal dimension of the Rips complex we want to compute.")(
+      "field-charac,p", po::value<int>(&p)->default_value(11),
+      "Characteristic p of the coefficient field Z/pZ for computing homology.")(
+      "min-persistence,m", po::value<Filtration_value>(&min_persistence),
+      "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length "
+      "intervals");
+
+  po::positional_options_description pos;
+  pos.add("input-file", 1);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+
+  po::variables_map vm;
+  po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  if (vm.count("help") || !vm.count("input-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::cout << "of a Rips complex defined on a set of distance matrix.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
+  }
+}
diff --git a/utilities/Rips_complex/rips_persistence.cpp b/utilities/Rips_complex/rips_persistence.cpp
new file mode 100644
index 00000000..8405c014
--- /dev/null
+++ b/utilities/Rips_complex/rips_persistence.cpp
@@ -0,0 +1,135 @@
+/*    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):       Clément Maria
+ *
+ *    Copyright (C) 2014  INRIA
+ *
+ *    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 <gudhi/Rips_complex.h>
+#include <gudhi/distance_functions.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Points_off_io.h>
+
+#include <boost/program_options.hpp>
+
+#include <string>
+#include <vector>
+#include <limits>  // infinity
+
+// Types definition
+using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
+using Filtration_value = Simplex_tree::Filtration_value;
+using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>;
+using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
+using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>;
+using Point = std::vector<double>;
+using Points_off_reader = Gudhi::Points_off_reader<Point>;
+
+void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag,
+                     Filtration_value& threshold, int& dim_max, int& p, Filtration_value& min_persistence);
+
+int main(int argc, char* argv[]) {
+  std::string off_file_points;
+  std::string filediag;
+  Filtration_value threshold;
+  int dim_max;
+  int p;
+  Filtration_value min_persistence;
+
+  program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence);
+
+  Points_off_reader off_reader(off_file_points);
+  Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance());
+
+  // Construct the Rips complex in a Simplex Tree
+  Simplex_tree simplex_tree;
+
+  rips_complex_from_file.create_complex(simplex_tree, dim_max);
+  std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n";
+  std::cout << "   and has dimension " << simplex_tree.dimension() << " \n";
+
+  // Sort the simplices in the order of the filtration
+  simplex_tree.initialize_filtration();
+
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(simplex_tree);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(p);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  // Output the diagram in filediag
+  if (filediag.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::ofstream out(filediag);
+    pcoh.output_diagram(out);
+    out.close();
+  }
+
+  return 0;
+}
+
+void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag,
+                     Filtration_value& threshold, int& dim_max, int& p, Filtration_value& min_persistence) {
+  namespace po = boost::program_options;
+  po::options_description hidden("Hidden options");
+  hidden.add_options()("input-file", po::value<std::string>(&off_file_points),
+                       "Name of an OFF file containing a point set.\n");
+
+  po::options_description visible("Allowed options", 100);
+  visible.add_options()("help,h", "produce help message")(
+      "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()),
+      "Name of file in which the persistence diagram is written. Default print in std::cout")(
+      "max-edge-length,r",
+      po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()),
+      "Maximal length of an edge for the Rips complex construction.")(
+      "cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),
+      "Maximal dimension of the Rips complex we want to compute.")(
+      "field-charac,p", po::value<int>(&p)->default_value(11),
+      "Characteristic p of the coefficient field Z/pZ for computing homology.")(
+      "min-persistence,m", po::value<Filtration_value>(&min_persistence),
+      "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length "
+      "intervals");
+
+  po::positional_options_description pos;
+  pos.add("input-file", 1);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+
+  po::variables_map vm;
+  po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  if (vm.count("help") || !vm.count("input-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::cout << "of a Rips complex defined on a set of input points.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
+  }
+}
diff --git a/utilities/Rips_complex/ripscomplex.md b/utilities/Rips_complex/ripscomplex.md
new file mode 100644
index 00000000..4291fae7
--- /dev/null
+++ b/utilities/Rips_complex/ripscomplex.md
@@ -0,0 +1,49 @@
+
+
+# Rips complex #
+
+## rips_persistence ##
+This program computes the persistent homology with coefficient field *Z/pZ* of a Rips complex defined on a set of input points, using Euclidean distance. The output diagram contains one bar per line, written with the convention:
+
+`p dim birth death`
+
+where `dim` is the dimension of the homological feature, `birth` and `death` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p` must be a prime number).
+
+**Usage**
+
+`rips_persistence [options] <OFF input file>`
+
+**Allowed options**
+
+* `-h [ --help ]` Produce help message
+* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output.
+* `-r [ --max-edge-length ]` (default = inf) Maximal length of an edge for the Rips complex construction.
+* `-d [ --cpx-dimension ]` (default = 1) Maximal dimension of the Rips complex we want to compute.
+* `-p [ --field-charac ]` (default = 11)     Characteristic p of the coefficient field Z/pZ for computing homology.
+* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals.
+
+Beware: this program may use a lot of RAM and take a lot of time if `max-edge-length` is set to a large value.
+
+**Example 1 with Z/2Z coefficients**
+
+`rips_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 2`
+
+**Example 2 with Z/3Z coefficients**
+
+`rips_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 3`
+
+
+## rips_distance_matrix_persistence ##
+
+Same as `rips_persistence` but taking a distance matrix as input.
+
+**Usage**
+
+`rips_persistence [options] <CSV input file>`
+
+where
+`<CSV input file>` is the path to the file containing a distance matrix. Can be square or lower triangular matrix. Separator is ';'.
+
+**Example**
+
+`rips_distance_matrix_persistence data/distance_matrix/full_square_distance_matrix.csv -r 15 -d 3 -p 3 -m 0`
diff --git a/utilities/Witness_complex/CMakeLists.txt b/utilities/Witness_complex/CMakeLists.txt
new file mode 100644
index 00000000..125a41ff
--- /dev/null
+++ b/utilities/Witness_complex/CMakeLists.txt
@@ -0,0 +1,28 @@
+cmake_minimum_required(VERSION 2.6)
+project(Witness_complex_utilities)
+
+# CGAL and Eigen3 are required for Euclidean version of Witness
+if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.6.0)
+  
+  add_executable ( Witness_complex_strong_witness_persistence strong_witness_persistence.cpp )
+  target_link_libraries(Witness_complex_strong_witness_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY})
+
+  add_executable ( Witness_complex_weak_witness_persistence weak_witness_persistence.cpp )
+  target_link_libraries(Witness_complex_weak_witness_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY})
+
+  if (TBB_FOUND)
+    target_link_libraries(Witness_complex_strong_witness_persistence ${TBB_LIBRARIES})
+    target_link_libraries(Witness_complex_weak_witness_persistence ${TBB_LIBRARIES})
+  endif()
+
+  add_test(NAME Witness_complex_strong_test_torus_persistence
+      COMMAND $<TARGET_FILE:Witness_complex_strong_witness_persistence>
+      "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-l" "20" "-a" "0.5")
+  add_test(NAME Witness_complex_weak_test_torus_persistence
+      COMMAND $<TARGET_FILE:Witness_complex_weak_witness_persistence>
+      "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-l" "20" "-a" "0.5")
+
+  install(TARGETS Witness_complex_strong_witness_persistence DESTINATION bin)
+  install(TARGETS Witness_complex_weak_witness_persistence DESTINATION bin)
+
+endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.6.0)
diff --git a/utilities/Witness_complex/strong_witness_persistence.cpp b/utilities/Witness_complex/strong_witness_persistence.cpp
new file mode 100644
index 00000000..2fba631b
--- /dev/null
+++ b/utilities/Witness_complex/strong_witness_persistence.cpp
@@ -0,0 +1,156 @@
+/*    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):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Simplex_tree.h>
+#include <gudhi/Euclidean_strong_witness_complex.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Points_off_io.h>
+#include <gudhi/pick_n_random_points.h>
+#include <gudhi/choose_n_farthest_points.h>
+
+#include <boost/program_options.hpp>
+
+#include <CGAL/Epick_d.h>
+
+#include <string>
+#include <vector>
+#include <limits>  // infinity
+
+using K = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
+using Point_d = K::Point_d;
+
+using Point_vector = std::vector<Point_d>;
+using Strong_witness_complex = Gudhi::witness_complex::Euclidean_strong_witness_complex<K>;
+using SimplexTree = Gudhi::Simplex_tree<>;
+
+using Filtration_value = SimplexTree::Filtration_value;
+
+using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
+using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<SimplexTree, Field_Zp>;
+
+void program_options(int argc, char* argv[], int& nbL, std::string& file_name, std::string& filediag,
+                     Filtration_value& max_squared_alpha, int& p, int& dim_max, Filtration_value& min_persistence);
+
+int main(int argc, char* argv[]) {
+  std::string file_name;
+  std::string filediag;
+  Filtration_value max_squared_alpha;
+  int p, nbL, lim_d;
+  Filtration_value min_persistence;
+  SimplexTree simplex_tree;
+
+  program_options(argc, argv, nbL, file_name, filediag, max_squared_alpha, p, lim_d, min_persistence);
+
+  // Extract the points from the file file_name
+  Point_vector witnesses, landmarks;
+  Gudhi::Points_off_reader<Point_d> off_reader(file_name);
+  if (!off_reader.is_valid()) {
+    std::cerr << "Witness complex - Unable to read file " << file_name << "\n";
+    exit(-1);  // ----- >>
+  }
+  witnesses = Point_vector(off_reader.get_point_cloud());
+  std::cout << "Successfully read " << witnesses.size() << " points.\n";
+  std::cout << "Ambient dimension is " << witnesses[0].dimension() << ".\n";
+
+  // Choose landmarks (decomment one of the following two lines)
+  // Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks));
+  Gudhi::subsampling::choose_n_farthest_points(K(), witnesses, nbL, Gudhi::subsampling::random_starting_point,
+                                               std::back_inserter(landmarks));
+
+  // Compute witness complex
+  Strong_witness_complex strong_witness_complex(landmarks, witnesses);
+
+  strong_witness_complex.create_complex(simplex_tree, max_squared_alpha, lim_d);
+
+  std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n";
+  std::cout << "   and has dimension " << simplex_tree.dimension() << " \n";
+
+  // Sort the simplices in the order of the filtration
+  simplex_tree.initialize_filtration();
+
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(simplex_tree);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(p);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  // Output the diagram in filediag
+  if (filediag.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::ofstream out(filediag);
+    pcoh.output_diagram(out);
+    out.close();
+  }
+
+  return 0;
+}
+
+void program_options(int argc, char* argv[], int& nbL, std::string& file_name, std::string& filediag,
+                     Filtration_value& max_squared_alpha, int& p, int& dim_max, Filtration_value& min_persistence) {
+  namespace po = boost::program_options;
+
+  po::options_description hidden("Hidden options");
+  hidden.add_options()("input-file", po::value<std::string>(&file_name),
+                       "Name of file containing a point set in off format.");
+
+  po::options_description visible("Allowed options", 100);
+  Filtration_value default_alpha = std::numeric_limits<Filtration_value>::infinity();
+  visible.add_options()("help,h", "produce help message")("landmarks,l", po::value<int>(&nbL),
+                                                          "Number of landmarks to choose from the point cloud.")(
+      "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()),
+      "Name of file in which the persistence diagram is written. Default print in std::cout")(
+      "max-sq-alpha,a", po::value<Filtration_value>(&max_squared_alpha)->default_value(default_alpha),
+      "Maximal squared relaxation parameter.")(
+      "field-charac,p", po::value<int>(&p)->default_value(11),
+      "Characteristic p of the coefficient field Z/pZ for computing homology.")(
+      "min-persistence,m", po::value<Filtration_value>(&min_persistence)->default_value(0),
+      "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length "
+      "intervals")("cpx-dimension,d", po::value<int>(&dim_max)->default_value(std::numeric_limits<int>::max()),
+                   "Maximal dimension of the strong witness complex we want to compute.");
+
+  po::positional_options_description pos;
+  pos.add("input-file", 1);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+  po::variables_map vm;
+
+  po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  if (vm.count("help") || !vm.count("input-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::cout << "of a Strong witness complex defined on a set of input points.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
+  }
+}
diff --git a/utilities/Witness_complex/weak_witness_persistence.cpp b/utilities/Witness_complex/weak_witness_persistence.cpp
new file mode 100644
index 00000000..23fa93aa
--- /dev/null
+++ b/utilities/Witness_complex/weak_witness_persistence.cpp
@@ -0,0 +1,156 @@
+/*    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):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2016  INRIA (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 <gudhi/Simplex_tree.h>
+#include <gudhi/Euclidean_witness_complex.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Points_off_io.h>
+#include <gudhi/pick_n_random_points.h>
+#include <gudhi/choose_n_farthest_points.h>
+
+#include <boost/program_options.hpp>
+
+#include <CGAL/Epick_d.h>
+
+#include <string>
+#include <vector>
+#include <limits>  // infinity
+
+using K = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
+using Point_d = K::Point_d;
+
+using Point_vector = std::vector<Point_d>;
+using Witness_complex = Gudhi::witness_complex::Euclidean_witness_complex<K>;
+using SimplexTree = Gudhi::Simplex_tree<>;
+
+using Filtration_value = SimplexTree::Filtration_value;
+
+using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
+using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<SimplexTree, Field_Zp>;
+
+void program_options(int argc, char* argv[], int& nbL, std::string& file_name, std::string& filediag,
+                     Filtration_value& max_squared_alpha, int& p, int& dim_max, Filtration_value& min_persistence);
+
+int main(int argc, char* argv[]) {
+  std::string file_name;
+  std::string filediag;
+  Filtration_value max_squared_alpha;
+  int p, nbL, lim_d;
+  Filtration_value min_persistence;
+  SimplexTree simplex_tree;
+
+  program_options(argc, argv, nbL, file_name, filediag, max_squared_alpha, p, lim_d, min_persistence);
+
+  // Extract the points from the file file_name
+  Point_vector witnesses, landmarks;
+  Gudhi::Points_off_reader<Point_d> off_reader(file_name);
+  if (!off_reader.is_valid()) {
+    std::cerr << "Witness complex - Unable to read file " << file_name << "\n";
+    exit(-1);  // ----- >>
+  }
+  witnesses = Point_vector(off_reader.get_point_cloud());
+  std::cout << "Successfully read " << witnesses.size() << " points.\n";
+  std::cout << "Ambient dimension is " << witnesses[0].dimension() << ".\n";
+
+  // Choose landmarks (decomment one of the following two lines)
+  // Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks));
+  Gudhi::subsampling::choose_n_farthest_points(K(), witnesses, nbL, Gudhi::subsampling::random_starting_point,
+                                               std::back_inserter(landmarks));
+
+  // Compute witness complex
+  Witness_complex witness_complex(landmarks, witnesses);
+
+  witness_complex.create_complex(simplex_tree, max_squared_alpha, lim_d);
+
+  std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n";
+  std::cout << "   and has dimension " << simplex_tree.dimension() << " \n";
+
+  // Sort the simplices in the order of the filtration
+  simplex_tree.initialize_filtration();
+
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(simplex_tree);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(p);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  // Output the diagram in filediag
+  if (filediag.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::ofstream out(filediag);
+    pcoh.output_diagram(out);
+    out.close();
+  }
+
+  return 0;
+}
+
+void program_options(int argc, char* argv[], int& nbL, std::string& file_name, std::string& filediag,
+                     Filtration_value& max_squared_alpha, int& p, int& dim_max, Filtration_value& min_persistence) {
+  namespace po = boost::program_options;
+
+  po::options_description hidden("Hidden options");
+  hidden.add_options()("input-file", po::value<std::string>(&file_name),
+                       "Name of file containing a point set in off format.");
+
+  Filtration_value default_alpha = std::numeric_limits<Filtration_value>::infinity();
+  po::options_description visible("Allowed options", 100);
+  visible.add_options()("help,h", "produce help message")("landmarks,l", po::value<int>(&nbL),
+                                                          "Number of landmarks to choose from the point cloud.")(
+      "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()),
+      "Name of file in which the persistence diagram is written. Default print in std::cout")(
+      "max-sq-alpha,a", po::value<Filtration_value>(&max_squared_alpha)->default_value(default_alpha),
+      "Maximal squared relaxation parameter.")(
+      "field-charac,p", po::value<int>(&p)->default_value(11),
+      "Characteristic p of the coefficient field Z/pZ for computing homology.")(
+      "min-persistence,m", po::value<Filtration_value>(&min_persistence)->default_value(0),
+      "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length "
+      "intervals")("cpx-dimension,d", po::value<int>(&dim_max)->default_value(std::numeric_limits<int>::max()),
+                   "Maximal dimension of the weak witness complex we want to compute.");
+
+  po::positional_options_description pos;
+  pos.add("input-file", 1);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+  po::variables_map vm;
+
+  po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  if (vm.count("help") || !vm.count("input-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::cout << "of a Weak witness complex defined on a set of input points.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
+  }
+}
diff --git a/utilities/Witness_complex/witnesscomplex.md b/utilities/Witness_complex/witnesscomplex.md
new file mode 100644
index 00000000..2341759b
--- /dev/null
+++ b/utilities/Witness_complex/witnesscomplex.md
@@ -0,0 +1,66 @@
+

+

+# Witness complex #

+

+

+For more details about the witness complex, please read the [user manual of the package](/doc/latest/group__witness__complex.html).

+

+## weak_witness_persistence ##

+This program computes the persistent homology with coefficient field *Z/pZ* of a Weak witness complex defined on a set of input points.

+The output diagram contains one bar per line, written with the convention:

+

+`p dim birth death`

+

+where `dim` is the dimension of the homological feature, `birth` and `death` are respectively the birth and death of the feature,

+and `p` is the characteristic of the field *Z/pZ* used for homology coefficients.

+

+**Usage**

+

+`weak_witness_persistence [options] <OFF input file>`

+

+**Allowed options**

+

+* `-h [ --help ]` Produce help message

+* `-l [ --landmarks ]` Number of landmarks to choose from the point cloud.

+* `-o [ --output-file ]` Name of file in which the persistence diagram is written. By default, print in std::cout.

+* `-a [ --max-sq-alpha ]` (default = inf) Maximal squared relaxation parameter.

+* `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology.

+* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals.

+* `-d [ --cpx-dimension ]` (default = 2147483647) Maximal dimension of the weak witness complex we want to compute.

+

+**Example**

+

+`weak_witness_persistence data/points/tore3D_1307.off -l 20 -a 0.5 -m 0.006`

+

+N.B.: output is random as the 20 landmarks are chosen randomly.

+

+

+## strong_witness_persistence ##

+

+This program computes the persistent homology with coefficient field *Z/pZ* of a Strong witness complex defined on a set of input points.

+The output diagram contains one bar per line, written with the convention:

+

+`p dim birth death`

+

+where `dim` is the dimension of the homological feature, `birth` and `death` are respectively the birth and death of the feature,

+and `p` is the characteristic of the field *Z/pZ* used for homology coefficients.

+

+**Usage**

+

+`strong_witness_persistence [options] <OFF input file>`

+

+**Allowed options**

+

+* `-h [ --help ]` Produce help message

+* `-l [ --landmarks ]` Number of landmarks to choose from the point cloud.

+* `-o [ --output-file ]` Name of file in which the persistence diagram is written. By default, print in std::cout.

+* `-a [ --max-sq-alpha ]` (default = inf) Maximal squared relaxation parameter.

+* `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology.

+* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals.

+* `-d [ --cpx-dimension ]` (default = 2147483647) Maximal dimension of the weak witness complex we want to compute.

+

+**Example**

+

+`strong_witness_persistence data/points/tore3D_1307.off -l 20 -a 0.5 -m 0.06`

+

+N.B.: output is random as the 20 landmarks are chosen randomly.

diff --git a/utilities/common/README b/utilities/common/README
deleted file mode 100644
index dc841521..00000000
--- a/utilities/common/README
+++ /dev/null
@@ -1,19 +0,0 @@
-======================= off_file_from_shape_generator ==================================
-
-Example of use :
-
-*** on|in sphere|cube|curve|torus|klein generator
-
-./off_file_from_shape_generator on sphere onSphere.off 1000 3 15.2
-
- => generates a onSphere.off file with 1000 points randomized on a sphere of dimension 3 and radius 15.2
-
-./off_file_from_shape_generator in sphere inSphere.off 100 2
-
- => generates a inSphere.off file with 100 points randomized in a sphere of dimension 2 (circle) and radius 1.0 (default)
-
-./off_file_from_shape_generator in cube inCube.off 10000 3 5.8
-
- => generates a inCube.off file with 10000 points randomized in a cube of dimension 3 and side 5.8
-
-!! Warning: hypegenerator on cube is not available !!
diff --git a/utilities/common/pointsetgenerator.md b/utilities/common/pointsetgenerator.md
new file mode 100644
index 00000000..284715d4
--- /dev/null
+++ b/utilities/common/pointsetgenerator.md
@@ -0,0 +1,33 @@
+
+
+# common #
+
+## off_file_from_shape_generator ##
+
+Generates a pointset and save it in an OFF file. Command-line is:
+
+```
+off_file_from_shape_generator on|in sphere|cube|curve|torus|klein <filename> <num_points> <dimension> <parameter1> <parameter2>...
+```
+
+Warning: "on cube" generator is not available!
+
+**Examples**
+
+```
+off_file_from_shape_generator on sphere onSphere.off 1000 3 15.2
+```
+
+* Generates an onSphere.off file with 1000 points randomized on a sphere of dimension 3 and radius 15.2.
+
+```
+off_file_from_shape_generator in sphere inSphere.off 100 2
+```
+
+* Generates an inSphere.off file with 100 points randomized in a sphere of dimension 2 (circle) and radius 1.0 (default).
+
+```
+off_file_from_shape_generator in cube inCube.off 10000 3 5.8
+```
+
+* Generates a inCube.off file with 10000 points randomized in a cube of dimension 3 and side 5.8.