Merge branch 'edge_collase_vincent' into edge_collapse_integration_vincent

5 files changed, 1495 insertions, 0 deletions

diff --git a/src/Collapse/include/gudhi/FlagComplexSpMatrix.h b/src/Collapse/include/gudhi/FlagComplexSpMatrix.h
new file mode 100644
index 00000000..ac02a46f
--- /dev/null
+++ b/src/Collapse/include/gudhi/FlagComplexSpMatrix.h
@@ -0,0 +1,852 @@
+/*    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):       Siddharth Pritam
+ *
+ *    Copyright (C) 2018 INRIA Sophia Antipolis (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/>.
+
+*/
+#pragma once
+
+#include <gudhi/Rips_edge_list.h>
+#include <boost/functional/hash.hpp>
+// #include <boost/graph/adjacency_list.hpp>
+
+#include <iostream>
+#include <utility>
+#include <vector>
+#include <queue>
+#include <unordered_map>
+#include <tuple>
+#include <list>
+#include <algorithm>
+#include <chrono>
+
+#include <ctime>
+#include <fstream>
+
+#include <Eigen/Sparse>
+
+typedef std::size_t Vertex;
+using Edge = std::pair<Vertex, Vertex>;  // This is an ordered pair, An edge is stored with convention of the first
+                                         // element being the smaller i.e {2,3} not {3,2}. However this is at the level
+                                         // of row indices on actual vertex lables
+using EdgeFilt = std::pair<Edge, double>;
+using edge_list = std::vector<Edge>;
+
+using MapVertexToIndex = std::unordered_map<Vertex, std::size_t>;
+using Map = std::unordered_map<Vertex, Vertex>;
+
+using sparseRowMatrix = Eigen::SparseMatrix<double, Eigen::RowMajor>;
+using rowInnerIterator = sparseRowMatrix::InnerIterator;
+
+using doubleVector = std::vector<double>;
+using vertexVector = std::vector<Vertex>;
+using boolVector = std::vector<bool>;
+
+using doubleQueue = std::queue<double>;
+
+using EdgeFiltQueue = std::queue<EdgeFilt>;
+using EdgeFiltVector = std::vector<EdgeFilt>;
+
+typedef std::vector<std::tuple<double, Vertex, Vertex>> Filtered_sorted_edge_list;
+typedef std::unordered_map<Edge, bool, boost::hash<Edge>> u_edge_map;
+typedef std::unordered_map<Edge, std::size_t, boost::hash<Edge>> u_edge_to_idx_map;
+
+//!  Class SparseMsMatrix
+/*!
+  The class for storing the Vertices v/s MaxSimplices Sparse Matrix and performing collapses operations using the N^2()
+  Algorithm.
+*/
+class FlagComplexSpMatrix {
+ private:
+  std::unordered_map<int, Vertex> rowToVertex;
+
+  // Vertices strored as an unordered_set
+  std::unordered_set<Vertex> vertices;
+  //! Stores the 1-simplices(edges) of the original Simplicial Complex.
+  edge_list oneSimplices;
+
+  // Unordered set of  removed edges. (to enforce removal from the matrix)
+  std::unordered_set<Edge, boost::hash<Edge>> u_set_removed_redges;
+
+  // Unordered set of  dominated edges. (to inforce removal from the matrix)
+  std::unordered_set<Edge, boost::hash<Edge>> u_set_dominated_redges;
+
+  // Map from egde to its index
+  u_edge_to_idx_map edge_to_index_map;
+  // Boolean vector to indicate if the index is critical or not.
+  boolVector critical_edge_indicator;  // critical indicator
+
+  // Boolean vector to indicate if the index is critical or not.
+  boolVector dominated_edge_indicator;  // domination indicator
+
+  //! Stores the Map between vertices<B>rowToVertex  and row indices <B>rowToVertex -> row-index</B>.
+  /*!
+    \code
+    MapVertexToIndex = std::unordered_map<Vertex,int>
+    \endcode
+    So, if the original simplex tree had vertices 0,1,4,5 <br>
+    <B>rowToVertex</B> would store : <br>
+    \verbatim
+    Values =  | 0 | 1 | 4 | 5 |
+    Indices =   0   1   2   3
+    \endverbatim
+    And <B>vertexToRow</B> would be a map like the following : <br>
+    \verbatim
+    0 -> 0
+    1 -> 1
+    4 -> 2
+    5 -> 3
+    \endverbatim
+  */
+  MapVertexToIndex vertexToRow;
+
+  //! Stores the number of vertices in the original Simplicial Complex.
+  /*!
+    This stores the count of vertices (which is also the number of rows in the Matrix).
+  */
+  std::size_t rows;
+
+  std::size_t numOneSimplices;
+
+  std::size_t numDomEdge;
+
+  //! Stores the Sparse matrix of double values representing the Original Simplicial Complex.
+  /*!
+    \code
+    sparseRowMatrix   = Eigen::SparseMatrix<double, Eigen::RowMajor> ;
+    \endcode
+    ;
+      */
+
+  sparseRowMatrix* sparse_colpsd_adj_Matrix;  // Stores the collapsed sparse matrix representaion.
+  sparseRowMatrix sparseRowAdjMatrix;  // This is row-major version of the same sparse-matrix, to facilitate easy access
+                                       // to elements when traversing the matrix row-wise.
+
+  //! Stores <I>true</I> for dominated rows and  <I>false</I> for undominated rows.
+  /*!
+    Initialised to a vector of length equal to the value of the variable <B>rows</B> with all <I>false</I> values.
+    Subsequent removal of dominated vertices is reflected by concerned entries changing to <I>true</I> in this vector.
+  */
+  boolVector vertDomnIndicator;  //(domination indicator)
+
+  boolVector contractionIndicator;  //(contraction indicator)
+
+  //! Stores the indices of the rows to-be checked for domination in the current iteration.
+  /*!
+    Initialised with all rows for the first iteration.
+    Subsequently once a dominated row is found, its non-dominated neighbhour indices are inserted.
+  */
+  // doubleQueue rowIterator;
+
+  doubleQueue rowIterator;
+
+  // Queue of filtered edges, for edge-collapse, the indices of the edges are the row-indices.
+  EdgeFiltQueue filteredEgdeIter;
+
+  // Vector of filtered edges, for edge-collapse, the indices of the edges are the row-indices.
+  EdgeFiltVector fEgdeVector;
+
+  // List of non-dominated edges,  the indices of the edges are the vertex lables!!.
+  Filtered_sorted_edge_list criticalCoreEdges;
+  // Stores the indices from the sorted filtered edge vector.
+  // std::set<std::size_t> recurCriticalCoreIndcs;
+
+  //! Stores <I>true</I> if the current row is inserted in the queue <B>rowIterator<B> otherwise its value is
+  //! <I>false<I>.
+  /*!
+    Initialised to a boolean vector of length equal to the value of the variable <B>rows</B> with all <I>true</I>
+    values. Subsequent removal/addition of a row from <B>rowIterator<B> is reflected by concerned entries changing to
+    <I>false</I>/<I>true</I> in this vector.
+  */
+  boolVector rowInsertIndicator;  //(current iteration row insertion indicator)
+
+  //! Map that stores the Reduction / Collapse of vertices.
+  /*!
+    \code
+    Map = std::unordered_map<Vertex,Vertex>
+    \endcode
+    This is empty to begin with. As and when collapses are done (let's say from dominated vertex <I>v</I> to dominating
+    vertex <I>v'</I>) : <br> <B>ReductionMap</B>[<I>v</I>] = <I>v'</I> is entered into the map. <br> <I>This does not
+    store uncollapsed vertices. What it means is that say vertex <I>x</I> was never collapsed onto any other vertex.
+    Then, this map <B>WILL NOT</B> have any entry like <I>x</I> -> <I>x</I>. Basically, it will have no entry
+    corresponding to vertex <I>x</I> at all. </I>
+  */
+  Map ReductionMap;
+
+  bool vertexCollapsed;
+  bool edgeCollapsed;
+  // Variable to indicate if filtered-edge-collapse has to be performed.
+  int expansion_limit;
+
+  void init() {
+    rowToVertex.clear();
+    vertexToRow.clear();
+    oneSimplices.clear();
+    ReductionMap.clear();
+
+    vertDomnIndicator.clear();
+    rowInsertIndicator.clear();
+    rowIterator.push(0);
+    rowIterator.pop();
+
+    filteredEgdeIter.push({{0, 0}, 0});
+    filteredEgdeIter.pop();
+    fEgdeVector.clear();
+
+    rows = 0;
+    numDomEdge = 0;
+
+    numOneSimplices = 0;
+    expansion_limit = 2;
+
+    vertexCollapsed = false;
+    edgeCollapsed = false;
+  }
+
+  //!	Function for computing the sparse-matrix corresponding to the core of the complex. It also prepares the working
+  //!list filteredEgdeIter for edge collapses
+  void after_vertex_collapse() {
+    sparse_colpsd_adj_Matrix = new sparseRowMatrix(rows, rows);  // Just for debugging purpose.
+    oneSimplices.clear();
+    if (not filteredEgdeIter.empty()) {
+      std::cout << "Working list for edge collapses are not empty before the edge-collapse." << std::endl;
+    }
+    for (std::size_t rw = 0; rw < rows; ++rw) {
+      if (not vertDomnIndicator[rw])  // If the current column is not dominated
+      {
+        auto nbhrs_to_insert = closed_neighbours_row_index(rw);  // returns row indices of the non-dominated vertices.
+        for (auto& v : nbhrs_to_insert) {
+          sparse_colpsd_adj_Matrix->insert(rw, v) = 1;  // This creates the full matrix
+          if (rw < v) {
+            oneSimplices.push_back({rowToVertex[rw], rowToVertex[v]});
+            filteredEgdeIter.push({{rw, v}, 1});
+            // if(rw == v)
+            // std::cout << "Pushed the edge {" << rw << ", " << v <<  "} " << std::endl;
+          }
+        }
+      }
+    }
+    // std::cout << "Total number of non-zero elements before domination check are: " <<
+    // sparse_colpsd_adj_Matrix->nonZeros() << std::endl; std::cout << "Total number of edges for domination check are:
+    // " << filteredEgdeIter.size() << std::endl;
+    // std::cout << *sparse_colpsd_adj_Matrix << std::endl;
+    return;
+  }
+
+  //! Function to fully compact a particular vertex of the ReductionMap.
+  /*!
+    It takes as argument the iterator corresponding to a particular vertex pair (key-value) stored in the ReductionMap.
+    <br> It then checks if the second element of this particular vertex pair is present as a first element of some other
+    key-value pair in the map. If no, then the first element of the vertex pair in consideration is fully compact. If
+    yes, then recursively call fully_compact_this_vertex() on the second element of the original pair in consideration
+    and assign its resultant image as the image of the first element of the original pair in consideration as well.
+  */
+  void fully_compact_this_vertex(Map::iterator iter) {
+    Map::iterator found = ReductionMap.find(iter->second);
+    if (found == ReductionMap.end()) return;
+
+    fully_compact_this_vertex(found);
+    iter->second = ReductionMap[iter->second];
+  }
+
+  //! Function to fully compact the Reduction Map.
+  /*!
+    While doing strong collapses, we store only the immediate collapse of a vertex. Which means that in one round,
+    vertex <I>x</I> may collapse to vertex <I>y</I>. And in some later round it may be possible that vertex <I>y</I>
+    collapses to <I>z</I>. In which case our map stores : <br> <I>x</I> -> <I>y</I> and also <I>y</I> -> <I>z</I>. But
+    it really should store : <I>x</I> -> <I>z</I> and <I>y</I> -> <I>z</I>. This function achieves the same. <br> It
+    basically calls fully_compact_this_vertex() for each entry in the map.
+  */
+  void fully_compact() {
+    Map::iterator it = ReductionMap.begin();
+    while (it != ReductionMap.end()) {
+      fully_compact_this_vertex(it);
+      it++;
+    }
+  }
+
+  void sparse_strong_vertex_collapse() {
+    complete_vertex_domination_check(
+        rowIterator, rowInsertIndicator,
+        vertDomnIndicator);  // Complete check for rows in rowIterator, rowInsertIndicator is a list of boolean
+                             // indicator if a vertex is already inserted in the working row_queue (rowIterator)
+    if (not rowIterator.empty())
+      sparse_strong_vertex_collapse();
+    else
+      return;
+  }
+
+  void complete_vertex_domination_check(doubleQueue& iterator, boolVector& insertIndicator, boolVector& domnIndicator) {
+    double k;
+    doubleVector nonZeroInnerIdcs;
+    while (not iterator.empty())  // "iterator" contains list(FIFO) of rows to be considered for domination check
+    {
+      k = iterator.front();
+      iterator.pop();
+      insertIndicator[k] = false;
+      if (not domnIndicator[k])  // Check if is  already dominated
+      {
+        nonZeroInnerIdcs = closed_neighbours_row_index(k);
+        for (doubleVector::iterator it = nonZeroInnerIdcs.begin(); it != nonZeroInnerIdcs.end(); it++) {
+          int checkDom = vertex_domination_check(k, *it);  // "true" for row domination comparison
+          if (checkDom == 1)                               // row k is dominated by *it, k <= *it;
+          {
+            setZero(k, *it);
+            break;
+          } else if (checkDom == -1)  // row *it is dominated by k, *it <= k;
+            setZero(*it, k);
+        }
+      }
+    }
+  }
+
+  bool check_edge_domination(Edge e)  // Edge e is the actual edge (u,v). Not the row ids in the matrixs
+  {
+    auto u = std::get<0>(e);
+    auto v = std::get<1>(e);
+
+    auto rw_u = vertexToRow[u];
+    auto rw_v = vertexToRow[v];
+    auto rw_e = std::make_pair(rw_u, rw_v);
+    // std::cout << "The edge {" << u << ", " << v <<  "} is going for domination check." << std::endl;
+    auto commonNeighbours = closed_common_neighbours_row_index(rw_e);
+    // std::cout << "And its common neighbours are." << std::endl;
+    // for (doubleVector::iterator it = commonNeighbours.begin(); it!=commonNeighbours.end(); it++) {
+    // std::cout << rowToVertex[*it] << ", " ;
+    // }
+    // std::cout<< std::endl;
+    if (commonNeighbours.size() > 2) {
+      if (commonNeighbours.size() == 3)
+        return true;
+      else
+        for (doubleVector::iterator it = commonNeighbours.begin(); it != commonNeighbours.end(); it++) {
+          auto rw_c = *it;  // Typecasting
+          if (rw_c != rw_u and rw_c != rw_v) {
+            auto neighbours_c = closed_neighbours_row_index(rw_c);
+            if (std::includes(neighbours_c.begin(), neighbours_c.end(), commonNeighbours.begin(),
+                              commonNeighbours.end()))  // If neighbours_c contains the common neighbours.
+              return true;
+          }
+        }
+    }
+    return false;
+  }
+
+  bool check_domination_indicator(Edge e)  // The edge should be sorted by the indices and indices are original
+  {
+    return dominated_edge_indicator[edge_to_index_map[e]];
+  }
+
+  std::set<std::size_t> three_clique_indices(std::size_t crit) {
+    std::set<std::size_t> edge_indices;
+
+    EdgeFilt fe = fEgdeVector.at(crit);
+    Edge e = std::get<0>(fe);
+    Vertex u = std::get<0>(e);
+    Vertex v = std::get<1>(e);
+
+    // std::cout << "The  current critical edge to re-check criticality with filt value is : {" << u << "," << v << "};
+    // "<< std::get<1>(fe) << std::endl;
+    auto rw_u = vertexToRow[u];
+    auto rw_v = vertexToRow[v];
+    auto rw_critical_edge = std::make_pair(rw_u, rw_v);
+
+    doubleVector commonNeighbours = closed_common_neighbours_row_index(rw_critical_edge);
+
+    if (commonNeighbours.size() > 2) {
+      for (doubleVector::iterator it = commonNeighbours.begin(); it != commonNeighbours.end(); it++) {
+        auto rw_c = *it;
+        if (rw_c != rw_u and rw_c != rw_v) {
+          auto e_with_new_nbhr_v = std::minmax(u, rowToVertex[rw_c]);
+          auto e_with_new_nbhr_u = std::minmax(v, rowToVertex[rw_c]);
+          edge_indices.emplace(edge_to_index_map[e_with_new_nbhr_v]);
+          edge_indices.emplace(edge_to_index_map[e_with_new_nbhr_u]);
+        }
+      }
+    }
+    return edge_indices;
+  }
+
+  void set_edge_critical(std::size_t indx, double filt) {
+    // std::cout << "The curent index  with filtration value " << indx << ", " << filt << " is primary critical" <<
+    // std::endl;
+    std::set<std::size_t> effectedIndcs = three_clique_indices(indx);
+    if (effectedIndcs.size() > 0) {
+      for (auto idx = indx - 1; idx > 0; idx--) {
+        EdgeFilt fec = fEgdeVector.at(idx);
+        Edge e = std::get<0>(fec);
+        Vertex u = std::get<0>(e);
+        Vertex v = std::get<1>(e);
+        if (not critical_edge_indicator.at(
+                idx)) {  // If idx is not critical so it should be proceses, otherwise it stays in the graph // prev
+                         // code : recurCriticalCoreIndcs.find(idx) == recurCriticalCoreIndcs.end()
+          if (effectedIndcs.find(idx) != effectedIndcs.end()) {  // If idx is affected
+            if (not check_edge_domination(e)) {
+              // std::cout << "The curent index is became critical " << idx  << std::endl;
+              critical_edge_indicator.at(idx) = true;
+              criticalCoreEdges.push_back({filt, u, v});
+              std::set<std::size_t> inner_effected_indcs = three_clique_indices(idx);
+              for (auto inr_idx = inner_effected_indcs.rbegin(); inr_idx != inner_effected_indcs.rend(); inr_idx++) {
+                if (*inr_idx < idx) effectedIndcs.emplace(*inr_idx);
+              }
+              inner_effected_indcs.clear();
+              // std::cout << "The following edge is critical with filt value: {" << std::get<0>(e) << "," <<
+              // std::get<1>(e) << "}; " << filt << std::endl;
+            } else
+              u_set_dominated_redges.emplace(std::minmax(vertexToRow[u], vertexToRow[v]));
+          } else  // Idx is not affected hence dominated.
+            u_set_dominated_redges.emplace(std::minmax(vertexToRow[u], vertexToRow[v]));
+        }
+      }
+    }
+    effectedIndcs.clear();
+    u_set_dominated_redges.clear();
+  }
+
+  void critical_core_edges() {
+    std::size_t totEdges = fEgdeVector.size();
+
+    std::size_t endIdx = 0;
+
+    u_set_removed_redges.clear();
+    u_set_dominated_redges.clear();
+    critical_edge_indicator.clear();
+
+    while (endIdx < totEdges) {
+      EdgeFilt fec = fEgdeVector.at(endIdx);
+
+      insert_new_edges(std::get<0>(std::get<0>(fec)), std::get<1>(std::get<0>(fec)),
+                       std::get<1>(fec));  // Inserts the edge in the sparse matrix to update the graph (G_i)
+
+      Edge e = std::get<0>(fec);
+      Vertex u = std::get<0>(e);
+      Vertex v = std::get<1>(e);
+      edge_to_index_map.emplace(std::minmax(u, v), endIdx);
+      critical_edge_indicator.push_back(false);
+      dominated_edge_indicator.push_back(false);
+
+      if (not check_edge_domination(e)) {
+        critical_edge_indicator.at(endIdx) = true;
+        dominated_edge_indicator.at(endIdx) = false;
+        criticalCoreEdges.push_back({std::get<1>(fec), u, v});
+        if (endIdx > 1) set_edge_critical(endIdx, std::get<1>(fec));
+
+      } else
+        dominated_edge_indicator.at(endIdx) = true;
+      endIdx++;
+    }
+
+    std::cout << "The total number of critical edges is: " << criticalCoreEdges.size() << std::endl;
+    std::cout << "The total number of non-critical edges is: " << totEdges - criticalCoreEdges.size() << std::endl;
+  }
+
+  int vertex_domination_check(double i, double j)  // True for row comparison, false for column comparison
+  {
+    if (i != j) {
+      doubleVector Listi = closed_neighbours_row_index(i);
+      doubleVector Listj = closed_neighbours_row_index(j);
+      if (Listj.size() <= Listi.size()) {
+        if (std::includes(Listi.begin(), Listi.end(), Listj.begin(), Listj.end()))  // Listj is a subset of Listi
+          return -1;
+      }
+
+      else if (std::includes(Listj.begin(), Listj.end(), Listi.begin(), Listi.end()))  // Listi is a subset of Listj
+        return 1;
+    }
+    return 0;
+  }
+
+  doubleVector closed_neighbours_row_index(double indx)  // Returns list of non-zero columns of the particular indx.
+  {
+    doubleVector nonZeroIndices;
+    Vertex u = indx;
+    Vertex v;
+    // std::cout << "The neighbours of the vertex: " << rowToVertex[u] << " are. " << std::endl;
+    if (not vertDomnIndicator[indx]) {
+      for (rowInnerIterator it(sparseRowAdjMatrix, indx); it; ++it) {  // Iterate over the non-zero columns
+        v = it.index();
+        if (not vertDomnIndicator[v] and u_set_removed_redges.find(std::minmax(u, v)) == u_set_removed_redges.end() and
+            u_set_dominated_redges.find(std::minmax(u, v)) ==
+                u_set_dominated_redges
+                    .end()) {  // If the vertex v is not dominated and the edge {u,v} is still in the matrix
+          nonZeroIndices.push_back(it.index());  // inner index, here it is equal to it.columns()
+                                                 // std::cout << rowToVertex[it.index()] << ", " ;
+        }
+      }
+      // std::cout << std::endl;
+    }
+    return nonZeroIndices;
+  }
+
+  doubleVector closed_common_neighbours_row_index(Edge e)  // Returns the list of closed neighbours of the edge :{u,v}.
+  {
+    doubleVector common;
+    doubleVector nonZeroIndices_u;
+    doubleVector nonZeroIndices_v;
+    double u = std::get<0>(e);
+    double v = std::get<1>(e);
+
+    nonZeroIndices_u = closed_neighbours_row_index(u);
+    nonZeroIndices_v = closed_neighbours_row_index(v);
+    std::set_intersection(nonZeroIndices_u.begin(), nonZeroIndices_u.end(), nonZeroIndices_v.begin(),
+                          nonZeroIndices_v.end(), std::inserter(common, common.begin()));
+
+    return common;
+  }
+
+  void setZero(double dominated, double dominating) {
+    for (auto& v : closed_neighbours_row_index(dominated))
+      if (not rowInsertIndicator[v])  // Checking if the row is already inserted
+      {
+        rowIterator.push(v);
+        rowInsertIndicator[v] = true;
+      }
+    vertDomnIndicator[dominated] = true;
+    ReductionMap[rowToVertex[dominated]] = rowToVertex[dominating];
+
+    vertexToRow.erase(rowToVertex[dominated]);
+    vertices.erase(rowToVertex[dominated]);
+    rowToVertex.erase(dominated);
+  }
+
+  vertexVector closed_neighbours_vertex_index(double rowIndx)  // Returns list of non-zero "vertices" of the particular
+                                                               // colIndx. the difference is in the return type
+  {
+    vertexVector colmns;
+    for (auto& v : closed_neighbours_row_index(rowIndx))  // Iterate over the non-zero columns
+      colmns.push_back(rowToVertex[v]);
+    std::sort(colmns.begin(), colmns.end());
+    return colmns;
+  }
+
+  vertexVector vertex_closed_active_neighbours(
+      double rowIndx)  // Returns list of all non-zero "vertices" of the particular colIndx which are currently active.
+                       // the difference is in the return type.
+  {
+    vertexVector colmns;
+    for (auto& v : closed_neighbours_row_index(rowIndx))  // Iterate over the non-zero columns
+      if (not contractionIndicator[v])                    // Check if the row corresponds to a contracted vertex
+        colmns.push_back(rowToVertex[v]);
+    std::sort(colmns.begin(), colmns.end());
+    return colmns;
+  }
+
+  vertexVector closed_all_neighbours_row_index(
+      double rowIndx)  // Returns list of all non-zero "vertices" of the particular colIndx whether dominated or not.
+                       // the difference is in the return type.
+  {
+    vertexVector colmns;
+    for (rowInnerIterator itCol(sparseRowAdjMatrix, rowIndx); itCol; ++itCol)  // Iterate over the non-zero columns
+      colmns.push_back(rowToVertex[itCol.index()]);  // inner index, here it is equal to it.row()
+    std::sort(colmns.begin(), colmns.end());
+    return colmns;
+  }
+
+  void swap_rows(const Vertex& v,
+                 const Vertex& w) {  // swap the rows of v and w. Both should be members of the skeleton
+    if (membership(v) && membership(w)) {
+      auto rw_v = vertexToRow[v];
+      auto rw_w = vertexToRow[w];
+      vertexToRow[v] = rw_w;
+      vertexToRow[w] = rw_v;
+      rowToVertex[rw_v] = w;
+      rowToVertex[rw_w] = v;
+    }
+  }
+
+ public:
+  //! Default Constructor
+  /*!
+    Only initialises all Data Members of the class to empty/Null values as appropriate.
+    One <I>WILL</I> have to create the matrix using the Constructor that has an object of the Simplex_tree class as
+    argument.
+  */
+
+  FlagComplexSpMatrix() { init(); }
+
+  FlagComplexSpMatrix(std::size_t expRows) {
+    init();
+    sparseRowAdjMatrix = sparseRowMatrix(
+        expansion_limit * expRows,
+        expansion_limit * expRows);  // Initializing sparseRowAdjMatrix, This is a row-major sparse matrix.
+  }
+
+  //! Main Constructor
+  /*!
+    Argument is an instance of Filtered_sorted_edge_list. <br>
+    This is THE function that initialises all data members to appropriate values. <br>
+    <B>rowToVertex</B>, <B>vertexToRow</B>, <B>rows</B>, <B>cols</B>, <B>sparseRowAdjMatrix</B> are initialised here.
+    <B>vertDomnIndicator</B>, <B>rowInsertIndicator</B> ,<B>rowIterator<B> are initialised by init() function which is
+    called at the begining of this. <br>
+  */
+  FlagComplexSpMatrix(const size_t& num_vertices, const Filtered_sorted_edge_list& edge_t) {
+    init();
+    sparseRowAdjMatrix = sparseRowMatrix(
+        num_vertices, num_vertices);  // Initializing sparseRowAdjMatrix, This is a row-major sparse matrix.
+
+    for (size_t bgn_idx = 0; bgn_idx < edge_t.size(); bgn_idx++) {
+      // std::vector<size_t>  s = {std::get<1>(edge_t.at(bgn_idx)), std::get<2>(edge_t.at(bgn_idx))};
+      // insert_new_edges(std::get<1>(edge_t.at(bgn_idx)), std::get<2>(edge_t.at(bgn_idx)), 1);
+      fEgdeVector.push_back(
+          {{std::get<1>(edge_t.at(bgn_idx)), std::get<2>(edge_t.at(bgn_idx))}, std::get<0>(edge_t.at(bgn_idx))});
+    }
+
+    // sparseRowAdjMatrix.makeCompressed();
+  }
+
+  //!	Destructor.
+  /*!
+    Frees up memory locations on the heap.
+  */
+  ~FlagComplexSpMatrix() {}
+
+  //!	Function for performing strong collapse.
+  /*!
+    calls sparse_strong_vertex_collapse(), and
+    Then, it compacts the ReductionMap by calling the function fully_compact().
+  */
+  double strong_vertex_collapse() {
+    auto begin_collapse = std::chrono::high_resolution_clock::now();
+    sparse_strong_vertex_collapse();
+    vertexCollapsed = true;
+    auto end_collapse = std::chrono::high_resolution_clock::now();
+
+    auto collapseTime = std::chrono::duration<double, std::milli>(end_collapse - begin_collapse).count();
+    // std::cout << "Time of Collapse : " << collapseTime << " ms\n" << std::endl;
+
+    // Now we complete the Reduction Map
+    fully_compact();
+    // Post processing...
+    after_vertex_collapse();
+    return collapseTime;
+  }
+
+  // Performs edge collapse in a decreasing sequence of the filtration value.
+  Filtered_sorted_edge_list filtered_edge_collapse() {
+    critical_core_edges();
+    vertexCollapsed = false;
+    edgeCollapsed = true;
+    return criticalCoreEdges;
+  }
+
+  // double strong_vertex_edge_collapse() {
+  // 	auto begin_collapse  = std::chrono::high_resolution_clock::now();
+  // 	strong_vertex_collapse();
+  // 	strong_edge_collapse();
+  // 	// strong_vertex_collapse();
+  // 	auto end_collapse = std::chrono::high_resolution_clock::now();
+
+  // 	auto collapseTime = std::chrono::duration<double, std::milli>(end_collapse- begin_collapse).count();
+  // 	return collapseTime;
+  // }
+
+  bool membership(const Vertex& v) {
+    auto rw = vertexToRow.find(v);
+    if (rw != vertexToRow.end())
+      return true;
+    else
+      return false;
+  }
+
+  bool membership(const Edge& e) {
+    auto u = std::get<0>(e);
+    auto v = std::get<1>(e);
+    if (membership(u) && membership(v)) {
+      auto rw_u = vertexToRow[u];
+      auto rw_v = vertexToRow[v];
+      if (rw_u <= rw_v)
+        for (auto x : closed_neighbours_row_index(rw_v)) {  // Taking advantage of sorted lists.
+          if (rw_u == x)
+            return true;
+          else if (rw_u < x)
+            return false;
+        }
+      else
+        for (auto x : closed_neighbours_row_index(rw_u)) {  // Taking advantage of sorted lists.
+          if (rw_v == x)
+            return true;
+          else if (rw_v < x)
+            return false;
+        }
+    }
+    return false;
+  }
+  void insert_vertex(const Vertex& vertex, double filt_val) {
+    auto rw = vertexToRow.find(vertex);
+    if (rw == vertexToRow.end()) {
+      sparseRowAdjMatrix.insert(rows, rows) =
+          filt_val;  // Initializing the diagonal element of the adjency matrix corresponding to rw_b.
+      vertDomnIndicator.push_back(false);
+      rowInsertIndicator.push_back(true);
+      contractionIndicator.push_back(false);
+      rowIterator.push(rows);
+      vertexToRow.insert(std::make_pair(vertex, rows));
+      rowToVertex.insert(std::make_pair(rows, vertex));
+      vertices.emplace(vertex);
+      rows++;
+    }
+  }
+
+  void insert_new_edges(const Vertex& u, const Vertex& v,
+                        double filt_val)  // The edge must not be added before, it should be a new edge.
+  {
+    insert_vertex(u, filt_val);
+    if (u != v) {
+      insert_vertex(v, filt_val);
+      // std::cout << "Insertion of the edge begins " << u <<", " << v << std::endl;
+
+      auto rw_u = vertexToRow.find(u);
+      auto rw_v = vertexToRow.find(v);
+      // std::cout << "Inserting the edge " << u <<", " << v << std::endl;
+      sparseRowAdjMatrix.insert(rw_u->second, rw_v->second) = filt_val;
+      sparseRowAdjMatrix.insert(rw_v->second, rw_u->second) = filt_val;
+      oneSimplices.emplace_back(u, v);
+      numOneSimplices++;
+    }
+    // else
+    // 	std::cout << "Already a member simplex,  skipping..." << std::endl;
+  }
+
+  std::size_t num_vertices() const { return vertices.size(); }
+
+  //!	Function for returning the ReductionMap.
+  /*!
+    This is the (stl's unordered) map that stores all the collapses of vertices. <br>
+    It is simply returned.
+  */
+
+  Map reduction_map() const { return ReductionMap; }
+  std::unordered_set<Vertex> vertex_set() const { return vertices; }
+  sparseRowMatrix collapsed_matrix() const { return *sparse_colpsd_adj_Matrix; }
+
+  sparseRowMatrix uncollapsed_matrix() const { return sparseRowAdjMatrix; }
+
+  edge_list all_edges() const { return oneSimplices; }
+
+  vertexVector active_neighbors(const Vertex& v) {
+    vertexVector nb;
+    auto rw_v = vertexToRow.find(v);
+    if (rw_v != vertexToRow.end()) {
+      nb = vertex_closed_active_neighbours(rw_v->second);
+    }
+    return nb;
+  }
+
+  vertexVector neighbors(const Vertex& v) {
+    vertexVector nb;
+    auto rw_v = vertexToRow.find(v);
+    if (rw_v != vertexToRow.end()) nb = closed_neighbours_vertex_index(rw_v->second);
+
+    return nb;
+  }
+
+  vertexVector active_relative_neighbors(const Vertex& v, const Vertex& w) {
+    std::vector<Vertex> diff;
+    if (membership(v) && membership(w)) {
+      auto nbhrs_v = active_neighbors(v);
+      auto nbhrs_w = active_neighbors(w);
+      std::set_difference(nbhrs_v.begin(), nbhrs_v.end(), nbhrs_w.begin(), nbhrs_w.end(),
+                          std::inserter(diff, diff.begin()));
+    }
+    return diff;
+  }
+
+  void contraction(const Vertex& del, const Vertex& keep) {
+    if (del != keep) {
+      bool del_mem = membership(del);
+      bool keep_mem = membership(keep);
+      if (del_mem && keep_mem) {
+        doubleVector del_indcs, keep_indcs, diff;
+        auto row_del = vertexToRow[del];
+        auto row_keep = vertexToRow[keep];
+        del_indcs = closed_neighbours_row_index(row_del);
+        keep_indcs = closed_neighbours_row_index(row_keep);
+        std::set_difference(del_indcs.begin(), del_indcs.end(), keep_indcs.begin(), keep_indcs.end(),
+                            std::inserter(diff, diff.begin()));
+        for (auto& v : diff) {
+          if (v != row_del) {
+            sparseRowAdjMatrix.insert(row_keep, v) = 1;
+            sparseRowAdjMatrix.insert(v, row_keep) = 1;
+          }
+        }
+        vertexToRow.erase(del);
+        vertices.erase(del);
+        rowToVertex.erase(row_del);
+        // setZero(row_del->second, row_keep->second);
+      } else if (del_mem && not keep_mem) {
+        vertexToRow.insert(std::make_pair(keep, vertexToRow.find(del)->second));
+        rowToVertex[vertexToRow.find(del)->second] = keep;
+        vertices.emplace(keep);
+        vertices.erase(del);
+        vertexToRow.erase(del);
+
+      } else {
+        std::cerr << "The first vertex entered in the method contraction() doesn't exist in the skeleton." << std::endl;
+        exit(-1);
+      }
+    }
+  }
+
+  void relable(const Vertex& v, const Vertex& w) {  // relable v as w.
+    if (membership(v) and v != w) {
+      auto rw_v = vertexToRow[v];
+      rowToVertex[rw_v] = w;
+      vertexToRow.insert(std::make_pair(w, rw_v));
+      vertices.emplace(w);
+      vertexToRow.erase(v);
+      vertices.erase(v);
+      // std::cout<< "Re-named the vertex " << v << " as " << w << std::endl;
+    }
+  }
+
+  // Returns the contracted edge. along with the contracted vertex in the begining of the list at {u,u} or {v,v}
+
+  void active_strong_expansion(const Vertex& v, const Vertex& w, double filt_val) {
+    if (membership(v) && membership(w) && v != w) {
+      // std::cout << "Strong expansion of the vertex " << v << " and " << w << " begins. " << std::endl;
+      auto active_list_v_w = active_relative_neighbors(v, w);
+      auto active_list_w_v = active_relative_neighbors(w, v);
+      if (active_list_w_v.size() <
+          active_list_v_w.size()) {  // simulate the contraction of w by expanding the star of v
+        for (auto& x : active_list_w_v) {
+          active_edge_insertion(v, x, filt_val);
+          // std::cout << "Inserted the edge " << v << " , " << x  << std::endl;
+        }
+        swap_rows(v, w);
+        // std::cout << "A swap of the vertex " << v << " and " << w << "took place." << std::endl;
+      } else {
+        for (auto& y : active_list_v_w) {
+          active_edge_insertion(w, y, filt_val);
+          // std::cout << "Inserted the edge " << w << ", " << y  << std::endl;
+        }
+      }
+      auto rw_v = vertexToRow.find(v);
+      contractionIndicator[rw_v->second] = true;
+    }
+    if (membership(v) && !membership(w)) {
+      relable(v, w);
+    }
+  }
+  void active_edge_insertion(const Vertex& v, const Vertex& w, double filt_val) {
+    insert_new_edges(v, w, filt_val);
+    // update_active_indicator(v,w);
+  }
+
+  void print_sparse_skeleton() { std::cout << sparseRowAdjMatrix << std::endl; }
+};
\ No newline at end of file
diff --git a/src/Collapse/include/gudhi/Rips_edge_list.h b/src/Collapse/include/gudhi/Rips_edge_list.h
new file mode 100644
index 00000000..b7c4dcff
--- /dev/null
+++ b/src/Collapse/include/gudhi/Rips_edge_list.h
@@ -0,0 +1,184 @@
+/*    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, Pawel Dlotko, Vincent Rouvreau Siddharth Pritam
+ *
+ *    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/>.
+ */
+
+#ifndef RIPS_edge_list_H_
+#define RIPS_edge_list_H_
+
+#include <gudhi/Debug_utils.h>
+#include <gudhi/graph_simplicial_complex.h>
+#include <boost/graph/adjacency_list.hpp>
+
+#include <iostream>
+#include <vector>
+#include <map>
+#include <string>
+#include <limits>  // for numeric_limits
+#include <utility>  // for pair<>
+
+
+namespace Gudhi {
+
+namespace rips_edge_list {
+
+/**
+ * \class Rips_complex
+ * \brief Rips complex data structure.
+ * 
+ * \ingroup rips_complex
+ * 
+ * \details
+ * The data structure is a one skeleton graph, or Rips graph, containing edges when the edge length is less or equal
+ * to a given threshold. Edge length is computed from a user given point cloud with a given distance function, or a
+ * distance matrix.
+ * 
+ * \tparam Filtration_value is the type used to store the filtration values of the simplicial complex.
+ */
+template<typename Filtration_value>
+class Rips_edge_list {
+ public:
+  /**
+   * \brief Type of the one skeleton graph stored inside the Rips complex structure.
+   */
+  // typedef typename boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS
+  // , boost::property < vertex_filtration_t, Filtration_value >
+  // , boost::property < edge_filtration_t, Filtration_value >> OneSkeletonGraph;
+
+ private:
+  typedef int Vertex_handle;
+
+ public:
+  /** \brief Rips_complex constructor from a list of points.
+   *
+   * @param[in] points Range of points.
+   * @param[in] threshold Rips value.
+   * @param[in] distance distance function that returns a `Filtration_value` from 2 given points.
+   * 
+   * \tparam ForwardPointRange must be a range for which `std::begin` and `std::end` return input iterators on a
+   * point.
+   *
+   * \tparam Distance furnishes `operator()(const Point& p1, const Point& p2)`, where
+   * `Point` is a point from the `ForwardPointRange`, and that returns a `Filtration_value`.
+   */
+  template<typename ForwardPointRange, typename Distance >
+  Rips_edge_list(const ForwardPointRange& points, Filtration_value threshold, Distance distance) {
+    compute_proximity_graph(points, threshold, distance);
+  }
+
+  /** \brief Rips_complex constructor from a distance matrix.
+   *
+   * @param[in] distance_matrix Range of distances.
+   * @param[in] threshold Rips value.
+   * 
+   * \tparam DistanceMatrix must have a `size()` method and on which `distance_matrix[i][j]` returns
+   * the distance between points \f$i\f$ and \f$j\f$ as long as \f$ 0 \leqslant i < j \leqslant
+   * distance\_matrix.size().\f$
+   */
+  template<typename DistanceMatrix>
+  Rips_edge_list(const DistanceMatrix& distance_matrix, Filtration_value threshold) {
+    compute_proximity_graph(boost::irange((size_t)0, distance_matrix.size()), threshold,
+                            [&](size_t i, size_t j){return distance_matrix[j][i];});
+  }
+
+  /** \brief Initializes the egde list (one skeleton) complex from the Rips graph 
+   *
+   * \tparam EdgeListForRips must meet `EdgeListForRips` concept.
+   * 
+   * @param[in] edges EdgeListForRips to be created.
+   * @param[in] dim_max graph expansion for Rips until this given maximal dimension.
+   * @exception std::invalid_argument In debug mode, if `edges.num_vertices()` does not return 0.
+   *
+   */
+  template <typename EdgeListForRips>
+  void create_edges(EdgeListForRips& edge_list) {
+    GUDHI_CHECK(edges.num_vertices() == 0,
+                std::invalid_argument("Rips_complex::create_complex - edge list is not empty"));
+
+    // sort the tuple (filteration_valuem, (v1,v2){edge}) 
+    //By default the sort is done on the first element, so here it's filteration value.
+    std::sort(edges.begin(),edges.end());
+    for(size_t i = 0; i < edges.size(); i++ )
+      edge_list.emplace_back(edges.at(i));
+
+  }
+
+ private:
+  /** \brief Computes the proximity graph of the points.
+   *
+   * If points contains n elements, the proximity graph is the graph with n vertices, and an edge [u,v] iff the
+   * distance function between points u and v is smaller than threshold.
+   *
+   * \tparam ForwardPointRange furnishes `.begin()` and `.end()`
+   * methods.
+   *
+   * \tparam Distance furnishes `operator()(const Point& p1, const Point& p2)`, where
+   * `Point` is a point from the `ForwardPointRange`, and that returns a `Filtration_value`.
+   */
+  template< typename ForwardPointRange, typename Distance >
+  void compute_proximity_graph(const ForwardPointRange& points, Filtration_value threshold,
+               Distance distance) {
+    edges.clear();
+    // Compute the proximity graph of the points.
+    // If points contains n elements, the proximity graph is the graph with n vertices, and an edge [u,v] iff the
+    // distance function between points u and v is smaller than threshold.
+    // --------------------------------------------------------------------------------------------
+    // Creates the vector of edges and its filtration values (returned by distance function)
+    Vertex_handle idx_u = 0;
+    for (auto it_u = std::begin(points); it_u != std::end(points); ++it_u, ++idx_u) {
+      Vertex_handle idx_v = idx_u + 1;
+      for (auto it_v = it_u + 1; it_v != std::end(points); ++it_v, ++idx_v) {
+        Filtration_value fil = distance(*it_u, *it_v);
+        if (fil <= threshold) {
+          edges.emplace_back(fil, idx_u, idx_v);
+        }
+      }
+    }
+
+    // --------------------------------------------------------------------------------------------
+    // Creates the proximity graph from edges and sets the property with the filtration value.
+    // Number of points is labeled from 0 to idx_u-1
+    // --------------------------------------------------------------------------------------------
+    // Do not use : rips_skeleton_graph_ = OneSkeletonGraph(...) -> deep copy of the graph (boost graph is not
+    // move-enabled)
+    // rips_skeleton_graph_.~OneSkeletonGraph();
+    // new(&rips_skeleton_graph_)OneSkeletonGraph(edges.begin(), edges.end(), edges_fil.begin(), idx_u);
+
+    // auto vertex_prop = boost::get(vertex_filtration_t(), rips_skeleton_graph_);
+
+    // using vertex_iterator = typename boost::graph_traits<OneSkeletonGraph>::vertex_iterator;
+    // vertex_iterator vi, vi_end;
+    // for (std::tie(vi, vi_end) = boost::vertices(rips_skeleton_graph_);
+    //      vi != vi_end; ++vi) {
+    //   boost::put(vertex_prop, *vi, 0.);
+    // }
+  }
+
+ private:
+  // OneSkeletonGraph rips_skeleton_graph_;
+  std::vector< std::tuple < Filtration_value, Vertex_handle, Vertex_handle > > edges;
+  // std::vector< Filtration_value > edges_fil;
+};
+
+}  // namespace rips_complex
+
+}  // namespace Gudhi
+
+#endif  // RIPS_COMPLEX_H_
diff --git a/src/Collapse/utilities/CMakeLists.txt b/src/Collapse/utilities/CMakeLists.txt
new file mode 100644
index 00000000..c742144b
--- /dev/null
+++ b/src/Collapse/utilities/CMakeLists.txt
@@ -0,0 +1,33 @@
+project(Collapse_utilities)
+
+# From a point cloud
+add_executable ( point_cloud_edge_collapse_rips_persistence point_cloud_edge_collapse_rips_persistence.cpp )
+target_link_libraries(point_cloud_edge_collapse_rips_persistence Boost::program_options)
+
+if (TBB_FOUND)
+  target_link_libraries(point_cloud_edge_collapse_rips_persistence ${TBB_LIBRARIES})
+endif()
+add_test(NAME Edge_collapse_utilities_point_cloud_rips_persistence COMMAND $<TARGET_FILE:point_cloud_edge_collapse_rips_persistence>
+         "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3" "-o" "off_results.pers")
+
+install(TARGETS point_cloud_edge_collapse_rips_persistence DESTINATION bin)
+
+# From a distance matrix
+add_executable ( distance_matrix_edge_collapse_rips_persistence distance_matrix_edge_collapse_rips_persistence.cpp )
+target_link_libraries(distance_matrix_edge_collapse_rips_persistence Boost::program_options)
+
+if (TBB_FOUND)
+  target_link_libraries(distance_matrix_edge_collapse_rips_persistence ${TBB_LIBRARIES})
+endif()
+add_test(NAME Edge_collapse_utilities_distance_matrix_rips_persistence COMMAND $<TARGET_FILE:distance_matrix_edge_collapse_rips_persistence>
+         "${CMAKE_SOURCE_DIR}/data/distance_matrix/tore3D_1307_distance_matrix.csv" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3" "-o" "csv_results.pers")
+
+install(TARGETS distance_matrix_edge_collapse_rips_persistence DESTINATION bin)
+
+if (DIFF_PATH)
+  add_test(Edge_collapse_utilities_diff_persistence ${DIFF_PATH}
+           "off_results.pers" "csv_results.pers")
+  set_tests_properties(Edge_collapse_utilities_diff_persistence PROPERTIES DEPENDS
+                       "Edge_collapse_utilities_point_cloud_rips_persistence;Edge_collapse_utilities_distance_matrix_rips_persistence")
+         
+endif()
diff --git a/src/Collapse/utilities/distance_matrix_edge_collapse_rips_persistence.cpp b/src/Collapse/utilities/distance_matrix_edge_collapse_rips_persistence.cpp
new file mode 100644
index 00000000..63c91ebc
--- /dev/null
+++ b/src/Collapse/utilities/distance_matrix_edge_collapse_rips_persistence.cpp
@@ -0,0 +1,221 @@
+#include <gudhi/FlagComplexSpMatrix.h>
+#include <gudhi/Rips_complex.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Rips_edge_list.h>
+#include <gudhi/distance_functions.h>
+#include <gudhi/reader_utils.h>
+#include <gudhi/Points_off_io.h>
+
+#include <CGAL/Epick_d.h>
+
+#include <boost/program_options.hpp>
+
+// Types definition
+using Point = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>::Point_d;
+using Vector_of_points = std::vector<Point>;
+
+using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
+using Filtration_value = double;
+using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>;
+using Rips_edge_list = Gudhi::rips_edge_list::Rips_edge_list<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* const argv[], double& min_persistence, double& end_thresold,
+                     int& dimension, int& dim_max, std::string& csv_matrix_file, std::string& filediag) {
+  namespace po = boost::program_options;
+  po::options_description visible("Allowed options", 100);
+      visible.add_options()
+        ("help,h", "produce help message")
+      	("min_persistence,m", po::value<double>(&min_persistence)->default_value(0.1),
+         "Minimum persistence interval length")
+        ("end_thresold,e", po::value<double>(&end_thresold)->default_value(1),
+         "Final threshold for rips complex.")
+        ("dimensions,D", po::value<int>(&dimension)->default_value(2),
+         "Dimension of the manifold.")
+        ("dim_max,k ", po::value<int>(&dim_max)->default_value(2),
+         "Maximum allowed dimension of the Rips complex.")
+        ("input_file_name,i", po::value<std::string>(&csv_matrix_file),
+         "The input file.")
+        ("filediag,o", po::value<std::string>(&filediag),
+         "The output file.");
+
+  po::options_description all;
+  all.add(visible);
+  po::variables_map vm;
+  po::store(po::command_line_parser(argc, argv).options(all).run(), vm);
+  po::notify(vm);
+  if (vm.count("help")) {
+    std::cout << std::endl;
+    std::cout << "Computes rips complexes of different threshold values, to 'end_thresold' from a n random uniform "
+                 "point_vector on a selected manifold, . \n";
+    std::cout << "Strongly collapses all the rips complexes and output the results in out_file. \n";
+    std::cout << "The experiments are repeted 'repete' num of times for each threshold value. \n";
+    std::cout << "type -m for manifold options, 's' for uni sphere, 'b' for unit ball, 'f' for file. \n";
+    std::cout << "type -i 'filename' for Input file option for exported point sample. \n";
+    std::cout << std::endl << std::endl;
+    std::cout << "Usage: " << argv[0] << " [options]" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
+  }
+}
+
+class filt_edge_to_dist_matrix {
+ public:
+  template <class Distance_matrix, class Filtered_sorted_edge_list>
+  filt_edge_to_dist_matrix(Distance_matrix& distance_mat, Filtered_sorted_edge_list& edge_filt,
+                           std::size_t number_of_points) {
+    double inf = std::numeric_limits<double>::max();
+    doubleVector distances;
+    std::pair<std::size_t, std::size_t> e;
+    for (std::size_t indx = 0; indx < number_of_points; indx++) {
+      for (std::size_t j = 0; j <= indx; j++) {
+        if (j == indx)
+          distances.push_back(0);
+        else
+          distances.push_back(inf);
+      }
+      distance_mat.push_back(distances);
+      distances.clear();
+    }
+
+    for (auto edIt = edge_filt.begin(); edIt != edge_filt.end(); edIt++) {
+      e = std::minmax(std::get<1>(*edIt), std::get<2>(*edIt));
+      distance_mat.at(std::get<1>(e)).at(std::get<0>(e)) = std::get<0>(*edIt);
+    }
+  }
+};
+
+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[]) {
+  auto the_begin = std::chrono::high_resolution_clock::now();
+
+  typedef size_t Vertex_handle;
+  typedef std::vector<std::tuple<Filtration_value, Vertex_handle, Vertex_handle>> Filtered_sorted_edge_list;
+
+  std::string csv_matrix_file;
+  std::string filediag;
+  double threshold;
+  int dim_max = 2;
+  int p;
+  double min_persistence;
+
+  program_options(argc, argv, csv_matrix_file, filediag, threshold, dim_max, p, min_persistence);
+
+  Map map_empty;
+
+  Distance_matrix distances;
+  Distance_matrix sparse_distances;
+
+  distances = Gudhi::read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_matrix_file);
+  std::size_t number_of_points = distances.size();
+  std::cout << "Read the distance matrix succesfully, of size: " << number_of_points << std::endl;
+
+  Filtered_sorted_edge_list edge_t;
+  std::cout << "Computing the one-skeleton for threshold: " << threshold << std::endl;
+
+  Rips_edge_list Rips_edge_list_from_file(distances, threshold);
+  Rips_edge_list_from_file.create_edges(edge_t);
+  std::cout<< "Sorted edge list computed" << std::endl;
+  
+  if (edge_t.size() <= 0) {
+    std::cerr << "Total number of egdes are zero." << std::endl;
+    exit(-1);
+  }
+
+  std::cout << "Total number of edges before collapse are: " << edge_t.size() << std::endl;
+
+  // Now we will perform filtered edge collapse to sparsify the edge list edge_t.
+  std::cout << "Filtered edge collapse begins" << std::endl;
+  FlagComplexSpMatrix mat_filt_edge_coll(number_of_points, edge_t);
+  std::cout << "Matrix instansiated" << std::endl;
+  Filtered_sorted_edge_list collapse_edges;
+  collapse_edges = mat_filt_edge_coll.filtered_edge_collapse();
+  filt_edge_to_dist_matrix(sparse_distances, collapse_edges, number_of_points);
+  std::cout << "Total number of vertices after collapse in the sparse matrix are: " << mat_filt_edge_coll.num_vertices()
+            << std::endl;
+
+  Rips_complex rips_complex_after_collapse(sparse_distances, threshold);
+
+  Simplex_tree stree;
+  rips_complex_after_collapse.create_complex(stree, dim_max);
+
+  std::cout << "The complex contains " << stree.num_simplices() << " simplices  after collapse. \n";
+  std::cout << "   and has dimension " << stree.dimension() << " \n";
+
+  // Sort the simplices in the order of the filtration
+  stree.initialize_filtration();
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(stree);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(3);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+  if (filediag.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::ofstream out(filediag);
+    pcoh.output_diagram(out);
+    out.close();
+  }
+
+  auto the_end = std::chrono::high_resolution_clock::now();
+
+  std::cout << "Total computation time : " << std::chrono::duration<double, std::milli>(the_end - the_begin).count()
+            << " ms\n"
+            << std::endl;
+  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 after edge collapse 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;
+    exit(-1);
+  }
+}
diff --git a/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp b/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp
new file mode 100644
index 00000000..af08477c
--- /dev/null
+++ b/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp
@@ -0,0 +1,205 @@
+#include <gudhi/FlagComplexSpMatrix.h>
+#include <gudhi/Rips_complex.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Rips_edge_list.h>
+#include <gudhi/distance_functions.h>
+#include <gudhi/reader_utils.h>
+#include <gudhi/Points_off_io.h>
+
+#include <CGAL/Epick_d.h>
+
+#include <boost/program_options.hpp>
+
+// Types definition
+using Point = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>::Point_d;
+using Vector_of_points = std::vector<Point>;
+
+using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
+using Filtration_value = double;
+using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>;
+using Rips_edge_list = Gudhi::rips_edge_list::Rips_edge_list<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>>;
+
+
+class filt_edge_to_dist_matrix {
+ public:
+  template <class Distance_matrix, class Filtered_sorted_edge_list>
+  filt_edge_to_dist_matrix(Distance_matrix& distance_mat, Filtered_sorted_edge_list& edge_filt,
+                           std::size_t number_of_points) {
+    double inf = std::numeric_limits<double>::max();
+    doubleVector distances;
+    std::pair<std::size_t, std::size_t> e;
+    for (std::size_t indx = 0; indx < number_of_points; indx++) {
+      for (std::size_t j = 0; j <= indx; j++) {
+        if (j == indx)
+          distances.push_back(0);
+        else
+          distances.push_back(inf);
+      }
+      distance_mat.push_back(distances);
+      distances.clear();
+    }
+
+    for (auto edIt = edge_filt.begin(); edIt != edge_filt.end(); edIt++) {
+      e = std::minmax(std::get<1>(*edIt), std::get<2>(*edIt));
+      distance_mat.at(std::get<1>(e)).at(std::get<0>(e)) = std::get<0>(*edIt);
+    }
+  }
+};
+
+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[]) {
+  typedef size_t Vertex_handle;
+  typedef std::vector<std::tuple<Filtration_value, Vertex_handle, Vertex_handle>> Filtered_sorted_edge_list;
+
+  auto the_begin = std::chrono::high_resolution_clock::now();
+  std::size_t number_of_points;
+
+
+  std::string off_file_points;
+  std::string filediag;
+  double threshold;
+  int dim_max;
+  int p;
+  double min_persistence;
+
+  program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence);
+
+  std::cout << "The current input values to run the program is: " << std::endl;
+  std::cout << "min_persistence, threshold, max_complex_dimension, off_file_points, filediag"
+            << std::endl;
+  std::cout << min_persistence << ", " << threshold << ", " << dim_max
+            << ", " << off_file_points << ", " << filediag << std::endl;
+
+  Map map_empty;
+
+  Distance_matrix sparse_distances;
+
+  Gudhi::Points_off_reader<Point> off_reader(off_file_points);
+  if (!off_reader.is_valid()) {
+    std::cerr << "Unable to read file " << off_file_points << "\n";
+    exit(-1);  // ----- >>
+  }
+
+  Vector_of_points point_vector = off_reader.get_point_cloud();
+  if (point_vector.size() <= 0) {
+    std::cerr << "Empty point cloud." << std::endl;
+    exit(-1);  // ----- >>
+  }
+
+  int dimension = point_vector[0].dimension();
+  number_of_points = point_vector.size();
+  std::cout << "Successfully read " << number_of_points << " point_vector.\n";
+  std::cout << "Ambient dimension is " << dimension << ".\n";
+
+  std::cout << "Point Set Generated." << std::endl;
+
+  Filtered_sorted_edge_list edge_t;
+  std::cout << "Computing the one-skeleton for threshold: " << threshold << std::endl;
+
+  Rips_edge_list Rips_edge_list_from_file(point_vector, threshold, Gudhi::Euclidean_distance());
+  Rips_edge_list_from_file.create_edges(edge_t);
+
+  std::cout << "Sorted edge list computed" << std::endl;
+  std::cout << "Total number of edges before collapse are: " << edge_t.size() << std::endl;
+
+  if (edge_t.size() <= 0) {
+    std::cerr << "Total number of egdes are zero." << std::endl;
+    exit(-1);
+  }
+
+  // Now we will perform filtered edge collapse to sparsify the edge list edge_t.
+  std::cout << "Filtered edge collapse begins" << std::endl;
+  FlagComplexSpMatrix mat_filt_edge_coll(number_of_points, edge_t);
+  std::cout << "Matrix instansiated" << std::endl;
+  Filtered_sorted_edge_list collapse_edges;
+  collapse_edges = mat_filt_edge_coll.filtered_edge_collapse();
+  filt_edge_to_dist_matrix(sparse_distances, collapse_edges, number_of_points);
+  std::cout << "Total number of vertices after collapse in the sparse matrix are: " << mat_filt_edge_coll.num_vertices()
+            << std::endl;
+
+  // Rips_complex rips_complex_before_collapse(distances, threshold);
+  Rips_complex rips_complex_after_collapse(sparse_distances, threshold);
+
+  Simplex_tree stree;
+  rips_complex_after_collapse.create_complex(stree, dim_max);
+
+  std::cout << "The complex contains " << stree.num_simplices() << " simplices  after collapse. \n";
+  std::cout << "   and has dimension " << stree.dimension() << " \n";
+
+  // Sort the simplices in the order of the filtration
+  stree.initialize_filtration();
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(stree);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(p);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+  if (filediag.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::ofstream out(filediag);
+    pcoh.output_diagram(out);
+    out.close();
+  }
+
+  auto the_end = std::chrono::high_resolution_clock::now();
+
+  std::cout << "Total computation time : " << std::chrono::duration<double, std::milli>(the_end - the_begin).count()
+            << " ms\n"
+            << std::endl;
+  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, after edge collapse, 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;
+    exit(-1);
+  }
+}
\ No newline at end of file