diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/Doxyfile | 3 | ||||
| -rw-r--r-- | src/Nerve_GIC/example/CMakeLists.txt | 8 | ||||
| -rw-r--r-- | src/Nerve_GIC/example/simple_GIC.cpp | 65 | ||||
| -rw-r--r-- | src/Nerve_GIC/include/gudhi/GIC.h | 380 | ||||
| -rw-r--r-- | src/cmake/modules/GUDHI_user_version_target.txt | 2 |
5 files changed, 456 insertions, 2 deletions
diff --git a/src/Doxyfile b/src/Doxyfile index d2d0a447..60959235 100644 --- a/src/Doxyfile +++ b/src/Doxyfile @@ -851,7 +851,8 @@ IMAGE_PATH = doc/Skeleton_blocker/ \ doc/Subsampling/ \ doc/Spatial_searching/ \ doc/Tangential_complex/ \ - doc/Bottleneck_distance/ + doc/Bottleneck_distance/ \ + doc/Nerve_GIC/ # The INPUT_FILTER tag can be used to specify a program that doxygen should # invoke to filter for each input file. Doxygen will invoke the filter program diff --git a/src/Nerve_GIC/example/CMakeLists.txt b/src/Nerve_GIC/example/CMakeLists.txt new file mode 100644 index 00000000..e4debf2d --- /dev/null +++ b/src/Nerve_GIC/example/CMakeLists.txt @@ -0,0 +1,8 @@ +cmake_minimum_required(VERSION 2.6) +project(Nerve_GIC_examples) + +add_executable ( GIC simple_GIC.cpp ) +target_link_libraries(GIC ${Boost_SYSTEM_LIBRARY}) +if (TBB_FOUND) + target_link_libraries(GIC ${TBB_LIBRARIES}) +endif()
\ No newline at end of file diff --git a/src/Nerve_GIC/example/simple_GIC.cpp b/src/Nerve_GIC/example/simple_GIC.cpp new file mode 100644 index 00000000..e3d19cc8 --- /dev/null +++ b/src/Nerve_GIC/example/simple_GIC.cpp @@ -0,0 +1,65 @@ +#include <gudhi/GIC.h> + +void usage(int nbArgs, char * const progName) { + std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; + std::cerr << "Usage: " << progName << " filename.off threshold function resolution gain [ouput_file.txt]\n"; + std::cerr << " i.e.: " << progName << " ../../data/points/test.off 1.5 test_cov \n"; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if ((argc != 6) && (argc != 7)) usage(argc, (argv[0] - 1)); + + std::string off_file_name(argv[1]); + double threshold = atof(argv[2]); + std::string function_file_name(argv[3]); + double resolution = atof(argv[4]); + double gain = atof(argv[5]); + + // Type definitions + using Graph_t = boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS,\ + boost::property < vertex_filtration_t, Filtration_value >,\ + boost::property < edge_filtration_t, Filtration_value > >; + + // ---------------------------------------------------------------------------- + // Init of a graph induced complex from an OFF file + // ---------------------------------------------------------------------------- + + Gudhi::graph_induced_complex::Graph_induced_complex GIC; + GIC.set_graph_from_rips(threshold, off_file_name); + GIC.set_cover_from_function(function_file_name,resolution,gain,0); + //GIC.find_GIC_simplices(); + GIC.find_GIC_simplices_with_functional_minimal_cover(resolution,gain); + Simplex_tree stree; GIC.create_complex(stree); + + std::streambuf* streambufffer; + std::ofstream ouput_file_stream; + + if (argc == 7) { + ouput_file_stream.open(std::string(argv[4])); + streambufffer = ouput_file_stream.rdbuf(); + } else { + streambufffer = std::cout.rdbuf(); + } + + std::ostream output_stream(streambufffer); + + // ---------------------------------------------------------------------------- + // Display information about the graph induced complex + // ---------------------------------------------------------------------------- + output_stream << "Graph induced complex is of dimension " << stree.dimension() << + " - " << stree.num_simplices() << " simplices - " << + stree.num_vertices() << " vertices." << std::endl; + + output_stream << "Iterator on graph induced complex simplices" << std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + output_stream << vertex << " "; + } + output_stream << std::endl; + } + + ouput_file_stream.close(); + + return 0; +} diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h new file mode 100644 index 00000000..c36a36c9 --- /dev/null +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -0,0 +1,380 @@ +/* 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: Mathieu Carriere + * + * Copyright (C) 2017 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef GIC_H_ +#define GIC_H_ + +#include <gudhi/Debug_utils.h> +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/reader_utils.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Rips_complex.h> +#include <gudhi/Points_off_io.h> +#include <gudhi/distance_functions.h> + +#include <boost/graph/adjacency_list.hpp> + +#include <iostream> +#include <vector> +#include <map> +#include <string> +#include <limits> // for numeric_limits +#include <utility> // for pair<> +#include <functional> // for greater<> +#include <stdexcept> +#include <initializer_list> +#include <algorithm> // for std::max +#include <cstdint> // for std::uint32_t + +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Point = std::vector<float>; + +std::map<int, double> func; + +namespace Gudhi { + +namespace graph_induced_complex { + + +/** + * \class Graph_induced_complex + * \brief Graph induced complex data structure. + * + * \ingroup graph_induced_complex + * + * \details + * + * + */ + +class Graph_induced_complex { + + private: + typedef int Cover_t; + + private: + std::vector<std::vector<Cover_t> > simplices; + + private: + std::map<int, std::vector<Cover_t> > cover; + + private: + int maximal_dim; + + private: + Simplex_tree<> st; + std::map<int,std::vector<int> > adjacency_matrix; + + + // Simplex comparator + private: + bool simplex_comp(const std::vector<Cover_t>& s1, const std::vector<Cover_t>& s2){ + if(s1.size() == s2.size()){ + return s1[0] < s2[0]; + } + else return s1.size() < s2.size(); + } + + // Point comparator + private: + static bool functional_comp(const int& a, const int& b){ + if(func[a] == func[b]) return a < b; + else return func[a] < func[b]; + } + + private: + void dfs(std::map<int,std::vector<int> >& G, const int& p, std::vector<int>& cc, std::map<int,bool>& visit){ + cc.push_back(p); + visit[p] = true; int neighb = G[p].size(); + for (int i = 0; i < neighb; i++) + if ( visit.find(G[p][i]) != visit.end() ) + if( !(visit[G[p][i]]) ) + dfs(G,G[p][i],cc,visit); + } + + public: + void set_cover_from_file(const std::string& cover_file_name){ + int vertex_id = 0; Cover_t cov; std::vector<Cover_t> cov_elts, cov_number; + std::ifstream input(cover_file_name); std::string line; + while(std::getline(input,line)){ + cov_elts.clear(); std::stringstream stream(line); + while(stream >> cov){cov_elts.push_back(cov); cov_number.push_back(cov);} + cover.insert(std::pair<int, std::vector<Cover_t> >(vertex_id, cov_elts)); vertex_id++; + } + std::vector<Cover_t>::iterator it; + std::sort(cov_number.begin(),cov_number.end()); it = std::unique(cov_number.begin(),cov_number.end()); + cov_number.resize(std::distance(cov_number.begin(),it)); maximal_dim = cov_number.size(); + } + + public: + void set_cover_from_function(const std::string& func_file_name, const double& resolution, const double& gain, const bool& token){ + + // Read function values and compute min and max + int vertex_id = 0; double f; std::vector<double> range; double maxf, minf; minf = std::numeric_limits<float>::max(); maxf = std::numeric_limits<float>::min(); + std::ifstream input(func_file_name); std::string line; + while(std::getline(input,line)){ + std::stringstream stream(line); + stream >> f; range.push_back(f); minf = std::min(minf, f); maxf = std::max(maxf, f); + func.insert(std::pair<int,double>(vertex_id, f)); vertex_id++; + } + int num_pts = func.size(); + + // Compute cover of im(f) + std::vector<std::pair<double,double> > intervals; int res; + if(!token){ + double incr = (maxf-minf)/resolution; double x = minf; double alpha = (incr*gain)/(2-2*gain); + double y = minf + incr + alpha; std::pair<double, double> interm(x,y); intervals.push_back(interm); + for(int i = 1; i < resolution-1; i++){ + x = minf + i*incr - alpha; + y = minf + (i+1)*incr + alpha; + std::pair<double, double> inter(x,y); intervals.push_back(inter); + } + x = minf + (resolution-1)*incr - alpha; y = maxf; + std::pair<double, double> interM(x,y); intervals.push_back(interM); res = intervals.size(); + } + else{ + double x = minf; double y = x + resolution; + while(y <= maxf && maxf - (y-gain*resolution) >= resolution){ + std::pair<double, double> inter(x,y); intervals.push_back(inter); + x = y - gain*resolution; + y = x + resolution; + } + std::pair<double, double> interM(x,maxf); intervals.push_back(interM); res = intervals.size(); + } + + // Sort points according to function values + std::vector<int> points(num_pts); for(int i = 0; i < num_pts; i++) points[i] = i; + std::sort(points.begin(),points.end(),functional_comp); + + // Build adjacency matrix + std::vector<int> empty; + for(int i = 0; i < num_pts; i++) adjacency_matrix.insert(std::pair<int,std::vector<int> >(points[i],empty)); + for (auto simplex : st.complex_simplex_range()) { + if(st.dimension(simplex) == 1){ + std::vector<int> vertices; + for(auto vertex : st.simplex_vertex_range(simplex)) vertices.push_back(vertex); + adjacency_matrix[vertices[0]].push_back(vertices[1]); adjacency_matrix[vertices[1]].push_back(vertices[0]); + } + } + + int id = 0; int pos = 0; + for(int i = 0; i < res; i++){ + + // Find points in the preimage + std::map<int,std::vector<int> > prop; prop.clear(); + std::pair<double, double> inter1 = intervals[i]; + int tmp = pos; + + if(i != res-1){ + if(i != 0){ + std::pair<double, double> inter3 = intervals[i-1]; + while(func[points[tmp]] < inter3.second && tmp != num_pts){ + prop.insert(std::pair<int,std::vector<int> >(points[tmp],adjacency_matrix[points[tmp]])); + tmp++; + } + } + std::pair<double, double> inter2 = intervals[i+1]; + while(func[points[tmp]] < inter2.first && tmp != num_pts){ + prop.insert(std::pair<int,std::vector<int> >(points[tmp],adjacency_matrix[points[tmp]])); + tmp++; + } + pos = tmp; + while(func[points[tmp]] < inter1.second && tmp != num_pts){ + prop.insert(std::pair<int,std::vector<int> >(points[tmp],adjacency_matrix[points[tmp]])); + tmp++; + } + + } + else{ + std::pair<double, double> inter3 = intervals[i-1]; + while(func[points[tmp]] < inter3.second && tmp != num_pts){ + prop.insert(std::pair<int,std::vector<int> >(points[tmp],adjacency_matrix[points[tmp]])); + tmp++; + } + while(tmp != num_pts){ + prop.insert(std::pair<int,std::vector<int> >(points[tmp],adjacency_matrix[points[tmp]])); + tmp++; + } + + } + + // Compute the connected components with DFS + std::map<int,bool> visit; + for(std::map<int, std::vector<int> >::iterator it = prop.begin(); it != prop.end(); it++) + visit.insert(std::pair<int,bool>(it->first, false)); + if (!(prop.empty())){ + for(std::map<int, std::vector<int> >::iterator it = prop.begin(); it != prop.end(); it++){ + if ( !(visit[it->first]) ){ + std::vector<int> cc; cc.clear(); + dfs(prop,it->first,cc,visit); int cci = cc.size(); + for(int i = 0; i < cci; i++) cover[cc[i]].push_back(id); + id++; + } + } + } + + } + + maximal_dim = id; + + } + + public: + void set_graph_from_rips(const double& threshold, const std::string& off_file_name){ + Points_off_reader<Point> off_reader(off_file_name); + Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), threshold, Euclidean_distance()); + rips_complex_from_points.create_complex(st, 1); + } + + public: + void set_graph_from_file(const std::string& graph_file_name){ + int neighb; int vid; std::ifstream input(graph_file_name); std::string line; std::vector<int> edge(2); + while(std::getline(input,line)){ + std::stringstream stream(line); stream >> vid; edge[0] = vid; + while(stream >> neighb){edge[1] = neighb; st.insert_simplex_and_subfaces(edge);} + } + } + + public: + template<typename SimplicialComplexForGIC> + void create_complex(SimplicialComplexForGIC & complex) { + size_t sz = simplices.size(); int dimension = 0; + for(int i = 0; i < sz; i++){ + complex.insert_simplex_and_subfaces(simplices[i]); + if(dimension < simplices[i].size()-1) dimension = simplices[i].size()-1; + } + complex.set_dimension(dimension); + } + + public: + void find_all_simplices(const std::vector<std::vector<Cover_t> > & cover_elts,\ + int token, std::vector<Cover_t> & simplex_tmp){ + int num_nodes = cover_elts.size(); + if(token == num_nodes-1){ + int num_clus = cover_elts[token].size(); + for(int i = 0; i < num_clus; i++){ + std::vector<Cover_t> simplex = simplex_tmp; simplex.push_back(cover_elts[token][i]); + std::vector<Cover_t>::iterator it; + std::sort(simplex.begin(),simplex.end()); it = std::unique(simplex.begin(),simplex.end()); + simplex.resize(std::distance(simplex.begin(),it)); + simplices.push_back(simplex); + } + } + else{ + int num_clus = cover_elts[token].size(); + for(int i = 0; i < num_clus; i++){ + std::vector<Cover_t> simplex = simplex_tmp; simplex.push_back(cover_elts[token][i]); + find_all_simplices(cover_elts, token+1, simplex); + } + } + } + + public: + void find_Nerve_simplices(){ + simplices.clear(); + for(std::map<int,std::vector<Cover_t> >::iterator it = cover.begin(); it!=cover.end(); it++){ + simplices.push_back(it->second); + } + std::vector<std::vector<Cover_t> >::iterator it; + std::sort(simplices.begin(),simplices.end()); it = std::unique(simplices.begin(),simplices.end()); + simplices.resize(std::distance(simplices.begin(),it)); + } + + public: + void find_GIC_simplices() { + + // Find IDs of edges to remove + std::vector<int> simplex_to_remove; int simplex_id = 0; + for (auto simplex : st.complex_simplex_range()) { + if(st.dimension(simplex) == 1){ + std::vector<std::vector<Cover_t> > comp; + for(auto vertex : st.simplex_vertex_range(simplex)) comp.push_back(cover[vertex]); + if(comp[0].size() == 1 && comp[1].size() == 1 && comp[0] == comp[1]) simplex_to_remove.push_back(simplex_id); + } + simplex_id++; + } + + // Remove edges + int current_id = 1; + auto simplex = st.complex_simplex_range().begin(); int num_rem = 0; + for(int i = 0; i < simplex_id-1; i++){ + int j = i+1; auto simplex_tmp = simplex; simplex_tmp++; + if(j == simplex_to_remove[current_id]){st.remove_maximal_simplex(*simplex_tmp); current_id++; num_rem++;} + else simplex++; + } simplex = st.complex_simplex_range().begin(); + for(int i = 0; i < simplex_to_remove[0]; i++) simplex++; st.remove_maximal_simplex(*simplex); + + // Build the Simplex Tree corresponding to the graph + st.expansion(maximal_dim); + + // Find simplices of GIC + simplices.clear(); + for (auto simplex : st.complex_simplex_range()) { + if(!st.has_children(simplex)){ + std::vector<std::vector<Cover_t> > cover_elts; std::vector<Cover_t> sim; + for (auto vertex : st.simplex_vertex_range(simplex)) cover_elts.push_back(cover[vertex]); + find_all_simplices(cover_elts,0,sim); + } + } + std::vector<std::vector<Cover_t> >::iterator it; + std::sort(simplices.begin(),simplices.end()); it = std::unique(simplices.begin(),simplices.end()); + simplices.resize(std::distance(simplices.begin(),it)); + } + + public: + void find_GIC_simplices_with_functional_minimal_cover(const double& resolution, const double& gain){ + for(std::map<int,std::vector<Cover_t> >::iterator it = cover.begin(); it != cover.end(); it++){ + int vid = it->first; std::vector<int> neighbors = adjacency_matrix[vid]; int num_neighb = neighbors.size(); + for(int i = 0; i < num_neighb; i++){ + int neighb = neighbors[i]; int v1, v2; + if(func[vid] > func[neighb]){ + if( func[vid]-func[neighb] <= resolution*(2-gain) ){ + if(cover[vid].size() == 2) v1 = std::min(cover[vid][0],cover[vid][1]); else v1 = cover[vid][0]; + if(cover[neighb].size() == 2) v2 = std::max(cover[neighb][0],cover[neighb][1]); else v2 = cover[neighb][0]; + std::vector<int> edge(2); edge[0] = v1; edge[1] = v2; + if(v1 > v2) simplices.push_back(edge); + } + } + else{ + if( func[neighb]-func[vid] <= resolution*(2-gain) ){ + if(cover[vid].size() == 2) v1 = std::max(cover[vid][0],cover[vid][1]); else v1 = cover[vid][0]; + if(cover[neighb].size() == 2) v2 = std::min(cover[neighb][0],cover[neighb][1]); else v2 = cover[neighb][0]; + std::vector<int> edge(2); edge[0] = v1; edge[1] = v2; + if(v2 > v1) simplices.push_back(edge); + } + } + } + } + std::vector<std::vector<Cover_t> >::iterator it; + std::sort(simplices.begin(),simplices.end()); it = std::unique(simplices.begin(),simplices.end()); + simplices.resize(std::distance(simplices.begin(),it)); + } + +}; + +} // namespace graph_induced_complex + +} // namespace Gudhi + +#endif // GIC_H_ diff --git a/src/cmake/modules/GUDHI_user_version_target.txt b/src/cmake/modules/GUDHI_user_version_target.txt index 13fccd32..e4c4ca93 100644 --- a/src/cmake/modules/GUDHI_user_version_target.txt +++ b/src/cmake/modules/GUDHI_user_version_target.txt @@ -49,7 +49,7 @@ if (NOT CMAKE_VERSION VERSION_LESS 2.8.11) add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_SOURCE_DIR}/src/GudhUI ${GUDHI_USER_VERSION_DIR}/GudhUI) - set(GUDHI_MODULES "common;Alpha_complex;Bitmap_cubical_complex;Bottleneck_distance;Contraction;Hasse_complex;Persistent_cohomology;Rips_complex;Simplex_tree;Skeleton_blocker;Spatial_searching;Subsampling;Tangential_complex;Witness_complex") + set(GUDHI_MODULES "common;Alpha_complex;Bitmap_cubical_complex;Bottleneck_distance;Contraction;Hasse_complex;Nerve_GIC;Persistent_cohomology;Rips_complex;Simplex_tree;Skeleton_blocker;Spatial_searching;Subsampling;Tangential_complex;Witness_complex") set(GUDHI_DIRECTORIES "doc;example;concept") set(GUDHI_INCLUDE_DIRECTORIES "include/gudhi;include/gudhi_patches") |