From ea53f69c13512465be3371f6378b5367f62164ab Mon Sep 17 00:00:00 2001 From: mcarrier Date: Wed, 20 Dec 2017 18:36:10 +0000 Subject: added extended persistence git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/Nerve_GIC@3092 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: a20ecf24407ce88654ad01257d8e44f8395c6b62 --- src/Nerve_GIC/include/gudhi/GIC.h | 168 +++++++++++++++++++++++--------------- 1 file changed, 102 insertions(+), 66 deletions(-) (limited to 'src/Nerve_GIC/include') diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h index 3cd8a92a..e72480b7 100644 --- a/src/Nerve_GIC/include/gudhi/GIC.h +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -31,7 +31,7 @@ #include #include #include -#include +//#include #include #include @@ -124,6 +124,8 @@ class Cover_complex { double rate_power = 0.001; // Power in the subsampling. int mask = 0; // Ignore nodes containing less than mask points. + std::map name2id, name2idinv; + std::string cover_name; std::string point_cloud_name; std::string color_name; @@ -600,66 +602,72 @@ class Cover_complex { std::vector points(n); for (int i = 0; i < n; i++) points[i] = i; std::sort(points.begin(), points.end(), Less(this->func)); - int id = 0; int pos = 0; int maxc = -1; IndexMap index = boost::get(boost::vertex_index, one_skeleton); + int id = 0; int pos = 0; IndexMap index = boost::get(boost::vertex_index, one_skeleton); //int maxc = -1; + std::map > preimages; std::map funcstd; - for (int i = 0; i < res; i++) { + if(verbose) std::cout << "Computing preimages..." << std::endl; + for (int i = 0; i < res; i++){ // Find points in the preimage - std::vector indices; std::pair inter1 = intervals[i]; - int tmp = pos; double u, v; Graph G = one_skeleton.create_subgraph(); + std::pair inter1 = intervals[i]; int tmp = pos; double u, v; if (i != res - 1) { if (i != 0) { std::pair inter3 = intervals[i - 1]; - while (func[points[tmp]] < inter3.second && tmp != n){ - boost::add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; - } + while (func[points[tmp]] < inter3.second && tmp != n){preimages[i].push_back(points[tmp]); tmp++;} u = inter3.second; } else u = inter1.first; std::pair inter2 = intervals[i + 1]; - while (func[points[tmp]] < inter2.first && tmp != n){ - boost::add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; - } - + while (func[points[tmp]] < inter2.first && tmp != n){preimages[i].push_back(points[tmp]); tmp++;} v = inter2.first; - pos = tmp; - while (func[points[tmp]] < inter1.second && tmp != n){ - boost::add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; - } + while (func[points[tmp]] < inter1.second && tmp != n){preimages[i].push_back(points[tmp]); tmp++;} } else { - std::pair inter3 = intervals[i - 1]; - while (func[points[tmp]] < inter3.second && tmp != n){ - boost::add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; - } - - while (tmp != n){ - boost::add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; - } + std::pair inter3 = intervals[i - 1]; + while (func[points[tmp]] < inter3.second && tmp != n){preimages[i].push_back(points[tmp]); tmp++;} + while (tmp != n){preimages[i].push_back(points[tmp]); tmp++;} u = inter3.second; v = inter1.second; } - int num = num_vertices(G); std::vector component(num); + funcstd[i] = 0.5*(u+v); + + } + + if(verbose) std::cout << "Computing connected components..." << std::endl; + for (int i = 0; i < res; i++){ // Compute connected components - boost::connected_components(G, &component[0]); int maxct = maxc + 1; + Graph G = one_skeleton.create_subgraph(); int num = preimages[i].size(); std::vector component(num); + for(int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G); + boost::connected_components(G, &component[0]); int max = 0; - // Update covers + // For each point in preimage for(int j = 0; j < num; j++){ - maxc = std::max(maxc, maxct + component[j]); - cover [indices[j]] .push_back(maxct + component[j]); - cover_back [maxct + component[j]] .push_back(indices[j]); - cover_fct [maxct + component[j]] = i; - cover_std [maxct + component[j]] = 0.5*(u+v); - cover_color [maxct + component[j]] .second += func_color[indices[j]]; - cover_color [maxct + component[j]] .first += 1; + + // Update number of components in preimage + if(component[j] > max) max = component[j]; + + // Identify component with Cantor polynomial N^2 -> N + int identifier = (std::pow(i + component[j], 2) + 3*i + component[j])/2; + + // Update covers + cover [preimages[i][j]] .push_back(identifier); + cover_back [identifier] .push_back(preimages[i][j]); + cover_fct [identifier] = i; + cover_std [identifier] = funcstd[i]; + cover_color [identifier] .second += func_color[preimages[i][j]]; + cover_color [identifier] .first += 1; } + + // Maximal dimension is total number of connected components + id += max + 1; + } maximal_dim = id - 1; @@ -747,7 +755,7 @@ class Cover_complex { * @result cover_back(c) vector of IDs of data points. * */ - const std::vector & subpopulation(int c) { return cover_back[c]; } + const std::vector & subpopulation(int c) { return cover_back[name2idinv[c]]; } // ******************************************************************************************************************* // Visualization. @@ -767,9 +775,7 @@ class Cover_complex { double f; while (std::getline(input, line)) { std::stringstream stream(line); - //stream >> one_skeleton[vertices[i]].color; - stream >> f; - func_color.emplace(i, f); + stream >> f; func_color.emplace(i, f); i++; } color_name = color_file_name; @@ -813,11 +819,11 @@ class Cover_complex { int k = 0; std::vector nodes; nodes.clear(); - graphic << "graph GIC {" << std::endl; + graphic << "graph GIC {" << std::endl; int id = 0; for (std::map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) { if (iit->second.first > mask) { - nodes.push_back(iit->first); - graphic << iit->first << "[shape=circle fontcolor=black color=black label=\"" << iit->first << ":" + nodes.push_back(iit->first); name2id[iit->first] = id; name2idinv[id] = iit->first; id++; + graphic << name2id[iit->first] << "[shape=circle fontcolor=black color=black label=\"" << name2id[iit->first] << ":" << iit->second.first << "\" style=filled fillcolor=\"" << (1 - (maxv - iit->second.second) / (maxv - minv)) * 0.6 << ", 1, 1\"]" << std::endl; k++; @@ -828,7 +834,7 @@ class Cover_complex { for (int i = 0; i < num_simplices; i++) if (simplices[i].size() == 2) { if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) { - graphic << " " << simplices[i][0] << " -- " << simplices[i][1] << " [weight=15];" << std::endl; + graphic << " " << name2id[simplices[i][0]] << " -- " << name2id[simplices[i][1]] << " [weight=15];" << std::endl; ke++; } } @@ -856,13 +862,14 @@ class Cover_complex { graphic << resolution_double << " " << gain << std::endl; graphic << cover_color.size() << " " << num_edges << std::endl; - for (std::map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) - graphic << iit->first << " " << iit->second.second << " " << iit->second.first << std::endl; + int id = 0; + for (std::map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++){ + graphic << id << " " << iit->second.second << " " << iit->second.first << std::endl; name2id[iit->first] = id; name2idinv[id] = iit->first; id++;} for (int i = 0; i < num_simplices; i++) if (simplices[i].size() == 2) if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) - graphic << simplices[i][0] << " " << simplices[i][1] << std::endl; + graphic << name2id[simplices[i][0]] << " " << name2id[simplices[i][1]] << std::endl; graphic.close(); std::cout << ".txt generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox." << std::endl; @@ -928,35 +935,64 @@ class Cover_complex { complex.insert_simplex_and_subfaces(simplex, filt); if (dimension < simplex.size() - 1) dimension = simplex.size() - 1; } - complex.set_dimension(dimension); + //complex.set_dimension(dimension); } public: /** \brief Computes the extended persistence diagram of the complex. * */ - template + //template void compute_PD(){ - SimplicialComplex streef, streeb; unsigned int dimension = 0; - for (auto const & simplex : simplices) { - int numvert = simplex.size(); double filtM = std::numeric_limits::lowest(); double filtm = filtM; - for(int i = 0; i < numvert; i++){filtM = std::max(cover_std[simplex[i]], filtM); filtm = std::max(-cover_std[simplex[i]], filtm);} - streef.insert_simplex_and_subfaces(simplex, filtM); streeb.insert_simplex_and_subfaces(simplex, filtm); - if (dimension < simplex.size() - 1) dimension = simplex.size() - 1; - } streef.set_dimension(dimension); streeb.set_dimension(dimension); + Simplex_tree sc; - streef.initialize_filtration(); - Gudhi::persistent_cohomology::Persistent_cohomology pcohf(streef); - pcohf.init_coefficients(2); pcohf.compute_persistent_cohomology(); - pcohf.output_diagram(); + // Add magic point + std::vector node = {-2}; sc.insert_simplex(node); + + // Compute max and min + double maxf = std::numeric_limits::lowest(); double minf = std::numeric_limits::max(); + for(std::map::iterator it = cover_std.begin(); it != cover_std.end(); it++){ + maxf = std::max(maxf, it->second); minf = std::min(minf, it->second); + } - streeb.initialize_filtration(); - Gudhi::persistent_cohomology::Persistent_cohomology pcohb(streeb); - pcohb.init_coefficients(2); pcohb.compute_persistent_cohomology(); - pcohb.output_diagram(); + // Build filtration + for (auto const & simplex : simplices) { + int numvert = simplex.size(); double filta = std::numeric_limits::lowest(); double filts = filta; + for(int i = 0; i < numvert; i++) filta = std::max( -2 + (cover_std[simplex[i]] - minf)/(maxf - minf), filta ); + for(int i = 0; i < numvert; i++) filts = std::max( 2 - (cover_std[simplex[i]] - minf)/(maxf - minf), filts ); + std::vector splx = simplex; splx.push_back(node[0]); + sc.insert_simplex_and_subfaces(simplex, filta); sc.insert_simplex_and_subfaces(splx, filts); + //for(int i = 0; i < numvert; i++) std::cout << simplex[i] << ", "; std::cout << std::endl << " with filter value " << filta << std::endl; + //for(int i = 0; i < numvert + 1; i++) std::cout << splx[i] << ", "; std::cout << std::endl << " with filter value " << filts << std::endl; + } - //PD = pcohf.get_persistent_pairs(); + // Compute PD + sc.initialize_filtration(); Gudhi::persistent_cohomology::Persistent_cohomology pcoh(sc); + pcoh.init_coefficients(2); pcoh.compute_persistent_cohomology(); + + // Output PD + std::vector > bars = pcoh.intervals_in_dimension(0); double birth_cc, death_cc; + int num_bars = bars.size(); std::cout << num_bars - 1 << " interval(s) in dimension 0:" << std::endl; + bool found_birth = false; bool found_death = false; + for(int j = 0; j < num_bars; j++){ + double birth = bars[j].first; double death = bars[j].second; + if(birth == 0){ std::cout << birth << " " << death << std::endl; if(found_birth) birth = birth_cc; else{death_cc = death; found_death = true; continue;} } + if(std::isinf(death)){ std::cout << birth << " " << death << std::endl; if(found_death) death = death_cc; else{birth_cc = birth; found_birth = true; continue;} } + if(birth < 0) birth = minf + (birth + 2)*(maxf - minf); else birth = minf + (2 - birth)*(maxf - minf); + if(death < 0) death = minf + (death + 2)*(maxf - minf); else death = minf + (2 - death)*(maxf - minf); + std::cout << " " << birth << ", " << death << std::endl; + } + + for(int i = 1; i < sc.dimension(); i++){ + bars = pcoh.intervals_in_dimension(i); num_bars = bars.size(); std::cout << num_bars << " interval(s) in dimension " << i << ":" << std::endl; + for(int j = 0; j < num_bars; j++){ + double birth = bars[j].first; double death = bars[j].second; + if(birth < 0) birth = minf + (birth + 2)*(maxf - minf); else birth = minf + (2 - birth)*(maxf - minf); + if(death < 0) death = minf + (death + 2)*(maxf - minf); else death = minf + (2 - death)*(maxf - minf); + std::cout << " " << birth << ", " << death << std::endl; + } + } } @@ -987,9 +1023,9 @@ class Cover_complex { Cboot.set_graph_from_automatic_rips(Gudhi::Euclidean_distance()); Cboot.set_automatic_resolution(); Cboot.set_gain(); Cboot.set_cover_from_function(); - Cboot.find_simplices(); Cboot.compute_PD >(); + Cboot.find_simplices(); Cboot.compute_PD(); - distribution.push_back(Gudhi::persistence_diagram::bottleneck_distance(this->PD,Cboot.PD)); + //distribution.push_back(Gudhi::persistence_diagram::bottleneck_distance(this->PD,Cboot.PD)); } @@ -1069,7 +1105,7 @@ class Cover_complex { for(int j = 0; j < numt; j++){ int vt = cover[index[boost::target(*ei, one_skeleton)]][j]; if(cover_fct[vs] == cover_fct[vt] + 1 || cover_fct[vt] == cover_fct[vs] + 1){ - std::vector edge(2); edge[0] = vs; edge[1] = vt; simplices.push_back(edge); goto afterLoop; + std::vector edge(2); edge[0] = std::min(vs,vt); edge[1] = std::max(vs,vt); simplices.push_back(edge); goto afterLoop; } } } -- cgit v1.2.3