Collapse first version from Siddharth

7 files changed, 2000 insertions, 0 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0b5f5144..8b82036b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -25,6 +25,7 @@ add_gudhi_module(common)
 add_gudhi_module(Alpha_complex)
 add_gudhi_module(Bitmap_cubical_complex)
 add_gudhi_module(Bottleneck_distance)
+add_gudhi_module(Collapse)
 add_gudhi_module(Contraction)
 add_gudhi_module(Cech_complex)
 add_gudhi_module(Hasse_complex)
diff --git a/src/Collapse/example/CMakeLists.txt b/src/Collapse/example/CMakeLists.txt
new file mode 100644
index 00000000..c9cba3fa
--- /dev/null
+++ b/src/Collapse/example/CMakeLists.txt
@@ -0,0 +1,9 @@
+project(Collapse_examples)
+
+add_executable ( Collapse_example_rips_persistence_with_collapse rips_persistence_with_sc.cpp )
+target_link_libraries(Collapse_example_rips_persistence_with_collapse Boost::program_options)
+
+if (TBB_FOUND)
+  target_link_libraries(Collapse_example_rips_persistence_with_collapse ${TBB_LIBRARIES})
+endif()
+add_test(NAME Collapse_example_rips_persistence_with_collapse COMMAND $<TARGET_FILE:Collapse_example_rips_persistence_with_collapse>)
diff --git a/src/Collapse/example/rips_persistence_with_sc.cpp b/src/Collapse/example/rips_persistence_with_sc.cpp
new file mode 100644
index 00000000..5336b202
--- /dev/null
+++ b/src/Collapse/example/rips_persistence_with_sc.cpp
@@ -0,0 +1,386 @@
+#include <gudhi/TowerAssembler_FlagComplex.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/PointSetGen.h>
+
+// Types definition
+using Vector_of_points         = std::vector<Point>;
+using Vector_of_SM_pointers    = std::vector<FlagComplexSpMatrix*>;
+
+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 extract_sub_one_skeleton
+{
+public:
+    template<class Filtered_sorted_edge_list,  class Fil_vector >
+    extract_sub_one_skeleton(double threshold, Filtered_sorted_edge_list & current_edge_t, Filtered_sorted_edge_list & edge_t, Fil_vector & edge_filt ) {
+ 
+        auto end_it = std::upper_bound(edge_filt.begin(), edge_filt.end(), threshold); // find_index(edge_t, threshold, 0, end_idx);
+        size_t end_idx = std::distance(edge_filt.begin(), end_it);
+
+        for( size_t idx = 0; idx < end_idx ; idx++) {
+           current_edge_t.push_back(*edge_t.begin()); 
+           edge_filt.erase(edge_filt.begin());
+           edge_t.erase(edge_t.begin());
+        }
+        
+    }
+};
+class extract_one_new_edge
+{
+public:
+    template<class Filtered_sorted_edge_list,  class Fil_vector >
+    extract_one_new_edge(Filtered_sorted_edge_list & current_edge_t, Filtered_sorted_edge_list & edge_t, Fil_vector & edge_filt ) {
+        current_edge_t.push_back(*edge_t.begin()); 
+        edge_filt.erase(edge_filt.begin());
+        edge_t.erase(edge_t.begin());
+      
+    }
+};
+class strong_filtered_vertex_collapse
+{
+public:
+    template<class Filtered_sorted_edge_list, class Distance_matrix>
+    strong_filtered_vertex_collapse(std::size_t number_of_points, Filtered_sorted_edge_list & edge_t, double steps, Distance_matrix & sparse_distances)
+    {
+         auto the_collapse_begin = std::chrono::high_resolution_clock::now();
+        //An additional vector <edge_filt> to perform binary search to find the index of given threshold
+        std::vector<Filtration_value> * edge_filt = new std::vector<Filtration_value>();
+        edge_filt->clear();
+        for(auto edIt = edge_t.begin(); edIt != edge_t.end(); edIt++) {
+            edge_filt->push_back(std::get<0>(*edIt));
+        }
+        double begin_thresold = edge_filt->at(0);
+        double end_threshold = edge_filt->at(edge_filt->size()-1);
+        double totAssembleTime = 0.0;
+        Filtered_sorted_edge_list * sub_skeleton  = new Filtered_sorted_edge_list();  
+        TowerAssembler_FlagComplex twr_assembler(number_of_points) ;
+
+        std::cout<< "Begin and end thresholds (filteration values) are , " << begin_thresold <<", " << end_threshold << "." <<std::endl;
+
+        // int iterations = (end_threshold - begin_thresold)/steps;
+        // std::cout << "Total number of iterations to be run are: " << iterations << std::endl;
+       
+        auto threshold =  begin_thresold;  
+
+        FlagComplexSpMatrix * mat_coll       = new FlagComplexSpMatrix(); 
+        FlagComplexSpMatrix * mat_prev_coll  = new FlagComplexSpMatrix(number_of_points); 
+
+        std::cout << "Going for vertex collapse and tower assembly" << std::endl;
+
+        int i = 1;
+        Map * redmap;
+        while(edge_filt->size() > 0) {
+           
+            if(steps > 0)
+                extract_sub_one_skeleton(threshold, *sub_skeleton, edge_t, *edge_filt);
+
+            else
+                extract_one_new_edge(*sub_skeleton, edge_t, *edge_filt);
+
+
+            mat_coll = new FlagComplexSpMatrix(number_of_points, *sub_skeleton);
+            // mat_coll->strong_vertex_edge_collapse();
+            mat_coll->strong_vertex_collapse();
+            // std::cout<< " Total number of current edges are " << sub_skeleton->size() << std::endl;
+                       
+            redmap = new Map();
+            *redmap = mat_coll->reduction_map(); 
+            
+            // std::cout << "Subcomplex #" << i << " of threshold "<< threshold << " Collapsed" << std::endl;
+            totAssembleTime += twr_assembler.build_tower_for_two_cmplxs(*mat_prev_coll, *mat_coll, *redmap, threshold, "./PersistenceOutput/CollapsedTowerRips.txt");
+            // std::cout << "Tower updated for subcomplex #" << i << std::endl; 
+            
+            delete mat_prev_coll;
+            mat_prev_coll = new FlagComplexSpMatrix();
+            mat_prev_coll = mat_coll;
+            mat_coll  = new FlagComplexSpMatrix();
+            i++;
+            delete redmap;
+            if(steps > 0) {
+                threshold = threshold+steps;
+                if(threshold > end_threshold)
+                    threshold = end_threshold; 
+            }
+            else {
+                if(edge_filt->size() > 0)
+                    threshold = edge_filt->at(0);  
+            }    
+        }
+        std::cout << "Tower updated for subcomplex #" << i << std::endl; 
+        auto the_collapse = std::chrono::high_resolution_clock::now();    
+        sparse_distances = twr_assembler.distance_matrix();
+        std::cout << "Collapse And assembly time : " <<  std::chrono::duration<double, std::milli>(the_collapse- the_collapse_begin).count() << " ms\n" << std::endl;
+    }    
+
+};
+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 (int 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);
+        } 
+    }
+};
+
+
+int main(int argc, char * const argv[]) {
+	
+    auto the_begin = std::chrono::high_resolution_clock::now();
+    PointSetGen point_generator;
+    std::string out_file_name   = "default";
+    std::string in_file_name    = "default";
+    std::size_t number_of_points;
+    
+    typedef size_t Vertex_handle;
+    typedef std::vector< std::tuple<Filtration_value, Vertex_handle, Vertex_handle > > Filtered_sorted_edge_list;
+
+    int     dimension;
+    double  begin_thresold;
+    double  end_threshold;
+    double  steps;
+    int     repetetions = 1;
+    char    manifold;
+
+    Vector_of_points * point_vector;
+    Vector_of_points file_all_points;
+
+    std::string manifold_full = "sphere";
+    
+    double radius  = 1;
+    double r_min  = 0.6;
+    double r_max = 1;
+    int dim_max  = 2;
+
+    point_generator.program_options(argc, argv, number_of_points, begin_thresold, steps, end_threshold, repetetions, manifold, dimension, dim_max, in_file_name, out_file_name);
+    
+    std::cout << "The current input values to run the program is: "<< std::endl;
+    std::cout << "number_of_points, begin_thresold, steps, end_threshold, repetetions, manifold, dimension, max_complex_dimension, in_file_name, out_file_name" << std::endl;
+    std::cout << number_of_points << ", " << begin_thresold << ", " << steps << ", " << end_threshold << ", " << repetetions << ", " << manifold << ", " << dimension << ", " << dim_max << ", " << in_file_name << ", " << out_file_name << std::endl;
+    
+    if(manifold == 'f' || manifold =='F') {
+        Gudhi::Points_off_reader<Point> off_reader(in_file_name);
+        if (!off_reader.is_valid()) {
+            std::cerr << "Unable to read file " << in_file_name << "\n";
+            exit(-1);  // ----- >>
+        }
+
+        file_all_points = Vector_of_points(off_reader.get_point_cloud());
+        dimension = file_all_points[0].dimension() ;
+        std::cout << "Successfully read " << file_all_points.size() << " point_vector.\n";
+        std::cout << "Ambient dimension is " << dimension << ".\n";
+    }
+   
+    Map map_empty;
+
+    std::string origFile ("./PersistenceOutput/original_tower_rips" );
+    std::string collFile  ("./PersistenceOutput/collapsed_tower_rips") ;
+    
+    std::string filediag_bfr ("./PersistenceOutput/uncollapsed_persistence_diags") ;
+    std::string filediag_aft ("./PersistenceOutput/collapsed_persistence_diags") ;
+    
+     std::string origPoints ("./PersistenceOutput/pointsamaple.off");
+    // std::string otherStats ("./PersistenceOutput/maximal_simplx_cnt");
+    // otherStats = otherStats+"_"+ out_file_name+ ".txt";
+    filediag_bfr = filediag_bfr+"_"+ out_file_name+ ".txt";
+    filediag_aft = filediag_aft+"_"+ out_file_name+ ".txt";
+
+    double currentCreationTime = 0.0;
+    double maxCreationTime     = 0.0;
+
+   
+    point_vector = new Vector_of_points();
+    Distance_matrix distances;
+    Distance_matrix *sparse_distances = new Distance_matrix();
+
+
+    if(manifold == 's' || manifold == 'S'){
+        // point_generator.generate_points_sphere(*point_vector, number_of_points, dimension, radius); 
+        point_generator.generate_grid_2sphere(*point_vector, number_of_points, radius); 
+        origFile = origFile+"_sphere_"+out_file_name+".txt";
+        collFile = collFile+"_sphere_"+out_file_name+".txt";
+        std::cout << number_of_points << " points successfully chosen randomly from "<< dimension <<"-sphere of radius " << radius << std::endl;
+    }
+    else if(manifold == 'b' || manifold == 'B'){
+        point_generator.generate_points_ball(*point_vector, number_of_points, dimension, radius); 
+        origFile = origFile+"_ball_"+out_file_name+".txt";
+        collFile = collFile+"_ball_"+out_file_name+".txt";
+        std::cout << number_of_points << " points successfully chosen randomly from "<< dimension <<"-ball of radius " << radius << std::endl;
+    
+    }
+    else if( (manifold == 'a' || manifold == 'A')&& dimension == 2){
+        point_generator.generate_points_2annulus(*point_vector, number_of_points, r_min, r_max); 
+        origFile = origFile+"_annulus_"+out_file_name+".txt";
+        collFile = collFile+"_annulus_"+out_file_name+".txt";
+        std::cout << number_of_points << " points successfully chosen randomly from "<< 2 <<"-annulus of radii (" << r_min << ',' << r_max << ") " << std::endl;
+    }
+    else if( (manifold == 'a' || manifold == 'A') && dimension == 3){
+        point_generator.generate_points_spherical_shell(*point_vector, number_of_points, r_min, r_max); 
+        origFile = origFile+"_annulus_"+out_file_name+".txt";
+        collFile = collFile+"_annulus_"+out_file_name+".txt";
+        std::cout << number_of_points << " points successfully chosen randomly from spherical shell of radii (" << r_min << ',' << r_max << ") " << std::endl;
+    }
+    
+    else if(manifold == 'f' || manifold =='f') {
+        // Subsampling from all points for each iterations
+        Gudhi::subsampling::pick_n_random_points(file_all_points, number_of_points, std::back_inserter(*point_vector));
+        origFile = origFile+"_"+ out_file_name+ ".txt";
+        collFile = collFile+"_"+ out_file_name+ ".txt";
+        std::cout << number_of_points << " points succesfully chosen randomly of dimension "<< dimension << " ." << std::endl;
+    }
+    else if (manifold == 'm'){
+        std::string csv_file_name(in_file_name);
+        distances = Gudhi::read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_file_name);
+        number_of_points = distances.size();
+        std::cout << "Read the distance matrix succesfully, of size: " << number_of_points << std::endl;
+        origFile = origFile+"_"+ out_file_name+ ".txt";
+        collFile = collFile+"_"+ out_file_name+ ".txt";
+    }
+    else {
+        std::cerr << "Wrong parameters for input manifold..." <<std::endl;  
+        exit(-1); 
+    }
+    
+    if( manifold != 'm')
+        number_of_points = point_vector->size();
+    std::cout << "Point Set Generated."  <<std::endl;
+ 
+    // for(int i = 0; i < number_of_points; i++ )
+    //     point_generator.print_point(point_vector->at(i));
+    // point_generator.output_points(*point_vector, origPoints);
+   
+    Filtered_sorted_edge_list * edge_t = new Filtered_sorted_edge_list();
+    std::cout << "Computing the one-skeleton for threshold: " << end_threshold << std::endl; 
+    
+    auto begin_full_cmplx = std::chrono::high_resolution_clock::now();
+    if(manifold == 'm') {                                           //Input is a distance 'm'atrix //Creating the edge list
+        Rips_edge_list Rips_edge_list_from_file(distances, end_threshold);
+        Rips_edge_list_from_file.create_edges(*edge_t);
+        std::cout<< "Sorted edge list computed" << std::endl;
+
+        //Creating the Rips Complex
+        //Rips_complex rips_complex_from_file(distances, end_threshold);
+        //rips_complex_from_file.create_complex(*subComplex, dim_max);
+        //std::cout<< "Rips complex computed" << std::endl;
+    }
+    else{ //Point cloud input  //Creating the edge list
+        Rips_edge_list Rips_edge_list_from_points(*point_vector, end_threshold, Gudhi::Euclidean_distance());
+        Rips_edge_list_from_points.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;
+
+        // Creating the Rips Complex
+        // Rips_complex rips_complex_after_collapse(*point_vector, end_threshold, Gudhi::Euclidean_distance());
+        // rips_complex_from_points.create_complex(*subComplex, dim_max);
+        // std::cout<< "Rips complex computed" << 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  = new FlagComplexSpMatrix(number_of_points,*edge_t,true);
+    std::cout<< "Matrix instansiated" << std::endl;
+    if(edge_t->size() >0){
+        delete edge_t;
+        edge_t = new Filtered_sorted_edge_list();
+        *edge_t = mat_filt_edge_coll->filtered_edge_collapse(0); 
+        filt_edge_to_dist_matrix(*sparse_distances, *edge_t, number_of_points); 
+        std::cout << "Total number of vertices after collapse in the sparse matrix are: " << mat_filt_edge_coll->num_vertices() << std::endl;
+    }
+    else
+    {
+        std::cerr << "Total number of egdes are zero." <<std::endl;
+        exit(-1) ;
+    }
+    
+
+    // std::cout<< "Vertex strong-collapse begins" << std::endl;
+    // strong_filtered_vertex_collapse(number_of_points, *edge_t, steps, *sparse_distances); 
+
+    auto end_full_cmplx = std::chrono::high_resolution_clock::now();
+    currentCreationTime = std::chrono::duration<double, std::milli>(end_full_cmplx - begin_full_cmplx).count();
+    maxCreationTime = currentCreationTime;
+   
+ 
+    // Rips_complex rips_complex_before_collapse(distances, end_threshold);
+    Rips_complex rips_complex_after_collapse(*sparse_distances, end_threshold);
+     // Rips_complex rips_complex_after_collapse(*point_vector, end_threshold, Gudhi::Euclidean_distance());
+    // Construct the Rips complex in a Simplex Tree
+    
+    Simplex_tree   simplex_tree_aft;
+    // Simplex_tree   simplex_tree_bfr;
+    // rips_complex_before_collapse.create_complex(simplex_tree_bfr, dim_max);
+    rips_complex_after_collapse.create_complex(simplex_tree_aft, dim_max);
+
+    // std::cout << "The complex contains " << simplex_tree_bfr.num_simplices() << " simplices before collapse. \n";
+    // std::cout << "   and has dimension " << simplex_tree_bfr.dimension() << " \n";
+
+    std::cout << "The complex contains " << simplex_tree_aft.num_simplices() << " simplices  after collapse. \n";
+    std::cout << "   and has dimension " << simplex_tree_aft.dimension() << " \n";
+
+    // Sort the simplices in the order of the filtration
+    // simplex_tree_bfr.initialize_filtration();
+    simplex_tree_aft.initialize_filtration();
+    // Compute the persistence diagram of the complex
+    // Persistent_cohomology pcoh_bfr(simplex_tree_bfr);
+    Persistent_cohomology pcoh_aft(simplex_tree_aft);
+    // initializes the coefficient field for homology
+    // pcoh_bfr.init_coefficients(2);
+    pcoh_aft.init_coefficients(2);
+
+    // pcoh_bfr.compute_persistent_cohomology(steps);
+    pcoh_aft.compute_persistent_cohomology(steps);
+    // Output the diagram in filediag
+    // if (filediag_bfr.empty()) {
+    //     pcoh_bfr.output_diagram();
+    // } 
+    // else {
+    //     std::ofstream out(filediag_bfr);
+    //     pcoh_bfr.output_diagram(out);
+    //     out.close();
+    //   }
+
+    if (filediag_aft.empty()) {
+        pcoh_aft.output_diagram();
+    } 
+    else {
+        std::ofstream out(filediag_aft);
+        pcoh_aft.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;
+
+}
+  
\ No newline at end of file
diff --git a/src/Collapse/include/gudhi/FlagComplexSpMatrix.h b/src/Collapse/include/gudhi/FlagComplexSpMatrix.h
new file mode 100644
index 00000000..d3e58d4f
--- /dev/null
+++ b/src/Collapse/include/gudhi/FlagComplexSpMatrix.h
@@ -0,0 +1,963 @@
+/*    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 intVector 	   = std::vector<int>;
+using doubleVector 	   = std::vector<double>;
+using vertexVector     = std::vector<Vertex>;
+using boolVector       = std::vector<bool>;
+
+using doubleQueue 	   = std::queue<double>;
+using edgeQueue		   = std::queue<Edge>;
+
+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 activeIndicator; 	// active 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;
+
+  	//! Stores the indices-pair of the edges to-be checked for domination in the current iteration. 
+    /*!
+      Initialised with all egdes for the first iteration.
+      Subsequently once a dominated row is found, its non-dominated neighbhour indices are inserted. // To be clarified.
+    */
+	//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 current status of the edges after the vertex-collapse has been performed. .
+    /*!
+      \code
+      u_edge_map = std::unordered_map<Edge, bool>
+      \endcode
+     The values an edge can take are true, false; 
+     true -> Inserted in filteredEgdeIter;
+     false -> Not inserted in filteredEgdeIter;
+    */
+	u_edge_map edgeStatusMap; 
+
+    //! 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. 
+	bool filtEdgeCol;
+	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;
+  		filtEdgeCol     = 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(int 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;
+          				edgeStatusMap[{rw,v}] = true;
+          			}
+        		}
+			}			
+		}
+		// 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)
+		    // cfiltVal = std::get<1>(fEgdeVector.at(endIdx));
+		    // std::cout << "The current processing index is " << endIdx << std::endl;
+
+	      	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>
+    */
+    //Filtered_sorted_edge_list * edge_t = new Filtered_sorted_edge_list();
+	FlagComplexSpMatrix(const size_t & num_vertices, const Filtered_sorted_edge_list  &edge_t) {
+		init();
+
+ 	  	filtEdgeCol = false; 
+ 		sparseRowAdjMatrix  = sparseRowMatrix(expansion_limit*num_vertices, expansion_limit*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);
+	  	}	
+  	 	sparseRowAdjMatrix.makeCompressed(); 
+       	
+ 	  	// std::cout << sparseRowAdjMatrix << std::endl;
+	}
+	FlagComplexSpMatrix(const size_t & num_vertices, const Filtered_sorted_edge_list  & edge_t,  const bool fEdgeCol) {
+		if(fEdgeCol)
+		{	
+			init();
+	 	  	filtEdgeCol = true;  
+	 		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(); 
+       	}
+       	else
+       		FlagComplexSpMatrix(num_vertices, edge_t);
+
+ 	  	// std::cout << sparseRowAdjMatrix << std::endl;
+	}
+
+	//!	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 of a given sparse-matrix(graph) without considering the filtration value.
+	// double strong_edge_collapse() {
+	// 	auto begin_collapse  = std::chrono::high_resolution_clock::now();
+	// 	critical_core_edges();
+	// 	vertexCollapsed = false;
+	// 	edgeCollapsed	= 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;
+	// 	//Post processing...
+	// 	after_edge_collapse();
+	// 	return collapseTime;
+	// }
+
+	
+	// Performs edge collapse in a decreasing sequence of the filtration value.
+	Filtered_sorted_edge_list filtered_edge_collapse(double steps) {
+		auto begin_collapse  = std::chrono::high_resolution_clock::now();
+		critical_core_edges();
+		vertexCollapsed = false;
+		edgeCollapsed	= 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 filtered edge Collapse : " << collapseTime << " ms\n" << std::endl;
+		//Post processing...
+		// after_edge_collapse();
+		// std::cout << sparseRowAdjMatrix << std::endl;
+		// for(auto idx = criticalCoreEdges.begin(); idx != criticalCoreEdges.end(); idx++ )
+		// {	
+		// 	std::cout << 
+		// }
+		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/PointSetGen.h b/src/Collapse/include/gudhi/PointSetGen.h
new file mode 100644
index 00000000..af78b63b
--- /dev/null
+++ b/src/Collapse/include/gudhi/PointSetGen.h
@@ -0,0 +1,273 @@
+#include <gudhi/Points_off_io.h>
+#include <gudhi/pick_n_random_points.h>
+
+#include <CGAL/Epick_d.h>
+#include <CGAL/point_generators_d.h>
+#include <CGAL/Random.h>
+
+#include <boost/program_options.hpp>
+
+using Point = CGAL::Epick_d< CGAL::Dynamic_dimension_tag>::Point_d;
+using Vector_of_points = std::vector<Point>;
+
+const double PI  = 3.141592653589793238463;
+#define _USE_MATH_DEFINES
+
+class PointSetGen {
+  private:
+  	double unirand(){return (double) rand()/(double) RAND_MAX;}
+  public:
+    void program_options(int argc, char * const argv[]
+                         , std::size_t & number_of_points
+                         , double & begin_thresold
+                         , double & steps
+                         , double & end_thresold
+                         , int 	  & repetetions
+                         , char   & manifold
+                         , int 	  & dimension
+                         , int    & dim_max
+                         , std::string & in_file_name
+                         , std::string & out_file_name
+                         ) {
+      namespace po = boost::program_options;
+
+      po::options_description visible("Allowed options", 100);
+      visible.add_options()
+        ("help,h", "produce help message")
+        ("number,n", po::value<std::size_t>(&number_of_points)->default_value(0),
+           "Number of generated point_vector.")
+        
+        ("begin_thresold,b", po::value<double>(&begin_thresold)->default_value(0),
+           "Initial threshold for rips complex.")
+      	("steps,s", po::value<double>(&steps)->default_value(0.1),
+       		"Steps of the threshold")
+        ("end_thresold,e", po::value<double>(&end_thresold)->default_value(1),
+      		"Final threshold for rips complex.")
+          
+        ("repetetions,r", po::value<int>(&repetetions)->default_value(1),
+        	"Num of repetetions of the experiments.")
+      	("manifold,m", po::value<char>(&manifold)->default_value('s'),
+       		"Type of manifold")
+           
+        ("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>(&in_file_name),
+         "The input file.")
+        ("out_file_name,o", po::value<std::string>(&out_file_name),
+         "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, from 'begin_thresold' to 'end_thresold', with priodic steps of 'steps' 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();
+      }
+    }
+
+   
+    void generate_points_sphere(Vector_of_points& W, int nbP, int dim, double radius) {
+    	CGAL::Random_points_on_sphere_d<Point> rp(dim+1, radius);
+    	for (int i = 0; i < nbP; i++)
+    		W.push_back(*rp++);
+    }
+    // void generate_fibonaci_grid_sphere(Vector_of_points& W, int nbP, int dim, double radius)
+    // {
+
+    // }
+    
+    void generate_grid_2sphere(Vector_of_points& W, int nbP, int r )
+    {
+    	std::vector<double> coords;
+    	int Ncount = 0;
+    	double p,v; //the angles phi and psi
+    	int M_p;
+		
+  		double a = (4*PI*pow(r,2))/nbP;
+  		double d = sqrt(a);
+  		
+  		int M_v = PI/d;
+  		
+  		double d_v = PI/M_v; 
+  		double d_p = a/d_v;
+  		
+  		for( int m = 0; m < M_v ; m++) {
+  			v = (PI*(m + 0.5))/M_v;
+  			M_p = ((2*PI*sin(v))/d_p);
+  			for(int n = 0;  n < M_p ; n++) {
+  				p = (2*PI*n)/M_p;
+  				coords = {r*sin(v)*cos(p), r*sin(v)*sin(p), r*cos(v)};
+  				W.push_back(Point(coords));
+  				Ncount += 1;
+  			}
+  		}
+    }
+
+    
+    /*
+      Generates point sets on <nbSpheres> spheres wedged at origin, sphere can have different radii from <init_radius> with steps of <multiplier_step>
+      Number of points on the sphere can also be different from <init_nbP> for the smallest sphere and then multiplied by <multiplier_step>
+    */
+    void generate_points_wedged_sphere(Vector_of_points& W, int init_nbP, int dim, double init_radius, int multiplier_step, int nbSpheres)  {
+      double radius = init_radius;
+      int nbP = init_nbP;
+      std::vector<double> translation;
+      for(int d = 0; d< dim; d++) {
+        translation.push_back(0);
+      }
+      for(int s = 0; s < nbSpheres; s++) {
+        CGAL::Random_points_on_sphere_d<Point> rp(dim+1, radius); 
+        for (int i = 0; i < nbP; i++) {
+          W.push_back(add_point(*rp++, translation, dim)); 
+        }
+        nbP = nbP*multiplier_step;
+        radius = radius*multiplier_step;
+        translation.at(dim-1) = (radius - init_radius);  
+      }
+    }
+
+    void generate_points_concentric_sphere(Vector_of_points& W, int init_nbP, int dim, int init_radius, int multiplier_step, int nbSpheres) {
+      double radius = init_radius;
+      int nbP = init_nbP;
+     
+      for(int s = 0; s < nbSpheres; s++) {
+        CGAL::Random_points_on_sphere_d<Point> rp(dim+1, radius); 
+        for (int i = 0; i < nbP; i++) {
+          W.push_back(*rp++); 
+        }
+        nbP = nbP*(pow(multiplier_step,2));
+        radius = radius*multiplier_step;
+      }
+      
+    }
+    void generate_points_2annulus(Vector_of_points& W, int nbP, double r_min, double r_max) {
+    	double rho, theta;
+        double x, y;
+	    std::vector<double> coords;
+	    double r_min_sq = pow(r_min,2);
+	    double r_max_sq = pow(r_max,2);
+	
+	    srand(time(NULL));
+	    for (int i=0; i<nbP; i++) {
+	        rho = sqrt((r_max_sq - r_min_sq)*unirand() + r_min_sq );
+	        theta = 2.*M_PI*unirand();
+	        x = rho*cos(theta);
+	        y = rho*sin(theta);
+
+	        coords = {x,y};
+	        W.push_back(Point(coords));
+	    }
+	}
+    void generate_points_spherical_shell(Vector_of_points& W, int nbP, double r_min, double r_max) {
+    	double rho, phi, theta;
+        double x, y, z;
+	    std::vector<double> coords;
+	    double r_min_cube = pow(r_min,3);
+	    double r_max_cube = pow(r_max,3);
+	
+	    srand(time(NULL));
+	    for (int i=0; i<nbP; i++) {
+	        rho = cbrt((r_max_cube - r_min_cube)*unirand() + r_min_cube );
+	        phi = 2.*M_PI*unirand();
+	        theta = acos(1. - 2.*unirand());
+	        	 
+	        x = rho*sin(theta)*cos(phi);
+	        y = rho*sin(theta)*sin(phi);
+	        z = rho*cos(theta);
+
+	        coords = {x,y,z};
+	        W.push_back(Point(coords));
+	    }
+	}   
+
+    void generate_points_ball(Vector_of_points& W, int nbP, int dim, double radius) {
+    	CGAL::Random_points_in_ball_d<Point> rp(dim, radius); 
+    	for (int i = 0; i < nbP; i++)
+    		W.push_back(*rp++);
+    }
+
+    void generate_points_cube(Vector_of_points& W, int nbP, int dim) {
+    	CGAL::Random_points_in_cube_d<Point> rp(dim, 6);
+    	for (int i = 0; i < nbP; i++)
+    		W.push_back(*rp++);
+    }
+
+    void add_point_vectors(Vector_of_points& V, Vector_of_points& U, int nbP, int dim) { // Adds two point vectors of the same size (nbP), by Modifying the first one, V = V+W.
+    	for (int i = 0; i < nbP; i++)
+    	{
+    		V[i] = add_point(V[i], U[i], dim);
+    	}
+    }
+
+    //returns x = x+y;
+    Point add_point(const Point & x, const Point & y, int dim) {
+      std::vector<double> coords;
+      for(int i =0; i< dim; i++)
+        coords.push_back(x[i]+y[i]); 
+      return Point(coords);
+    }
+    void print_point(const Point & x) { 
+      std::cout<< "(";
+      for(auto & p : x){
+        std::cout<< p << ", " ;
+      }
+      std::cout<< ")" << std::endl;  
+    }
+    void output_points(Vector_of_points & W, std::string outFile) {
+      std::ofstream myfile (outFile, std::ios::app);
+        if (myfile.is_open()) {    
+          myfile << "OFF"  << "  " << W.size() << " " << W.at(0).size() << std::endl;
+          for(auto & v : W){
+            for(auto & x : v){
+                myfile<< x << " " ;
+            }
+            myfile << std::endl;
+          }
+          myfile << "# Tower updated for the additional subcomplex.\n";
+          myfile.close();
+        }
+        else {
+            std::cerr << "Unable to open file";
+            exit(-1) ;
+        } 
+    }
+
+    Point noise_point(double noise_param, int dim, double radius) {
+	    std::vector<double> noise;
+	    for(int d = 0; d< dim; d++){
+	    	if(d % 2)
+	    		noise.push_back(-noise_param*radius);
+	    	else
+	    		noise.push_back(noise_param*radius);
+	    }
+	    return Point(noise);
+    }
+    //add noise to the points in W.
+    void add_noise(Vector_of_points& W, int nbP, int dim, double radius, double noise_param) {
+    	Point noise = noise_point(noise_param, dim, radius);
+    	for(Vector_of_points::iterator it = W.begin(); it != W.end(); it++ ) {
+    		*it = add_point(*it,noise,dim);
+    	}
+    }
+    
+    PointSetGen(){}
+    ~PointSetGen(){}
+};
\ 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/include/gudhi/TowerAssembler_FlagComplex.h b/src/Collapse/include/gudhi/TowerAssembler_FlagComplex.h
new file mode 100644
index 00000000..01d86642
--- /dev/null
+++ b/src/Collapse/include/gudhi/TowerAssembler_FlagComplex.h
@@ -0,0 +1,184 @@
+#pragma once
+#include <gudhi/FlagComplexSpMatrix.h>
+#include <gudhi/Rips_edge_list.h>
+
+#include <set>
+#include <fstream>
+#include <string>
+#include <algorithm>
+#include <limits>
+#include <cmath>
+
+
+typedef std::size_t Vertex;
+using Edge                  = std::pair<Vertex,Vertex>;
+using edge_list             = std::vector<Edge>;
+using Simplex               = std::vector<Vertex>;
+
+using vectorVertex          = std::vector<Vertex>;
+using vert_unSet            = std::unordered_set<Vertex>;
+using Map                   = std::unordered_map<Vertex,Vertex>;
+using Distance_matrix       = std::vector<std::vector<double>>;
+// using infinity              = std::numeric_limits<double>::max();
+// assumptions : (1) K1 and K2 have the same vertex set
+//               (2) The set of simplices of K1 is a subset of set of simplices of K2
+// K1  ->  K2    [Original Simplicial Complexes]
+// |       |
+// |       |
+// K1c ->  K2c   [Strongly Collapsed Flag Complexes]
+
+class TowerAssembler_FlagComplex
+{
+  private:
+	Map renamedVertices; 
+    Map representative_map;
+	std::size_t current_rename_counter;
+    Distance_matrix * distance_mat;
+    size_t total_vertices;
+
+    // Filtered_sorted_edge_list * edge_t = new Filtered_sorted_edge_list();
+    FlagComplexSpMatrix * flag_Filtration; 
+    typedef std::vector< std::tuple< double, Vertex, Vertex > > Filtered_sorted_edge_list;
+
+	
+  public:
+    
+    TowerAssembler_FlagComplex(std::size_t numVert)
+    {
+    	for (std::size_t i = 0; i <= numVert; ++i){
+    		renamedVertices[i] = i;
+        }
+    	total_vertices = 0;
+        current_rename_counter = numVert+1;
+        
+        flag_Filtration = new FlagComplexSpMatrix((numVert*log2(numVert)) +1);
+        distance_mat    = new Distance_matrix();
+        // distance_mat->push_back({0});
+    }
+    
+    ~TowerAssembler_FlagComplex(){};
+    double build_tower_for_two_cmplxs(FlagComplexSpMatrix mat_1, const FlagComplexSpMatrix & mat_2,  Map redmap_2,  double filtration_value, std::string outFile) // mat_1 and mat_2 are simplex_trees of K1c and K2c (the collapsed ones), redmap_2 is the map of K2 -> K2c
+    {
+        auto begin_print  = std::chrono::high_resolution_clock::now();
+        std::ofstream myfile (outFile, std::ios::app);
+        if (myfile.is_open())
+        {   
+            for (auto & v : mat_1.vertex_set()) {
+                auto collapsed_to = redmap_2.find(v);  // If v collapsed to something?
+                if(collapsed_to != redmap_2.end()) {  // Collapse happened, because there is a vertex in the map
+                    if(mat_1.membership(collapsed_to->second)) { // Collapsed to an existing vertex in mat_1.
+
+                    	// myfile << filtration_value  << " c " << renamedVertices.at(v) << " " << renamedVertices.at(collapsed_to->second) << std::endl; 
+                    	std::cout << filtration_value << " c " << renamedVertices.at(v) << " " << renamedVertices.at(collapsed_to->second) << std::endl;
+                        flag_Filtration->active_strong_expansion(renamedVertices.at(v), renamedVertices.at(collapsed_to->second), filtration_value);
+                        renamedVertices.at(v) = current_rename_counter;
+                        current_rename_counter++;                        
+                    }
+                    else {
+	                    // myfile << filtration_value << " i " << renamedVertices.at(collapsed_to->second) << std::endl;
+	                    // myfile  << filtration_value << " c " << renamedVertices.at(v) << " " << renamedVertices.at(collapsed_to->second) << std::endl; 
+	                    flag_Filtration->active_strong_expansion(renamedVertices.at(v), renamedVertices.at(collapsed_to->second),filtration_value);
+                        std::cout << filtration_value << " i " << renamedVertices.at(collapsed_to->second) << std::endl;
+                        std::cout  << filtration_value << " c " << renamedVertices.at(v) << " " << renamedVertices.at(collapsed_to->second) << std::endl; 
+                        renamedVertices.at(v) = current_rename_counter;
+	                    current_rename_counter++;
+                    }
+                     // If the vertex "collapsed_to->second" is not a member of mat_1, the contraction function will simply add and then collapse
+                     mat_1.contraction(v, collapsed_to->second);
+                     // std::cout << " Contraction Done " << std::endl;
+
+                }
+            }
+
+            //The core K1c (mat_1) has gone through the transformation(re-labeling)/collapse and it is now a subcomplex of K2c, the remaining simplices need to be included
+            // Writing the inclusion of all remaining simplices...
+            // std::cout << "Begining the inclusion of edges " << std::endl;
+            for( const Edge & e  : mat_2.all_edges()) {
+                auto u = std::get<0>(e);
+                auto v = std::get<1>(e);
+                // std::cout << "Going to insert the edge :" << renamedVertices.at(u) << ", " << renamedVertices.at(v) << std::endl;
+                // std::cout << "Going to insert the vertex :" << u << std::endl;
+                if(!mat_1.membership(u)) {
+                    flag_Filtration->insert_vertex(renamedVertices.at(u),filtration_value);
+                    // std::cout << "Inserted the vertex :" << renamedVertices.at(u) <<  " in the new distance matrix"<< std::endl;
+                    // myfile << filtration_value << " i";
+                    // myfile  << " " << renamedVertices.at(u);
+                    // myfile  << std::endl;
+                    mat_1.insert_vertex(u,1);
+                    // std::cout << "Inserted the vertex :" << renamedVertices.at(u) <<  " in the old skeleton matrix"<< std::endl;
+
+
+                }
+
+                // std::cout << "Going to insert the vertex :" << v << std::endl ;
+                if(!mat_1.membership(v)) {
+                    // std::cout << "Begining the insertion the vertex :" << renamedVertices.at(v) <<  " in the new distance matrix"<< std::endl;
+
+                    flag_Filtration->insert_vertex(renamedVertices.at(v),filtration_value);
+                    // std::cout << "Inserted the vertex :" << renamedVertices.at(v) <<  " in the new distance matrix"<< std::endl;
+
+                    // myfile << filtration_value << " i";
+                    // myfile  << " " << renamedVertices.at(v);
+                    // myfile  << std::endl;
+                    mat_1.insert_vertex(v,1);
+                    // std::cout << "Inserted the vertex :" << v <<  " in the old skeleton matrix"<< std::endl;
+
+                }
+                // std::cout << "Going to insert the edge :" << u << ", " << v << std::endl;
+                if(!mat_1.membership(e)){
+                    flag_Filtration->insert_new_edges(renamedVertices.at(u),renamedVertices.at(v), filtration_value);
+                    // std::cout << "Inserted the edge :" << renamedVertices.at(u) << ","<< renamedVertices.at(v) <<  " in the new distance matrix"<< std::endl;
+
+                    // myfile << filtration_value << " i";
+                    // myfile  << " " <<  renamedVertices.at(u) << ", " << renamedVertices.at(v);
+                    // myfile  << std::endl;
+                    // std::cout << "Going to insert the edge :" << u << ","<< v <<  " in the old skeleton"<< std::endl;
+                    mat_1.insert_new_edges(u,v,1);
+                    // std::cout << "Inserted the edge :" << u << ","<< v <<  " in the old skeleton"<< std::endl;
+                } 
+                 // std::cout << " Insertion Done " << std::endl;        
+            }
+
+            // myfile << "# Tower updated for the additional subcomplex.\n";
+            myfile.close();
+        }
+        else {
+            std::cerr << "Unable to open file";
+            exit(-1) ;
+        }
+        
+        auto end_print  = std::chrono::high_resolution_clock::now();
+        auto printTime = std::chrono::duration<double, std::milli>(end_print- begin_print).count();
+        // std::cout << " Time to print the tower : " << printTime << " ms\n" << std::endl;
+        return printTime;
+    } 
+    Distance_matrix distance_matrix()
+    {
+        size_t non_zero_rw = flag_Filtration->num_vertices();
+        double inf = std::numeric_limits<double>::max();
+        sparseRowMatrix mat = flag_Filtration->uncollapsed_matrix();
+        doubleVector distances ;  
+        for(size_t indx = 0; indx < non_zero_rw; indx++){ 
+            rowInnerIterator it(mat, indx);   
+            for(size_t j = 0; j <= indx; j++) {             // Iterate over the non-zero columns
+                //std::cout << "j = " << j << " - it.index = " << it.index() << " - indx = " << indx << std::endl;
+                if(it.index() == j && j != indx){
+                    distances.push_back(it.value());  // inner index, here it is equal to it.columns()
+                    ++it;
+                }
+                else if( j == indx)
+                    distances.push_back(0);
+                else
+                     distances.push_back(inf);
+
+            }
+            distance_mat->push_back(distances);
+            distances.clear();
+        }        
+        return *distance_mat;
+    }      
+   void print_sparse_matrix(){
+        flag_Filtration->print_sparse_skeleton();
+   }
+    
+};