summaryrefslogtreecommitdiff
path: root/src/Collapse
diff options
context:
space:
mode:
authorROUVREAU Vincent <vincent.rouvreau@inria.fr>2020-03-24 11:01:01 +0100
committerROUVREAU Vincent <vincent.rouvreau@inria.fr>2020-03-24 11:01:01 +0100
commit726329e4df9eb085da353ae25a949cbc7f66f00f (patch)
tree544dff85e571682dcfa32d0e73fbb99c906a602e /src/Collapse
parent1ad9a246611a33418720aa3c7753c858f55a2ae9 (diff)
clang format and remove useless variables and types
Diffstat (limited to 'src/Collapse')
-rw-r--r--src/Collapse/include/gudhi/FlagComplexSpMatrix.h1627
1 files changed, 785 insertions, 842 deletions
diff --git a/src/Collapse/include/gudhi/FlagComplexSpMatrix.h b/src/Collapse/include/gudhi/FlagComplexSpMatrix.h
index a1174084..ac02a46f 100644
--- a/src/Collapse/include/gudhi/FlagComplexSpMatrix.h
+++ b/src/Collapse/include/gudhi/FlagComplexSpMatrix.h
@@ -2,7 +2,7 @@
* (Geometric Understanding in Higher Dimensions) is a generic C++
* library for computational topology.
*
- * Author(s): Siddharth Pritam
+ * Author(s): Siddharth Pritam
*
* Copyright (C) 2018 INRIA Sophia Antipolis (France)
*
@@ -26,7 +26,6 @@
#include <boost/functional/hash.hpp>
// #include <boost/graph/adjacency_list.hpp>
-
#include <iostream>
#include <utility>
#include <vector>
@@ -42,868 +41,812 @@
#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 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 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 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 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 doubleQueue = std::queue<double>;
-using EdgeFiltQueue = std::queue<EdgeFilt>;
-using EdgeFiltVector = std::vector<EdgeFilt>;
+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;
+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
+//! Class SparseMsMatrix
/*!
- The class for storing the Vertices v/s MaxSimplices Sparse Matrix and performing collapses operations using the N^2() Algorithm.
+ 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;
+class FlagComplexSpMatrix {
+ private:
+ std::unordered_map<int, Vertex> rowToVertex;
- // Vertices strored as an unordered_set
- std::unordered_set<Vertex> vertices;
+ // 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.
- 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;
- 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)
-
- 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;
+ 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;
}
- sparseRowMatrix collapsed_matrix() const {
- return *sparse_colpsd_adj_Matrix;
+ 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;
+ }
+ }
+ }
}
-
- sparseRowMatrix uncollapsed_matrix() const {
- return sparseRowAdjMatrix;
+ // 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]);
+ }
+ }
}
-
- edge_list all_edges() const {
- return oneSimplices;
+ 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();
+ }
- 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;
+ 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++;
}
- 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;
+ 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))});
}
- 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;
+ // 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; }
-
- 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;
- }
+ 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