diff options
Diffstat (limited to 'src/Nerve_GIC')
29 files changed, 2722 insertions, 0 deletions
diff --git a/src/Nerve_GIC/doc/COPYRIGHT b/src/Nerve_GIC/doc/COPYRIGHT new file mode 100644 index 00000000..61f17f6d --- /dev/null +++ b/src/Nerve_GIC/doc/COPYRIGHT @@ -0,0 +1,12 @@ +The files of this directory are part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. +See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + +Author(s): Vincent Rouvreau + +Copyright (C) 2015 Inria + +This gives everyone the freedoms to use openFrameworks in any context: +commercial or non-commercial, public or private, open or closed source. + +You should have received a copy of the MIT License along with this program. +If not, see https://opensource.org/licenses/MIT.
\ No newline at end of file diff --git a/src/Nerve_GIC/doc/GIC.jpg b/src/Nerve_GIC/doc/GIC.jpg Binary files differnew file mode 100644 index 00000000..cb1b9b7f --- /dev/null +++ b/src/Nerve_GIC/doc/GIC.jpg diff --git a/src/Nerve_GIC/doc/GIC.pdf b/src/Nerve_GIC/doc/GIC.pdf Binary files differnew file mode 100644 index 00000000..30525745 --- /dev/null +++ b/src/Nerve_GIC/doc/GIC.pdf diff --git a/src/Nerve_GIC/doc/Intro_graph_induced_complex.h b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h new file mode 100644 index 00000000..f9441b24 --- /dev/null +++ b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h @@ -0,0 +1,173 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Mathieu Carriere + * + * Copyright (C) 2017 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef DOC_COVER_COMPLEX_INTRO_COVER_COMPLEX_H_ +#define DOC_COVER_COMPLEX_INTRO_COVER_COMPLEX_H_ + +namespace Gudhi { + +namespace cover_complex { + +/** \defgroup cover_complex Cover complex + * + * \author Mathieu Carrière + * + * @{ + * + * Visualizations of the simplicial complexes can be done with either + * neato (from <a target="_blank" href="http://www.graphviz.org/">graphviz</a>), + * <a target="_blank" href="http://www.geomview.org/">geomview</a>, + * <a target="_blank" href="https://github.com/MLWave/kepler-mapper">KeplerMapper</a>. + * Input point clouds are assumed to be \ref FileFormatsOFF "OFF files" + * + * \section covers Covers + * + * Nerves and Graph Induced Complexes require a cover C of the input point cloud P, + * that is a set of subsets of P whose union is P itself. + * Very often, this cover is obtained from the preimage of a family of intervals covering + * the image of some scalar-valued function f defined on P. This family is parameterized + * by its resolution, which can be either the number or the length of the intervals, + * and its gain, which is the overlap percentage between consecutive intervals (ordered by their first values). + * + * \section nerves Nerves + * + * \subsection nervedefinition Nerve definition + * + * Assume you are given a cover C of your point cloud P. Then, the Nerve of this cover + * is the simplicial complex that has one k-simplex per k-fold intersection of cover elements. + * See also <a target="_blank" href="https://en.wikipedia.org/wiki/Nerve_of_a_covering"> Wikipedia </a>. + * + * \image html "nerve.png" "Nerve of a double torus" + * + * \subsection nerveexample Example + * + * This example builds the Nerve of a point cloud sampled on a 3D human shape (human.off). + * The cover C comes from the preimages of intervals (10 intervals with gain 0.3) + * covering the height function (coordinate 2), + * which are then refined into their connected components using the triangulation of the .OFF file. + * + * \include Nerve_GIC/Nerve.cpp + * + * When launching: + * + * \code $> ./Nerve ../../data/points/human.off 2 10 0.3 -v + * \endcode + * + * the program output is: + * + * \include Nerve_GIC/Nerve.txt + * + * The program also writes a file ../../data/points/human_sc.txt. The first three lines in this file are the location + * of the input point cloud and the function used to compute the cover. + * The fourth line contains the number of vertices nv and edges ne of the Nerve. + * The next nv lines represent the vertices. Each line contains the vertex ID, + * the number of data points it contains, and their average color function value. + * Finally, the next ne lines represent the edges, characterized by the ID of their vertices. + * + * Using KeplerMapper, one can obtain the following visualization: + * + * \image html "nervevisu.jpg" "Visualization with KeplerMapper" + * + * \section gic Graph Induced Complexes (GIC) + * + * \subsection gicdefinition GIC definition + * + * Again, assume you are given a cover C of your point cloud P. Moreover, assume + * you are also given a graph G built on top of P. Then, for any clique in G + * whose nodes all belong to different elements of C, the GIC includes a corresponding + * simplex, whose dimension is the number of nodes in the clique minus one. + * See \cite Dey13 for more details. + * + * \image html "GIC.jpg" "GIC of a point cloud." + * + * \subsection gicexamplevor Example with cover from Voronoï + * + * This example builds the GIC of a point cloud sampled on a 3D human shape (human.off). + * We randomly subsampled 100 points in the point cloud, which act as seeds of + * a geodesic Voronoï diagram. Each cell of the diagram is then an element of C. + * The graph G (used to compute both the geodesics for Voronoï and the GIC) + * comes from the triangulation of the human shape. Note that the resulting simplicial complex is in dimension 3 + * in this example. + * + * \include Nerve_GIC/VoronoiGIC.cpp + * + * When launching: + * + * \code $> ./VoronoiGIC ../../data/points/human.off 700 -v + * \endcode + * + * the program outputs SC.off. Using e.g. + * + * \code $> geomview ../../data/points/human_sc.off + * \endcode + * + * one can obtain the following visualization: + * + * \image html "gicvoronoivisu.jpg" "Visualization with Geomview" + * + * \subsection functionalGICdefinition Functional GIC + * + * If one restricts to the cliques in G whose nodes all belong to preimages of consecutive + * intervals (assuming the cover of the height function is minimal, i.e. no more than + * two intervals can intersect at a time), the GIC is of dimension one, i.e. a graph. + * We call this graph the functional GIC. See \cite Carriere16 for more details. + * + * \subsection functionalGICexample Example + * + * Functional GIC comes with automatic selection of the Rips threshold, + * the resolution and the gain of the function cover. See \cite Carriere17c for more details. In this example, + * we compute the functional GIC of a Klein bottle embedded in R^5, + * where the graph G comes from a Rips complex with automatic threshold, + * and the cover C comes from the preimages of intervals covering the first coordinate, + * with automatic resolution and gain. Note that automatic threshold, resolution and gain + * can be computed as well for the Nerve. + * + * \include Nerve_GIC/CoordGIC.cpp + * + * When launching: + * + * \code $> ./CoordGIC ../../data/points/KleinBottle5D.off 0 -v + * \endcode + * + * the program outputs SC.dot. Using e.g. + * + * \code $> neato SC.dot -Tpdf -o SC.pdf + * \endcode + * + * one can obtain the following visualization: + * + * \image html "coordGICvisu2.jpg" "Visualization with Neato" + * + * where nodes are colored by the filter function values and, for each node, the first number is its ID + * and the second is the number of data points that its contain. + * + * We also provide an example on a set of 72 pictures taken around the same object (lucky_cat.off). + * The function is now the first eigenfunction given by PCA, whose values + * are written in a file (lucky_cat_PCA1). Threshold, resolution and gain are automatically selected as before. + * + * \include Nerve_GIC/FuncGIC.cpp + * + * When launching: + * + * \code $> ./FuncGIC ../../data/points/COIL_database/lucky_cat.off ../../data/points/COIL_database/lucky_cat_PCA1 -v + * \endcode + * + * the program outputs again SC.dot which gives the following visualization after using neato: + * + * \image html "funcGICvisu.jpg" "Visualization with neato" + * + */ +/** @} */ // end defgroup cover_complex + +} // namespace cover_complex + +} // namespace Gudhi + +#endif // DOC_COVER_COMPLEX_INTRO_COVER_COMPLEX_H_ diff --git a/src/Nerve_GIC/doc/coordGICvisu.pdf b/src/Nerve_GIC/doc/coordGICvisu.pdf Binary files differnew file mode 100644 index 00000000..313aa1b5 --- /dev/null +++ b/src/Nerve_GIC/doc/coordGICvisu.pdf diff --git a/src/Nerve_GIC/doc/coordGICvisu2.jpg b/src/Nerve_GIC/doc/coordGICvisu2.jpg Binary files differnew file mode 100644 index 00000000..046feb2a --- /dev/null +++ b/src/Nerve_GIC/doc/coordGICvisu2.jpg diff --git a/src/Nerve_GIC/doc/funcGICvisu.jpg b/src/Nerve_GIC/doc/funcGICvisu.jpg Binary files differnew file mode 100644 index 00000000..36b64dde --- /dev/null +++ b/src/Nerve_GIC/doc/funcGICvisu.jpg diff --git a/src/Nerve_GIC/doc/funcGICvisu.pdf b/src/Nerve_GIC/doc/funcGICvisu.pdf Binary files differnew file mode 100644 index 00000000..d7456cd3 --- /dev/null +++ b/src/Nerve_GIC/doc/funcGICvisu.pdf diff --git a/src/Nerve_GIC/doc/gicvisu.jpg b/src/Nerve_GIC/doc/gicvisu.jpg Binary files differnew file mode 100644 index 00000000..576dae47 --- /dev/null +++ b/src/Nerve_GIC/doc/gicvisu.jpg diff --git a/src/Nerve_GIC/doc/gicvoronoivisu.jpg b/src/Nerve_GIC/doc/gicvoronoivisu.jpg Binary files differnew file mode 100644 index 00000000..cd86c411 --- /dev/null +++ b/src/Nerve_GIC/doc/gicvoronoivisu.jpg diff --git a/src/Nerve_GIC/doc/nerve.png b/src/Nerve_GIC/doc/nerve.png Binary files differnew file mode 100644 index 00000000..b66da4a4 --- /dev/null +++ b/src/Nerve_GIC/doc/nerve.png diff --git a/src/Nerve_GIC/doc/nervevisu.jpg b/src/Nerve_GIC/doc/nervevisu.jpg Binary files differnew file mode 100644 index 00000000..67ae1d7e --- /dev/null +++ b/src/Nerve_GIC/doc/nervevisu.jpg diff --git a/src/Nerve_GIC/example/CMakeLists.txt b/src/Nerve_GIC/example/CMakeLists.txt new file mode 100644 index 00000000..1667472f --- /dev/null +++ b/src/Nerve_GIC/example/CMakeLists.txt @@ -0,0 +1,28 @@ +project(Nerve_GIC_examples) + +if (NOT CGAL_VERSION VERSION_LESS 4.11.0) + + add_executable ( CoordGIC CoordGIC.cpp ) + add_executable ( FuncGIC FuncGIC.cpp ) + + if (TBB_FOUND) + target_link_libraries(CoordGIC ${TBB_LIBRARIES}) + target_link_libraries(FuncGIC ${TBB_LIBRARIES}) + endif() + + # Copy files for not to pollute sources when testing + file(COPY "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + file(COPY "${CMAKE_SOURCE_DIR}/data/points/COIL_database/lucky_cat.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + file(COPY "${CMAKE_SOURCE_DIR}/data/points/COIL_database/lucky_cat_PCA1" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + + add_test(NAME Nerve_GIC_example_CoordGIC COMMAND $<TARGET_FILE:CoordGIC> + "${CMAKE_CURRENT_BINARY_DIR}/tore3D_1307.off" "0") + + add_test(NAME Nerve_GIC_example_FuncGIC COMMAND $<TARGET_FILE:FuncGIC> + "${CMAKE_CURRENT_BINARY_DIR}/lucky_cat.off" + "${CMAKE_CURRENT_BINARY_DIR}/lucky_cat_PCA1") + + install(TARGETS CoordGIC DESTINATION bin) + install(TARGETS FuncGIC DESTINATION bin) + +endif (NOT CGAL_VERSION VERSION_LESS 4.11.0) diff --git a/src/Nerve_GIC/example/CoordGIC.cpp b/src/Nerve_GIC/example/CoordGIC.cpp new file mode 100644 index 00000000..fd9c224a --- /dev/null +++ b/src/Nerve_GIC/example/CoordGIC.cpp @@ -0,0 +1,84 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Mathieu Carrière + * + * Copyright (C) 2017 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/GIC.h> + +#include <string> +#include <vector> + +void usage(int nbArgs, char *const progName) { + std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; + std::cerr << "Usage: " << progName << " filename.off coordinate [-v] \n"; + std::cerr << " i.e.: " << progName << " ../../data/points/human.off 2 -v \n"; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if ((argc != 3) && (argc != 4)) usage(argc, argv[0]); + + using Point = std::vector<float>; + + std::string off_file_name(argv[1]); + int coord = atoi(argv[2]); + bool verb = 0; + if (argc == 4) verb = 1; + + // ----------------------------------------- + // Init of a functional GIC from an OFF file + // ----------------------------------------- + + Gudhi::cover_complex::Cover_complex<Point> GIC; + GIC.set_verbose(verb); + + bool check = GIC.read_point_cloud(off_file_name); + + if (!check) { + std::cout << "Incorrect OFF file." << std::endl; + } else { + GIC.set_type("GIC"); + + GIC.set_color_from_coordinate(coord); + GIC.set_function_from_coordinate(coord); + + GIC.set_graph_from_automatic_rips(Gudhi::Euclidean_distance()); + GIC.set_automatic_resolution(); + GIC.set_gain(); + GIC.set_cover_from_function(); + + GIC.find_simplices(); + + GIC.compute_distribution(10); + GIC.compute_p_value(); + + GIC.plot_DOT(); + + Gudhi::Simplex_tree<> stree; + GIC.create_complex(stree); + + // -------------------------------------------- + // Display information about the functional GIC + // -------------------------------------------- + + if (verb) { + std::cout << "Coordinate GIC is of dimension " << stree.dimension() << " - " << stree.num_simplices() + << " simplices - " << stree.num_vertices() << " vertices." << std::endl; + + std::cout << "Iterator on coordinate GIC simplices" << std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << std::endl; + } + } + } + + return 0; +} diff --git a/src/Nerve_GIC/example/FuncGIC.cpp b/src/Nerve_GIC/example/FuncGIC.cpp new file mode 100644 index 00000000..5a323795 --- /dev/null +++ b/src/Nerve_GIC/example/FuncGIC.cpp @@ -0,0 +1,82 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Mathieu Carrière + * + * Copyright (C) 2017 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/GIC.h> + +#include <string> +#include <vector> + +void usage(int nbArgs, char *const progName) { + std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; + std::cerr << "Usage: " << progName << " filename.off function [-v] \n"; + std::cerr << " i.e.: " << progName << " ../../data/points/COIL_database/lucky_cat.off " + "../../data/points/COIL_database/lucky_cat_PCA1 -v \n"; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if ((argc != 3) && (argc != 4)) usage(argc, argv[0]); + + using Point = std::vector<float>; + + std::string off_file_name(argv[1]); + std::string func_file_name = argv[2]; + bool verb = 0; + if (argc == 4) verb = 1; + + // ----------------------------------------- + // Init of a functional GIC from an OFF file + // ----------------------------------------- + + Gudhi::cover_complex::Cover_complex<Point> GIC; + GIC.set_verbose(verb); + + bool check = GIC.read_point_cloud(off_file_name); + + if (!check) { + std::cout << "Incorrect OFF file." << std::endl; + } else { + GIC.set_type("GIC"); + + GIC.set_color_from_file(func_file_name); + GIC.set_function_from_file(func_file_name); + + GIC.set_graph_from_automatic_rips(Gudhi::Euclidean_distance()); + GIC.set_automatic_resolution(); + GIC.set_gain(); + GIC.set_cover_from_function(); + + GIC.find_simplices(); + + GIC.plot_DOT(); + + Gudhi::Simplex_tree<> stree; + GIC.create_complex(stree); + + // -------------------------------------------- + // Display information about the functional GIC + // -------------------------------------------- + + if (verb) { + std::cout << "Functional GIC is of dimension " << stree.dimension() << " - " << stree.num_simplices() + << " simplices - " << stree.num_vertices() << " vertices." << std::endl; + + std::cout << "Iterator on functional GIC simplices" << std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << std::endl; + } + } + } + + 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..fc6a2a91 --- /dev/null +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -0,0 +1,1420 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author: Mathieu Carriere + * + * Copyright (C) 2017 Inria + * + * Modification(s): + * - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL + * - YYYY/MM Author: Description of the modification + */ + +#ifndef GIC_H_ +#define GIC_H_ + +#ifdef GUDHI_USE_TBB +#include <tbb/parallel_for.h> +#include <tbb/mutex.h> +#endif + +#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 <gudhi/Persistent_cohomology.h> +#include <gudhi/Bottleneck.h> + +#include <boost/config.hpp> +#include <boost/graph/graph_traits.hpp> +#include <boost/graph/adjacency_list.hpp> +#include <boost/graph/connected_components.hpp> +#include <boost/graph/dijkstra_shortest_paths.hpp> +#include <boost/graph/subgraph.hpp> +#include <boost/graph/graph_utility.hpp> + +#include <CGAL/version.h> // for CGAL_VERSION_NR + +#include <iostream> +#include <vector> +#include <map> +#include <string> +#include <limits> // for numeric_limits +#include <utility> // for std::pair<> +#include <algorithm> // for (std::max) +#include <random> +#include <cassert> +#include <cmath> + +// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10 +#if CGAL_VERSION_NR < 1041101000 +# error Alpha_complex_3d is only available for CGAL >= 4.11 +#endif + +namespace Gudhi { + +namespace cover_complex { + +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Persistence_diagram = std::vector<std::pair<double, double> >; +using Graph = boost::subgraph< + boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS, boost::no_property, + boost::property<boost::edge_index_t, int, boost::property<boost::edge_weight_t, double> > > >; +using Vertex_t = boost::graph_traits<Graph>::vertex_descriptor; +using Index_map = boost::property_map<Graph, boost::vertex_index_t>::type; +using Weight_map = boost::property_map<Graph, boost::edge_weight_t>::type; + +/** + * \class Cover_complex + * \brief Cover complex data structure. + * + * \ingroup cover_complex + * + * \details + * The data structure is a simplicial complex, representing a + * Graph Induced simplicial Complex (GIC) or a Nerve, + * and whose simplices are computed with a cover C of a point + * cloud P, which often comes from the preimages of intervals + * covering the image of a function f defined on P. + * These intervals are parameterized by their resolution + * (either their length or their number) + * and their gain (percentage of overlap). + * To compute a GIC, one also needs a graph G built on top of P, + * whose cliques with vertices belonging to different elements of C + * correspond to the simplices of the GIC. + * + */ +template <typename Point> +class Cover_complex { + private: + bool verbose = false; // whether to display information. + std::string type; // Nerve or GIC + + std::vector<Point> point_cloud; // input point cloud. + std::vector<std::vector<double> > distances; // all pairwise distances. + int maximal_dim; // maximal dimension of output simplicial complex. + int data_dimension; // dimension of input data. + int n; // number of points. + + std::vector<double> func; // function used to compute the output simplicial complex. + std::vector<double> func_color; // function used to compute the colors of the nodes of the output simplicial complex. + bool functional_cover = false; // whether we use a cover with preimages of a function or not. + + Graph one_skeleton_OFF; // one-skeleton given by the input OFF file (if it exists). + Graph one_skeleton; // one-skeleton used to compute the connected components. + std::vector<Vertex_t> vertices; // vertices of one_skeleton. + + std::vector<std::vector<int> > simplices; // simplices of output simplicial complex. + std::vector<int> voronoi_subsamples; // Voronoi germs (in case of Voronoi cover). + + Persistence_diagram PD; + std::vector<double> distribution; + + std::vector<std::vector<int> > + cover; // function associating to each data point the vector of cover elements to which it belongs. + std::map<int, std::vector<int> > + cover_back; // inverse of cover, in order to get the data points associated to a specific cover element. + std::map<int, double> cover_std; // standard function (induced by func) used to compute the extended persistence + // diagram of the output simplicial complex. + std::map<int, int> + cover_fct; // integer-valued function that allows to state if two elements of the cover are consecutive or not. + std::map<int, std::pair<int, double> > + cover_color; // size and coloring (induced by func_color) of the vertices of the output simplicial complex. + + int resolution_int = -1; + double resolution_double = -1; + double gain = -1; + double rate_constant = 10; // Constant in the subsampling. + double rate_power = 0.001; // Power in the subsampling. + int mask = 0; // Ignore nodes containing less than mask points. + + std::map<int, int> name2id, name2idinv; + + std::string cover_name; + std::string point_cloud_name; + std::string color_name; + + // Remove all edges of a graph. + void remove_edges(Graph& G) { + boost::graph_traits<Graph>::edge_iterator ei, ei_end; + for (boost::tie(ei, ei_end) = boost::edges(G); ei != ei_end; ++ei) boost::remove_edge(*ei, G); + } + + // Thread local is not available on XCode version < V.8 + // If not available, random engine is a class member. +#ifndef GUDHI_CAN_USE_CXX11_THREAD_LOCAL + std::default_random_engine re; +#endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL + + // Find random number in [0,1]. + double GetUniform() { + // Thread local is not available on XCode version < V.8 + // If available, random engine is defined for each thread. +#ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL + thread_local std::default_random_engine re; +#endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL + std::uniform_real_distribution<double> Dist(0, 1); + return Dist(re); + } + + // Subsample points. + void SampleWithoutReplacement(int populationSize, int sampleSize, std::vector<int>& samples) { + int t = 0; + int m = 0; + double u; + while (m < sampleSize) { + u = GetUniform(); + if ((populationSize - t) * u >= sampleSize - m) { + t++; + } else { + samples[m] = t; + t++; + m++; + } + } + } + + // ******************************************************************************************************************* + // Utils. + // ******************************************************************************************************************* + + public: + /** \brief Specifies whether the type of the output simplicial complex. + * + * @param[in] t std::string (either "GIC" or "Nerve"). + * + */ + void set_type(const std::string& t) { type = t; } + + public: + /** \brief Specifies whether the program should display information or not. + * + * @param[in] verb boolean (true = display info, false = do not display info). + * + */ + void set_verbose(bool verb = false) { verbose = verb; } + + public: + /** \brief Sets the constants used to subsample the data set. These constants are + * explained in \cite Carriere17c. + * + * @param[in] constant double. + * @param[in] power double. + * + */ + void set_subsampling(double constant, double power) { + rate_constant = constant; + rate_power = power; + } + + public: + /** \brief Sets the mask, which is a threshold integer such that nodes in the complex that contain a number of data + * points which is less than or equal to + * this threshold are not displayed. + * + * @param[in] nodemask integer. + * + */ + void set_mask(int nodemask) { mask = nodemask; } + + public: + + + /** \brief Reads and stores the input point cloud from vector stored in memory. + * + * @param[in] point_cloud input vector representing the point cloud. Each row is a point and each coordinate is a vector. + * + */ + void set_point_cloud_from_range(const std::vector<std::vector<double> > & point_cloud) { + n = point_cloud.size(); data_dimension = point_cloud[0].size(); + point_cloud_name = "matrix"; cover.resize(n); + for(int i = 0; i < n; i++){ + boost::add_vertex(one_skeleton_OFF); + vertices.push_back(boost::add_vertex(one_skeleton)); + } + this->point_cloud = point_cloud; + } + + /** \brief Reads and stores the input point cloud from .(n)OFF file. + * + * @param[in] off_file_name name of the input .OFF or .nOFF file. + * + */ + bool read_point_cloud(const std::string& off_file_name) { + point_cloud_name = off_file_name; + std::ifstream input(off_file_name); + std::string line; + + char comment = '#'; + while (comment == '#') { + std::getline(input, line); + if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace)) + comment = line[line.find_first_not_of(' ')]; + } + if (strcmp((char*)line.c_str(), "nOFF") == 0) { + comment = '#'; + while (comment == '#') { + std::getline(input, line); + if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace)) + comment = line[line.find_first_not_of(' ')]; + } + std::stringstream stream(line); + stream >> data_dimension; + } else { + data_dimension = 3; + } + + comment = '#'; + int numedges, numfaces, i, dim; + while (comment == '#') { + std::getline(input, line); + if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace)) + comment = line[line.find_first_not_of(' ')]; + } + std::stringstream stream(line); + stream >> n; + stream >> numfaces; + stream >> numedges; + + i = 0; + while (i < n) { + std::getline(input, line); + if (!line.empty() && line[line.find_first_not_of(' ')] != '#' && + !all_of(line.begin(), line.end(), (int (*)(int))isspace)) { + std::stringstream iss(line); + std::vector<double> point; + point.assign(std::istream_iterator<double>(iss), std::istream_iterator<double>()); + point_cloud.emplace_back(point.begin(), point.begin() + data_dimension); + boost::add_vertex(one_skeleton_OFF); + vertices.push_back(boost::add_vertex(one_skeleton)); + cover.emplace_back(); + i++; + } + } + + i = 0; + while (i < numfaces) { + std::getline(input, line); + if (!line.empty() && line[line.find_first_not_of(' ')] != '#' && + !all_of(line.begin(), line.end(), (int (*)(int))isspace)) { + std::vector<int> simplex; + std::stringstream iss(line); + simplex.assign(std::istream_iterator<int>(iss), std::istream_iterator<int>()); + dim = simplex[0]; + for (int j = 1; j <= dim; j++) + for (int k = j + 1; k <= dim; k++) + boost::add_edge(vertices[simplex[j]], vertices[simplex[k]], one_skeleton_OFF); + i++; + } + } + + return input.is_open(); + } + + // ******************************************************************************************************************* + // Graphs. + // ******************************************************************************************************************* + + public: // Set graph from file. + /** \brief Creates a graph G from a file containing the edges. + * + * @param[in] graph_file_name name of the input graph file. + * The graph file contains one edge per line, + * each edge being represented by the IDs of its two nodes. + * + */ + void set_graph_from_file(const std::string& graph_file_name) { + remove_edges(one_skeleton); + int neighb; + std::ifstream input(graph_file_name); + std::string line; + int source; + while (std::getline(input, line)) { + std::stringstream stream(line); + stream >> source; + while (stream >> neighb) boost::add_edge(vertices[source], vertices[neighb], one_skeleton); + } + } + + public: // Set graph from OFF file. + /** \brief Creates a graph G from the triangulation given by the input .OFF file. + * + */ + void set_graph_from_OFF() { + remove_edges(one_skeleton); + if (num_edges(one_skeleton_OFF)) + one_skeleton = one_skeleton_OFF; + else + std::cout << "No triangulation read in OFF file!" << std::endl; + } + + public: // Set graph from Rips complex. + /** \brief Creates a graph G from a Rips complex. + * + * @param[in] threshold threshold value for the Rips complex. + * @param[in] distance distance used to compute the Rips complex. + * + */ + template <typename Distance> + void set_graph_from_rips(double threshold, Distance distance) { + remove_edges(one_skeleton); + if (distances.size() == 0) compute_pairwise_distances(distance); + for (int i = 0; i < n; i++) { + for (int j = i + 1; j < n; j++) { + if (distances[i][j] <= threshold) { + boost::add_edge(vertices[i], vertices[j], one_skeleton); + boost::put(boost::edge_weight, one_skeleton, boost::edge(vertices[i], vertices[j], one_skeleton).first, + distances[i][j]); + } + } + } + } + + public: + void set_graph_weights() { + Index_map index = boost::get(boost::vertex_index, one_skeleton); + Weight_map weight = boost::get(boost::edge_weight, one_skeleton); + boost::graph_traits<Graph>::edge_iterator ei, ei_end; + for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei) + boost::put(weight, *ei, + distances[index[boost::source(*ei, one_skeleton)]][index[boost::target(*ei, one_skeleton)]]); + } + + public: + /** \brief Reads and stores the distance matrices from vector stored in memory. + * + * @param[in] distance_matrix input vector representing the distance matrix. + * + */ + void set_distances_from_range(const std::vector<std::vector<double> > & distance_matrix) { + n = distance_matrix.size(); data_dimension = 0; point_cloud_name = "matrix"; + cover.resize(n); point_cloud.resize(n); + for(int i = 0; i < n; i++){ + boost::add_vertex(one_skeleton_OFF); + vertices.push_back(boost::add_vertex(one_skeleton)); + } + distances = distance_matrix; + } + + public: // Pairwise distances. + /** \private \brief Computes all pairwise distances. + */ + template <typename Distance> + void compute_pairwise_distances(Distance ref_distance) { + double d; + std::vector<double> zeros(n); + for (int i = 0; i < n; i++) distances.push_back(zeros); + std::string distance = point_cloud_name + "_dist"; + std::ifstream input(distance, std::ios::out | std::ios::binary); + + if (input.good()) { + if (verbose) std::cout << "Reading distances..." << std::endl; + for (int i = 0; i < n; i++) { + for (int j = i; j < n; j++) { + input.read((char*)&d, 8); + distances[i][j] = d; + distances[j][i] = d; + } + } + input.close(); + } else { + if (verbose) std::cout << "Computing distances..." << std::endl; + input.close(); + std::ofstream output(distance, std::ios::out | std::ios::binary); + for (int i = 0; i < n; i++) { + int state = (int)floor(100 * (i * 1.0 + 1) / n) % 10; + if (state == 0 && verbose) std::cout << "\r" << state << "%" << std::flush; + for (int j = i; j < n; j++) { + double dis = ref_distance(point_cloud[i], point_cloud[j]); + distances[i][j] = dis; + distances[j][i] = dis; + output.write((char*)&dis, 8); + } + } + output.close(); + if (verbose) std::cout << std::endl; + } + } + + public: // Automatic tuning of Rips complex. + /** \brief Creates a graph G from a Rips complex whose threshold value is automatically tuned with subsampling---see + * \cite Carriere17c. + * + * @param[in] distance distance between data points. + * @param[in] N number of subsampling iteration (the default reasonable value is 100, but there is no guarantee on + * how to choose it). + * @result delta threshold used for computing the Rips complex. + * + */ + template <typename Distance> + double set_graph_from_automatic_rips(Distance distance, int N = 100) { + int m = floor(n / std::exp((1 + rate_power) * std::log(std::log(n) / std::log(rate_constant)))); + m = (std::min)(m, n - 1); + double delta = 0; + + if (verbose) std::cout << n << " points in R^" << data_dimension << std::endl; + if (verbose) std::cout << "Subsampling " << m << " points" << std::endl; + + if (distances.size() == 0) compute_pairwise_distances(distance); + + // This cannot be parallelized if thread_local is not defined + // thread_local is not defined for XCode < v.8 + #if defined(GUDHI_USE_TBB) && defined(GUDHI_CAN_USE_CXX11_THREAD_LOCAL) + tbb::mutex deltamutex; + tbb::parallel_for(0, N, [&](int i){ + std::vector<int> samples(m); + SampleWithoutReplacement(n, m, samples); + double hausdorff_dist = 0; + for (int j = 0; j < n; j++) { + double mj = distances[j][samples[0]]; + for (int k = 1; k < m; k++) mj = (std::min)(mj, distances[j][samples[k]]); + hausdorff_dist = (std::max)(hausdorff_dist, mj); + } + deltamutex.lock(); + delta += hausdorff_dist / N; + deltamutex.unlock(); + }); + #else + for (int i = 0; i < N; i++) { + std::vector<int> samples(m); + SampleWithoutReplacement(n, m, samples); + double hausdorff_dist = 0; + for (int j = 0; j < n; j++) { + double mj = distances[j][samples[0]]; + for (int k = 1; k < m; k++) mj = (std::min)(mj, distances[j][samples[k]]); + hausdorff_dist = (std::max)(hausdorff_dist, mj); + } + delta += hausdorff_dist / N; + } + #endif + + if (verbose) std::cout << "delta = " << delta << std::endl; + set_graph_from_rips(delta, distance); + return delta; + } + + // ******************************************************************************************************************* + // Functions. + // ******************************************************************************************************************* + + public: // Set function from file. + /** \brief Creates the function f from a file containing the function values. + * + * @param[in] func_file_name name of the input function file. + * + */ + void set_function_from_file(const std::string& func_file_name) { + int i = 0; + std::ifstream input(func_file_name); + std::string line; + double f; + while (std::getline(input, line)) { + std::stringstream stream(line); + stream >> f; + func.push_back(f); + i++; + } + functional_cover = true; + cover_name = func_file_name; + } + + public: // Set function from kth coordinate + /** \brief Creates the function f from the k-th coordinate of the point cloud P. + * + * @param[in] k coordinate to use (start at 0). + * + */ + void set_function_from_coordinate(int k) { + if(point_cloud[0].size() > 0){ + for (int i = 0; i < n; i++) func.push_back(point_cloud[i][k]); + functional_cover = true; + cover_name = "coordinate " + std::to_string(k); + } + else{ + std::cout << "Only pairwise distances provided---cannot access " << k << "th coordinate; returning null vector instead" << std::endl; + for (int i = 0; i < n; i++) func.push_back(0.0); + functional_cover = true; + cover_name = "null"; + } + } + + public: // Set function from vector. + /** \brief Creates the function f from a vector stored in memory. + * + * @param[in] function input vector of values. + * + */ + template <class InputRange> + void set_function_from_range(InputRange const& function) { + for (int i = 0; i < n; i++) func.push_back(function[i]); + functional_cover = true; + } + + // ******************************************************************************************************************* + // Covers. + // ******************************************************************************************************************* + + public: // Automatic tuning of resolution. + /** \brief Computes the optimal length of intervals + * (i.e. the smallest interval length avoiding discretization artifacts---see \cite Carriere17c) for a functional + * cover. + * + * @result reso interval length used to compute the cover. + * + */ + double set_automatic_resolution() { + if (!functional_cover) { + std::cout << "Cover needs to come from the preimages of a function." << std::endl; + return 0; + } + if (type != "Nerve" && type != "GIC") { + std::cout << "Type of complex needs to be specified." << std::endl; + return 0; + } + + double reso = 0; + Index_map index = boost::get(boost::vertex_index, one_skeleton); + + if (type == "GIC") { + boost::graph_traits<Graph>::edge_iterator ei, ei_end; + for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei) + reso = (std::max)(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] - + func[index[boost::target(*ei, one_skeleton)]])); + if (verbose) std::cout << "resolution = " << reso << std::endl; + resolution_double = reso; + } + + if (type == "Nerve") { + boost::graph_traits<Graph>::edge_iterator ei, ei_end; + for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei) + reso = (std::max)(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] - + func[index[boost::target(*ei, one_skeleton)]]) / + gain); + if (verbose) std::cout << "resolution = " << reso << std::endl; + resolution_double = reso; + } + + return reso; + } + + public: + /** \brief Sets a length of intervals from a value stored in memory. + * + * @param[in] reso length of intervals. + * + */ + void set_resolution_with_interval_length(double reso) { resolution_double = reso; } + /** \brief Sets a number of intervals from a value stored in memory. + * + * @param[in] reso number of intervals. + * + */ + void set_resolution_with_interval_number(int reso) { resolution_int = reso; } + /** \brief Sets a gain from a value stored in memory (default value 0.3). + * + * @param[in] g gain. + * + */ + void set_gain(double g = 0.3) { gain = g; } + + public: // Set cover with preimages of function. + /** \brief Creates a cover C from the preimages of the function f. + * + */ + void set_cover_from_function() { + if (resolution_double == -1 && resolution_int == -1) { + std::cout << "Number and/or length of intervals not specified" << std::endl; + return; + } + if (gain == -1) { + std::cout << "Gain not specified" << std::endl; + return; + } + + // Read function values and compute min and max + double minf = (std::numeric_limits<float>::max)(); + double maxf = std::numeric_limits<float>::lowest(); + for (int i = 0; i < n; i++) { + minf = (std::min)(minf, func[i]); + maxf = (std::max)(maxf, func[i]); + } + if (verbose) std::cout << "Min function value = " << minf << " and Max function value = " << maxf << std::endl; + + // Compute cover of im(f) + std::vector<std::pair<double, double> > intervals; + int res; + + if (resolution_double == -1) { // Case we use an integer for the number of intervals. + double incr = (maxf - minf) / resolution_int; + 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_int - 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_int - 1) * incr - alpha; + y = maxf; + std::pair<double, double> interM(x, y); + intervals.push_back(interM); + res = intervals.size(); + if (verbose) { + for (int i = 0; i < res; i++) + std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" + << std::endl; + } + } else { + if (resolution_int == -1) { // Case we use a double for the length of the intervals. + double x = minf; + double y = x + resolution_double; + while (y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) { + std::pair<double, double> inter(x, y); + intervals.push_back(inter); + x = y - gain * resolution_double; + y = x + resolution_double; + } + std::pair<double, double> interM(x, maxf); + intervals.push_back(interM); + res = intervals.size(); + if (verbose) { + for (int i = 0; i < res; i++) + std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" + << std::endl; + } + } else { // Case we use an integer and a double for the length of the intervals. + double x = minf; + double y = x + resolution_double; + int count = 0; + while (count < resolution_int && y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) { + std::pair<double, double> inter(x, y); + intervals.push_back(inter); + count++; + x = y - gain * resolution_double; + y = x + resolution_double; + } + res = intervals.size(); + if (verbose) { + for (int i = 0; i < res; i++) + std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" + << std::endl; + } + } + } + + // Sort points according to function values + std::vector<int> points(n); + for (int i = 0; i < n; i++) points[i] = i; + std::sort(points.begin(), points.end(), [=](const int & p1, const int & p2){return (this->func[p1] < this->func[p2]);}); + + int id = 0; + int pos = 0; + Index_map index = boost::get(boost::vertex_index, one_skeleton); // int maxc = -1; + std::map<int, std::vector<int> > preimages; + std::map<int, double> funcstd; + + if (verbose) std::cout << "Computing preimages..." << std::endl; + for (int i = 0; i < res; i++) { + // Find points in the preimage + std::pair<double, double> inter1 = intervals[i]; + int tmp = pos; + double u, v; + + if (i != res - 1) { + if (i != 0) { + std::pair<double, double> inter3 = intervals[i - 1]; + while (func[points[tmp]] < inter3.second && tmp != n) { + preimages[i].push_back(points[tmp]); + tmp++; + } + u = inter3.second; + } else { + u = inter1.first; + } + + std::pair<double, double> inter2 = intervals[i + 1]; + 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) { + preimages[i].push_back(points[tmp]); + tmp++; + } + + } else { + std::pair<double, double> 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; + } + + funcstd[i] = 0.5 * (u + v); + } + + #ifdef GUDHI_USE_TBB + if (verbose) std::cout << "Computing connected components (parallelized)..." << std::endl; + tbb::mutex covermutex, idmutex; + tbb::parallel_for(0, res, [&](int i){ + // Compute connected components + Graph G = one_skeleton.create_subgraph(); + int num = preimages[i].size(); + std::vector<int> 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; + + // For each point in preimage + for (int j = 0; j < num; j++) { + // Update number of components in preimage + if (component[j] > max) max = component[j]; + + // Identify component with Cantor polynomial N^2 -> N + int identifier = ((i + component[j])*(i + component[j]) + 3 * i + component[j]) / 2; + + // Update covers + covermutex.lock(); + 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; + covermutex.unlock(); + } + + // Maximal dimension is total number of connected components + idmutex.lock(); + id += max + 1; + idmutex.unlock(); + }); + #else + if (verbose) std::cout << "Computing connected components..." << std::endl; + for (int i = 0; i < res; i++) { + // Compute connected components + Graph G = one_skeleton.create_subgraph(); + int num = preimages[i].size(); + std::vector<int> 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; + + // For each point in preimage + for (int j = 0; j < num; j++) { + // 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; + } + #endif + + maximal_dim = id - 1; + for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) + iit->second.second /= iit->second.first; + } + + public: // Set cover from file. + /** \brief Creates the cover C from a file containing the cover elements of each point (the order has to be the same + * as in the input file!). + * + * @param[in] cover_file_name name of the input cover file. + * + */ + void set_cover_from_file(const std::string& cover_file_name) { + int i = 0; + int cov; + std::vector<int> 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_fct[cov] = cov; + cover_color[cov].second += func_color[i]; + cover_color[cov].first++; + cover_back[cov].push_back(i); + } + cover[i] = cov_elts; + i++; + } + + std::sort(cov_number.begin(), cov_number.end()); + std::vector<int>::iterator it = std::unique(cov_number.begin(), cov_number.end()); + cov_number.resize(std::distance(cov_number.begin(), it)); + + maximal_dim = cov_number.size() - 1; + for (int i = 0; i <= maximal_dim; i++) cover_color[i].second /= cover_color[i].first; + cover_name = cover_file_name; + } + + public: // Set cover from Voronoi + /** \brief Creates the cover C from the Voronoï cells of a subsampling of the point cloud. + * + * @param[in] distance distance between the points. + * @param[in] m number of points in the subsample. + * + */ + template <typename Distance> + void set_cover_from_Voronoi(Distance distance, int m = 100) { + voronoi_subsamples.resize(m); + SampleWithoutReplacement(n, m, voronoi_subsamples); + if (distances.size() == 0) compute_pairwise_distances(distance); + set_graph_weights(); + Weight_map weight = boost::get(boost::edge_weight, one_skeleton); + Index_map index = boost::get(boost::vertex_index, one_skeleton); + std::vector<double> mindist(n); + for (int j = 0; j < n; j++) mindist[j] = (std::numeric_limits<double>::max)(); + + // Compute the geodesic distances to subsamples with Dijkstra + #ifdef GUDHI_USE_TBB + if (verbose) std::cout << "Computing geodesic distances (parallelized)..." << std::endl; + tbb::mutex coverMutex; tbb::mutex mindistMutex; + tbb::parallel_for(0, m, [&](int i){ + int seed = voronoi_subsamples[i]; + std::vector<double> dmap(n); + boost::dijkstra_shortest_paths( + one_skeleton, vertices[seed], + boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index))); + + coverMutex.lock(); mindistMutex.lock(); + for (int j = 0; j < n; j++) + if (mindist[j] > dmap[j]) { + mindist[j] = dmap[j]; + if (cover[j].size() == 0) + cover[j].push_back(i); + else + cover[j][0] = i; + } + coverMutex.unlock(); mindistMutex.unlock(); + }); + #else + for (int i = 0; i < m; i++) { + if (verbose) std::cout << "Computing geodesic distances to seed " << i << "..." << std::endl; + int seed = voronoi_subsamples[i]; + std::vector<double> dmap(n); + boost::dijkstra_shortest_paths( + one_skeleton, vertices[seed], + boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index))); + + for (int j = 0; j < n; j++) + if (mindist[j] > dmap[j]) { + mindist[j] = dmap[j]; + if (cover[j].size() == 0) + cover[j].push_back(i); + else + cover[j][0] = i; + } + } + #endif + + for (int i = 0; i < n; i++) { + cover_back[cover[i][0]].push_back(i); + cover_color[cover[i][0]].second += func_color[i]; + cover_color[cover[i][0]].first++; + } + for (int i = 0; i < m; i++) cover_color[i].second /= cover_color[i].first; + maximal_dim = m - 1; + cover_name = "Voronoi"; + } + + public: // return subset of data corresponding to a node + /** \brief Returns the data subset corresponding to a specific node of the created complex. + * + * @param[in] c ID of the node. + * @result cover_back(c) vector of IDs of data points. + * + */ + const std::vector<int>& subpopulation(int c) { return cover_back[name2idinv[c]]; } + + // ******************************************************************************************************************* + // Visualization. + // ******************************************************************************************************************* + + public: // Set color from file. + /** \brief Computes the function used to color the nodes of the simplicial complex from a file containing the function + * values. + * + * @param[in] color_file_name name of the input color file. + * + */ + void set_color_from_file(const std::string& color_file_name) { + int i = 0; + std::ifstream input(color_file_name); + std::string line; + double f; + while (std::getline(input, line)) { + std::stringstream stream(line); + stream >> f; + func_color.push_back(f); + i++; + } + color_name = color_file_name; + } + + public: // Set color from kth coordinate + /** \brief Computes the function used to color the nodes of the simplicial complex from the k-th coordinate. + * + * @param[in] k coordinate to use (start at 0). + * + */ + void set_color_from_coordinate(int k = 0) { + if(point_cloud[0].size() > 0){ + for (int i = 0; i < n; i++) func_color.push_back(point_cloud[i][k]); + color_name = "coordinate "; + color_name.append(std::to_string(k)); + } + else{ + std::cout << "Only pairwise distances provided---cannot access " << k << "th coordinate; returning null vector instead" << std::endl; + for (int i = 0; i < n; i++) func.push_back(0.0); + functional_cover = true; + cover_name = "null"; + } + } + + public: // Set color from vector. + /** \brief Computes the function used to color the nodes of the simplicial complex from a vector stored in memory. + * + * @param[in] color input vector of values. + * + */ + void set_color_from_range(std::vector<double> color) { + for (unsigned int i = 0; i < color.size(); i++) func_color.push_back(color[i]); + } + + public: // Create a .dot file that can be compiled with neato to produce a .pdf file. + /** \brief Creates a .dot file called SC.dot for neato (part of the graphviz package) once the simplicial complex is + * computed to get a visualization + * of its 1-skeleton in a .pdf file. + */ + void plot_DOT() { + std::string mapp = point_cloud_name + "_sc.dot"; + std::ofstream graphic(mapp); + + double maxv = std::numeric_limits<double>::lowest(); + double minv = (std::numeric_limits<double>::max)(); + for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) { + maxv = (std::max)(maxv, iit->second.second); + minv = (std::min)(minv, iit->second.second); + } + + int k = 0; + std::vector<int> nodes; + nodes.clear(); + + graphic << "graph GIC {" << std::endl; + int id = 0; + for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) { + if (iit->second.first > mask) { + 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++; + } + } + int ke = 0; + int num_simplices = simplices.size(); + 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 << " " << name2id[simplices[i][0]] << " -- " << name2id[simplices[i][1]] << " [weight=15];" + << std::endl; + ke++; + } + } + graphic << "}"; + graphic.close(); + std::cout << mapp << " file generated. It can be visualized with e.g. neato." << std::endl; + } + + public: // Create a .txt file that can be compiled with KeplerMapper. + /** \brief Creates a .txt file called SC.txt describing the 1-skeleton, which can then be plotted with e.g. + * KeplerMapper. + */ + void write_info() { + int num_simplices = simplices.size(); + int num_edges = 0; + std::string mapp = point_cloud_name + "_sc.txt"; + std::ofstream graphic(mapp); + + 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) num_edges++; + + graphic << point_cloud_name << std::endl; + graphic << cover_name << std::endl; + graphic << color_name << std::endl; + graphic << resolution_double << " " << gain << std::endl; + graphic << cover_color.size() << " " << num_edges << std::endl; + + int id = 0; + for (std::map<int, std::pair<int, double> >::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 << name2id[simplices[i][0]] << " " << name2id[simplices[i][1]] << std::endl; + graphic.close(); + std::cout << mapp + << " generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox." + << std::endl; + } + + public: // Create a .off file that can be visualized (e.g. with Geomview). + /** \brief Creates a .off file called SC.off for 3D visualization, which contains the 2-skeleton of the GIC. + * This function assumes that the cover has been computed with Voronoi. If data points are in 1D or 2D, + * the remaining coordinates of the points embedded in 3D are set to 0. + */ + void plot_OFF() { + assert(cover_name == "Voronoi"); + + int m = voronoi_subsamples.size(); + int numedges = 0; + int numfaces = 0; + std::vector<std::vector<int> > edges, faces; + int numsimplices = simplices.size(); + + std::string mapp = point_cloud_name + "_sc.off"; + std::ofstream graphic(mapp); + + graphic << "OFF" << std::endl; + for (int i = 0; i < numsimplices; i++) { + if (simplices[i].size() == 2) { + numedges++; + edges.push_back(simplices[i]); + } + if (simplices[i].size() == 3) { + numfaces++; + faces.push_back(simplices[i]); + } + } + graphic << m << " " << numedges + numfaces << std::endl; + for (int i = 0; i < m; i++) { + if (data_dimension <= 3) { + for (int j = 0; j < data_dimension; j++) graphic << point_cloud[voronoi_subsamples[i]][j] << " "; + for (int j = data_dimension; j < 3; j++) graphic << 0 << " "; + graphic << std::endl; + } else { + for (int j = 0; j < 3; j++) graphic << point_cloud[voronoi_subsamples[i]][j] << " "; + } + } + for (int i = 0; i < numedges; i++) graphic << 2 << " " << edges[i][0] << " " << edges[i][1] << std::endl; + for (int i = 0; i < numfaces; i++) + graphic << 3 << " " << faces[i][0] << " " << faces[i][1] << " " << faces[i][2] << std::endl; + graphic.close(); + std::cout << mapp << " generated. It can be visualized with e.g. geomview." << std::endl; + } + + // ******************************************************************************************************************* + // Extended Persistence Diagrams. + // ******************************************************************************************************************* + + public: + /** \brief Computes the extended persistence diagram of the complex. + * + */ + Persistence_diagram compute_PD() { + Simplex_tree st; + + // Compute max and min + double maxf = std::numeric_limits<double>::lowest(); + double minf = (std::numeric_limits<double>::max)(); + for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) { + maxf = (std::max)(maxf, it->second); + minf = (std::min)(minf, it->second); + } + + // Build filtration + for (auto const& simplex : simplices) { + std::vector<int> splx = simplex; + splx.push_back(-2); + st.insert_simplex_and_subfaces(splx, -3); + } + + for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) { + int vertex = it->first; float val = it->second; + int vert[] = {vertex}; int edge[] = {vertex, -2}; + if(st.find(vert) != st.null_simplex()){ + st.assign_filtration(st.find(vert), -2 + (val - minf)/(maxf - minf)); + st.assign_filtration(st.find(edge), 2 - (val - minf)/(maxf - minf)); + } + } + st.make_filtration_non_decreasing(); + + // Compute PD + Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(st); pcoh.init_coefficients(2); + pcoh.compute_persistent_cohomology(); + + // Output PD + int max_dim = st.dimension(); + for (int i = 0; i < max_dim; i++) { + std::vector<std::pair<double, double> > bars = pcoh.intervals_in_dimension(i); + int num_bars = bars.size(); if(i == 0) num_bars -= 1; + if(verbose) 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 (i == 0 && std::isinf(death)) 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); + PD.push_back(std::pair<double, double>(birth, death)); + if (verbose) std::cout << " [" << birth << ", " << death << "]" << std::endl; + } + } + return PD; + } + + public: + /** \brief Computes bootstrapped distances distribution. + * + * @param[in] N number of bootstrap iterations. + * + */ + void compute_distribution(unsigned int N = 100) { + unsigned int sz = distribution.size(); + if (sz >= N) { + std::cout << "Already done!" << std::endl; + } else { + for (unsigned int i = 0; i < N - sz; i++) { + if (verbose) std::cout << "Computing " << i << "th bootstrap, bottleneck distance = "; + + Cover_complex Cboot; Cboot.n = this->n; Cboot.data_dimension = this->data_dimension; Cboot.type = this->type; Cboot.functional_cover = true; + + std::vector<int> boot(this->n); + for (int j = 0; j < this->n; j++) { + double u = GetUniform(); + int id = std::floor(u * (this->n)); boot[j] = id; + Cboot.point_cloud.push_back(this->point_cloud[id]); Cboot.cover.emplace_back(); Cboot.func.push_back(this->func[id]); + boost::add_vertex(Cboot.one_skeleton_OFF); Cboot.vertices.push_back(boost::add_vertex(Cboot.one_skeleton)); + } + Cboot.set_color_from_range(Cboot.func); + + for (int j = 0; j < n; j++) { + std::vector<double> dist(n); + for (int k = 0; k < n; k++) dist[k] = distances[boot[j]][boot[k]]; + Cboot.distances.push_back(dist); + } + + Cboot.set_graph_from_automatic_rips(Gudhi::Euclidean_distance()); + Cboot.set_gain(); + Cboot.set_automatic_resolution(); + Cboot.set_cover_from_function(); + Cboot.find_simplices(); + Cboot.compute_PD(); + double db = Gudhi::persistence_diagram::bottleneck_distance(this->PD, Cboot.PD); + if (verbose) std::cout << db << std::endl; + distribution.push_back(db); + } + + std::sort(distribution.begin(), distribution.end()); + } + } + + public: + /** \brief Computes the bottleneck distance threshold corresponding to a specific confidence level. + * + * @param[in] alpha Confidence level. + * + */ + double compute_distance_from_confidence_level(double alpha) { + unsigned int N = distribution.size(); + double d = distribution[std::floor(alpha * N)]; + if (verbose) std::cout << "Distance corresponding to confidence " << alpha << " is " << d << std::endl; + return d; + } + + public: + /** \brief Computes the confidence level of a specific bottleneck distance threshold. + * + * @param[in] d Bottleneck distance. + * + */ + double compute_confidence_level_from_distance(double d) { + unsigned int N = distribution.size(); + double level = 1; + for (unsigned int i = 0; i < N; i++) + if (distribution[i] >= d){ level = i * 1.0 / N; break; } + if (verbose) std::cout << "Confidence level of distance " << d << " is " << level << std::endl; + return level; + } + + public: + /** \brief Computes the p-value, i.e. the opposite of the confidence level of the largest bottleneck + * distance preserving the points in the persistence diagram of the output simplicial complex. + * + */ + double compute_p_value() { + double distancemin = (std::numeric_limits<double>::max)(); int N = PD.size(); + for (int i = 0; i < N; i++) distancemin = (std::min)(distancemin, 0.5 * std::abs(PD[i].second - PD[i].first)); + double p_value = 1 - compute_confidence_level_from_distance(distancemin); + if (verbose) std::cout << "p value = " << p_value << std::endl; + return p_value; + } + + // ******************************************************************************************************************* + // Computation of simplices. + // ******************************************************************************************************************* + + public: + /** \brief Creates the simplicial complex. + * + * @param[in] complex SimplicialComplex to be created. + * + */ + template <typename SimplicialComplex> + void create_complex(SimplicialComplex& complex) { + unsigned int dimension = 0; + for (auto const& simplex : simplices) { + int numvert = simplex.size(); + double filt = std::numeric_limits<double>::lowest(); + for (int i = 0; i < numvert; i++) filt = (std::max)(cover_color[simplex[i]].second, filt); + complex.insert_simplex_and_subfaces(simplex, filt); + if (dimension < simplex.size() - 1) dimension = simplex.size() - 1; + } + } + + public: + /** \brief Computes the simplices of the simplicial complex. + */ + void find_simplices() { + if (type != "Nerve" && type != "GIC") { + std::cout << "Type of complex needs to be specified." << std::endl; + return; + } + + if (type == "Nerve") { + for(int i = 0; i < n; i++) simplices.push_back(cover[i]); + std::sort(simplices.begin(), simplices.end()); + std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end()); + simplices.resize(std::distance(simplices.begin(), it)); + } + + if (type == "GIC") { + Index_map index = boost::get(boost::vertex_index, one_skeleton); + + if (functional_cover) { + // Computes the simplices in the GIC by looking at all the edges of the graph and adding the + // corresponding edges in the GIC if the images of the endpoints belong to consecutive intervals. + + if (gain >= 0.5) + throw std::invalid_argument( + "the output of this function is correct ONLY if the cover is minimal, i.e. the gain is less than 0.5."); + + // Loop on all edges. + boost::graph_traits<Graph>::edge_iterator ei, ei_end; + for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei) { + int nums = cover[index[boost::source(*ei, one_skeleton)]].size(); + for (int i = 0; i < nums; i++) { + int vs = cover[index[boost::source(*ei, one_skeleton)]][i]; + int numt = cover[index[boost::target(*ei, one_skeleton)]].size(); + 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<int> edge(2); + edge[0] = (std::min)(vs, vt); + edge[1] = (std::max)(vs, vt); + simplices.push_back(edge); + goto afterLoop; + } + } + } + afterLoop:; + } + std::sort(simplices.begin(), simplices.end()); + std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end()); + simplices.resize(std::distance(simplices.begin(), it)); + + } else { + // Find edges to keep + Simplex_tree st; + boost::graph_traits<Graph>::edge_iterator ei, ei_end; + for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei) + if (!(cover[index[boost::target(*ei, one_skeleton)]].size() == 1 && + cover[index[boost::target(*ei, one_skeleton)]] == cover[index[boost::source(*ei, one_skeleton)]])) { + std::vector<int> edge(2); + edge[0] = index[boost::source(*ei, one_skeleton)]; + edge[1] = index[boost::target(*ei, one_skeleton)]; + st.insert_simplex_and_subfaces(edge); + } + + // st.insert_graph(one_skeleton); + + // 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<int> simplx; + for (auto vertex : st.simplex_vertex_range(simplex)) { + unsigned int sz = cover[vertex].size(); + for (unsigned int i = 0; i < sz; i++) { + simplx.push_back(cover[vertex][i]); + } + } + std::sort(simplx.begin(), simplx.end()); + std::vector<int>::iterator it = std::unique(simplx.begin(), simplx.end()); + simplx.resize(std::distance(simplx.begin(), it)); + simplices.push_back(simplx); + } + } + std::sort(simplices.begin(), simplices.end()); + std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end()); + simplices.resize(std::distance(simplices.begin(), it)); + } + } + } +}; + +} // namespace cover_complex + +} // namespace Gudhi + +#endif // GIC_H_ diff --git a/src/Nerve_GIC/test/CMakeLists.txt b/src/Nerve_GIC/test/CMakeLists.txt new file mode 100644 index 00000000..b89c18a2 --- /dev/null +++ b/src/Nerve_GIC/test/CMakeLists.txt @@ -0,0 +1,16 @@ +project(Graph_induced_complex_tests) + +if (NOT CGAL_VERSION VERSION_LESS 4.11.0) + include(GUDHI_test_coverage) + + add_executable ( Nerve_GIC_test_unit test_GIC.cpp ) + target_link_libraries(Nerve_GIC_test_unit ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(Nerve_GIC_test_unit ${TBB_LIBRARIES}) + endif() + + file(COPY data DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + + gudhi_add_coverage_test(Nerve_GIC_test_unit) + +endif (NOT CGAL_VERSION VERSION_LESS 4.11.0) diff --git a/src/Nerve_GIC/test/data/cloud b/src/Nerve_GIC/test/data/cloud new file mode 100644 index 00000000..4a0c170d --- /dev/null +++ b/src/Nerve_GIC/test/data/cloud @@ -0,0 +1,6 @@ +nOFF +3 +3 0 0 +0 0 0 +2 1 0 +4 0 0
\ No newline at end of file diff --git a/src/Nerve_GIC/test/data/cover b/src/Nerve_GIC/test/data/cover new file mode 100644 index 00000000..5f5fbe75 --- /dev/null +++ b/src/Nerve_GIC/test/data/cover @@ -0,0 +1,3 @@ +1 +2 +3
\ No newline at end of file diff --git a/src/Nerve_GIC/test/data/graph b/src/Nerve_GIC/test/data/graph new file mode 100644 index 00000000..37142800 --- /dev/null +++ b/src/Nerve_GIC/test/data/graph @@ -0,0 +1,3 @@ +0 1 +0 2 +1 2
\ No newline at end of file diff --git a/src/Nerve_GIC/test/test_GIC.cpp b/src/Nerve_GIC/test/test_GIC.cpp new file mode 100644 index 00000000..9d326234 --- /dev/null +++ b/src/Nerve_GIC/test/test_GIC.cpp @@ -0,0 +1,81 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Mathieu Carrière + * + * Copyright (C) 2017 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "graph_induced_complex" + +#include <boost/test/unit_test.hpp> + +#include <limits> +#include <string> +#include <vector> + +#include <gudhi/GIC.h> +#include <gudhi/distance_functions.h> +#include <gudhi/reader_utils.h> + +BOOST_AUTO_TEST_CASE(check_nerve) { + using Point = std::vector<float>; + Gudhi::cover_complex::Cover_complex<Point> N; + N.set_type("Nerve"); + std::string cloud_file_name("data/cloud"); + N.read_point_cloud(cloud_file_name); + N.set_color_from_coordinate(); + std::string graph_file_name("data/graph"); + N.set_graph_from_file(graph_file_name); + std::string cover_file_name("data/cover"); + N.set_cover_from_file(cover_file_name); + N.find_simplices(); + Gudhi::Simplex_tree<> stree; + N.create_complex(stree); + + BOOST_CHECK(stree.num_vertices() == 3); + BOOST_CHECK((stree.num_simplices() - stree.num_vertices()) == 0); + BOOST_CHECK(stree.dimension() == 0); +} + +BOOST_AUTO_TEST_CASE(check_GIC) { + using Point = std::vector<float>; + Gudhi::cover_complex::Cover_complex<Point> GIC; + GIC.set_type("GIC"); + std::string cloud_file_name("data/cloud"); + GIC.read_point_cloud(cloud_file_name); + GIC.set_color_from_coordinate(); + std::string graph_file_name("data/graph"); + GIC.set_graph_from_file(graph_file_name); + std::string cover_file_name("data/cover"); + GIC.set_cover_from_file(cover_file_name); + GIC.find_simplices(); + Gudhi::Simplex_tree<> stree; + GIC.create_complex(stree); + + BOOST_CHECK(stree.num_vertices() == 3); + BOOST_CHECK((stree.num_simplices() - stree.num_vertices()) == 4); + BOOST_CHECK(stree.dimension() == 2); +} + +BOOST_AUTO_TEST_CASE(check_voronoiGIC) { + using Point = std::vector<float>; + Gudhi::cover_complex::Cover_complex<Point> GIC; + GIC.set_type("GIC"); + std::string cloud_file_name("data/cloud"); + GIC.read_point_cloud(cloud_file_name); + GIC.set_color_from_coordinate(); + std::string graph_file_name("data/graph"); + GIC.set_graph_from_file(graph_file_name); + GIC.set_cover_from_Voronoi(Gudhi::Euclidean_distance(), 2); + GIC.find_simplices(); + Gudhi::Simplex_tree<> stree; + GIC.create_complex(stree); + + BOOST_CHECK(stree.num_vertices() == 2); + BOOST_CHECK((stree.num_simplices() - stree.num_vertices()) == 1); + BOOST_CHECK(stree.dimension() == 1); +} diff --git a/src/Nerve_GIC/utilities/CMakeLists.txt b/src/Nerve_GIC/utilities/CMakeLists.txt new file mode 100644 index 00000000..65a08d9a --- /dev/null +++ b/src/Nerve_GIC/utilities/CMakeLists.txt @@ -0,0 +1,27 @@ +project(Nerve_GIC_examples) + +if (NOT CGAL_VERSION VERSION_LESS 4.11.0) + + add_executable ( Nerve Nerve.cpp ) + add_executable ( VoronoiGIC VoronoiGIC.cpp ) + + if (TBB_FOUND) + target_link_libraries(Nerve ${TBB_LIBRARIES}) + target_link_libraries(VoronoiGIC ${TBB_LIBRARIES}) + endif() + + file(COPY KeplerMapperVisuFromTxtFile.py km.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + # Copy files for not to pollute sources when testing + file(COPY "${CMAKE_SOURCE_DIR}/data/points/human.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + + add_test(NAME Nerve_GIC_utilities_nerve COMMAND $<TARGET_FILE:Nerve> + "human.off" "2" "10" "0.3") + + add_test(NAME Nerve_GIC_utilities_VoronoiGIC COMMAND $<TARGET_FILE:VoronoiGIC> + "human.off" "100") + + install(TARGETS Nerve DESTINATION bin) + install(TARGETS VoronoiGIC DESTINATION bin) + install(FILES KeplerMapperVisuFromTxtFile.py km.py km.py.COPYRIGHT DESTINATION bin) + +endif (NOT CGAL_VERSION VERSION_LESS 4.11.0) diff --git a/src/Nerve_GIC/utilities/KeplerMapperVisuFromTxtFile.py b/src/Nerve_GIC/utilities/KeplerMapperVisuFromTxtFile.py new file mode 100755 index 00000000..e3101549 --- /dev/null +++ b/src/Nerve_GIC/utilities/KeplerMapperVisuFromTxtFile.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python + +import km +import numpy as np +from collections import defaultdict +import argparse + +"""This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + Author(s): Mathieu Carriere + + Copyright (C) 2017 Inria + + Modification(s): + - YYYY/MM Author: Description of the modification +""" + +__author__ = "Mathieu Carriere" +__copyright__ = "Copyright (C) 2017 Inria" +__license__ = "GPL v3" + +parser = argparse.ArgumentParser(description='Creates an html Keppler Mapper ' + 'file to visualize a SC.txt file.', + epilog='Example: ' + './KeplerMapperVisuFromTxtFile.py ' + '-f ../../data/points/human.off_sc.txt' + '- Constructs an human.off_sc.html file.') +parser.add_argument("-f", "--file", type=str, required=True) + +args = parser.parse_args() + +with open(args.file, 'r') as f: + network = {} + mapper = km.KeplerMapper(verbose=0) + data = np.zeros((3,3)) + projected_data = mapper.fit_transform( data, projection="sum", scaler=None ) + + nodes = defaultdict(list) + links = defaultdict(list) + custom = defaultdict(list) + + dat = f.readline() + lens = f.readline() + color = f.readline(); + param = [float(i) for i in f.readline().split(" ")] + + nums = [int(i) for i in f.readline().split(" ")] + num_nodes = nums[0] + num_edges = nums[1] + + for i in range(0,num_nodes): + point = [float(j) for j in f.readline().split(" ")] + nodes[ str(int(point[0])) ] = [ int(point[0]), point[1], int(point[2]) ] + links[ str(int(point[0])) ] = [] + custom[ int(point[0]) ] = point[1] + + m = min([custom[i] for i in range(0,num_nodes)]) + M = max([custom[i] for i in range(0,num_nodes)]) + + for i in range(0,num_edges): + edge = [int(j) for j in f.readline().split(" ")] + links[ str(edge[0]) ].append( str(edge[1]) ) + links[ str(edge[1]) ].append( str(edge[0]) ) + + network["nodes"] = nodes + network["links"] = links + network["meta"] = lens + + html_output_filename = args.file.rsplit('.', 1)[0] + '.html' + mapper.visualize(network, color_function = color, path_html=html_output_filename, title=dat, + graph_link_distance=30, graph_gravity=0.1, graph_charge=-120, custom_tooltips=custom, width_html=0, + height_html=0, show_tooltips=True, show_title=True, show_meta=True, res=param[0],gain=param[1], minimum=m,maximum=M) + message = repr(html_output_filename) + " is generated. You can now use your favorite web browser to visualize it." + print(message) + + + f.close() diff --git a/src/Nerve_GIC/utilities/Nerve.cpp b/src/Nerve_GIC/utilities/Nerve.cpp new file mode 100644 index 00000000..d34e922c --- /dev/null +++ b/src/Nerve_GIC/utilities/Nerve.cpp @@ -0,0 +1,84 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Mathieu Carrière + * + * Copyright (C) 2017 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/GIC.h> + +#include <string> +#include <vector> + +void usage(int nbArgs, char *const progName) { + std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; + std::cerr << "Usage: " << progName << " filename.off coordinate resolution gain [-v] \n"; + std::cerr << " i.e.: " << progName << " ../../data/points/human.off 2 10 0.3 -v \n"; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if ((argc != 5) && (argc != 6)) usage(argc, argv[0]); + + using Point = std::vector<float>; + + std::string off_file_name(argv[1]); + int coord = atoi(argv[2]); + int resolution = atoi(argv[3]); + double gain = atof(argv[4]); + bool verb = 0; + if (argc == 6) verb = 1; + + // -------------------------------- + // Init of a Nerve from an OFF file + // -------------------------------- + + Gudhi::cover_complex::Cover_complex<Point> SC; + SC.set_verbose(verb); + + bool check = SC.read_point_cloud(off_file_name); + + if (!check) { + std::cout << "Incorrect OFF file." << std::endl; + } else { + SC.set_type("Nerve"); + + SC.set_color_from_coordinate(coord); + SC.set_function_from_coordinate(coord); + + SC.set_graph_from_OFF(); + SC.set_resolution_with_interval_number(resolution); + SC.set_gain(gain); + SC.set_cover_from_function(); + + SC.find_simplices(); + + SC.write_info(); + + Gudhi::Simplex_tree<> stree; + SC.create_complex(stree); + SC.compute_PD(); + + // ---------------------------------------------------------------------------- + // Display information about the graph induced complex + // ---------------------------------------------------------------------------- + + if (verb) { + std::cout << "Nerve is of dimension " << stree.dimension() << " - " << stree.num_simplices() << " simplices - " + << stree.num_vertices() << " vertices." << std::endl; + + std::cout << "Iterator on Nerve simplices" << std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << std::endl; + } + } + } + + return 0; +} diff --git a/src/Nerve_GIC/utilities/Nerve.txt b/src/Nerve_GIC/utilities/Nerve.txt new file mode 100644 index 00000000..839ff45e --- /dev/null +++ b/src/Nerve_GIC/utilities/Nerve.txt @@ -0,0 +1,63 @@ +Min function value = -0.979672 and Max function value = 0.816414 +Interval 0 = [-0.979672, -0.761576] +Interval 1 = [-0.838551, -0.581967] +Interval 2 = [-0.658942, -0.402359] +Interval 3 = [-0.479334, -0.22275] +Interval 4 = [-0.299725, -0.0431415] +Interval 5 = [-0.120117, 0.136467] +Interval 6 = [0.059492, 0.316076] +Interval 7 = [0.239101, 0.495684] +Interval 8 = [0.418709, 0.675293] +Interval 9 = [0.598318, 0.816414] +Computing preimages... +Computing connected components... +.txt generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox. +5 interval(s) in dimension 0: + [-0.909111, 0.00817529] + [-0.171433, 0.367392] + [-0.171433, 0.367392] + [-0.909111, 0.745853] +0 interval(s) in dimension 1: +Nerve is of dimension 1 - 41 simplices - 21 vertices. +Iterator on Nerve simplices +1 +0 +4 +4 0 +2 +2 1 +8 +8 2 +5 +5 4 +9 +9 8 +13 +13 5 +14 +14 9 +19 +19 13 +25 +32 +20 +32 20 +33 +33 25 +26 +26 14 +26 19 +42 +42 26 +34 +34 33 +27 +27 20 +35 +35 27 +35 34 +42 35 +44 +44 35 +54 +54 44
\ No newline at end of file diff --git a/src/Nerve_GIC/utilities/VoronoiGIC.cpp b/src/Nerve_GIC/utilities/VoronoiGIC.cpp new file mode 100644 index 00000000..0182c948 --- /dev/null +++ b/src/Nerve_GIC/utilities/VoronoiGIC.cpp @@ -0,0 +1,78 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Mathieu Carrière + * + * Copyright (C) 2017 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/GIC.h> + +#include <string> +#include <vector> + +void usage(int nbArgs, char *const progName) { + std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; + std::cerr << "Usage: " << progName << " filename.off N [-v] \n"; + std::cerr << " i.e.: " << progName << " ../../data/points/human.off 100 -v \n"; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if ((argc != 3) && (argc != 4)) usage(argc, argv[0]); + + using Point = std::vector<float>; + + std::string off_file_name(argv[1]); + int m = atoi(argv[2]); + bool verb = 0; + if (argc == 4) verb = 1; + + // ---------------------------------------------------------------------------- + // Init of a graph induced complex from an OFF file + // ---------------------------------------------------------------------------- + + Gudhi::cover_complex::Cover_complex<Point> GIC; + GIC.set_verbose(verb); + + bool check = GIC.read_point_cloud(off_file_name); + + if (!check) { + std::cout << "Incorrect OFF file." << std::endl; + } else { + GIC.set_type("GIC"); + + GIC.set_color_from_coordinate(); + + GIC.set_graph_from_OFF(); + GIC.set_cover_from_Voronoi(Gudhi::Euclidean_distance(), m); + + GIC.find_simplices(); + + GIC.plot_OFF(); + + Gudhi::Simplex_tree<> stree; + GIC.create_complex(stree); + + // ---------------------------------------------------------------------------- + // Display information about the graph induced complex + // ---------------------------------------------------------------------------- + + if (verb) { + std::cout << "Graph induced complex is of dimension " << stree.dimension() << " - " << stree.num_simplices() + << " simplices - " << stree.num_vertices() << " vertices." << std::endl; + + std::cout << "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)) { + std::cout << vertex << " "; + } + std::cout << std::endl; + } + } + } + + return 0; +} diff --git a/src/Nerve_GIC/utilities/covercomplex.md b/src/Nerve_GIC/utilities/covercomplex.md new file mode 100644 index 00000000..683c1b75 --- /dev/null +++ b/src/Nerve_GIC/utilities/covercomplex.md @@ -0,0 +1,69 @@ +--- +layout: page +title: "Cover complex" +meta_title: "Cover complex" +teaser: "" +permalink: /covercomplex/ +--- +{::comment} +Leave the lines above as it is required by the web site generator 'Jekyll' +{:/comment} + + +## Nerve ## +This program builds the Nerve of a point cloud sampled on an OFF file. +The cover C comes from the preimages of intervals covering a coordinate function, +which are then refined into their connected components using the triangulation of the .OFF file. + +The program also writes a file SC.txt. +The first three lines in this file are the location of the input point cloud and the function used to compute the cover. +The fourth line contains the number of vertices nv and edges ne of the Nerve. The next nv lines represent the vertices. +Each line contains the vertex ID, the number of data points it contains, and their average color function value. +Finally, the next ne lines represent the edges, characterized by the ID of their vertices. + +**Usage** + +`Nerve <OFF input file> coordinate resolution gain [-v]` + +where + +* `coordinate` is the coordinate function to cover +* `resolution` is the number of the intervals +* `gain` is the gain for each interval +* `-v` is optional, it activates verbose mode. + +**Example** + +`Nerve ../../data/points/human.off 2 10 0.3` + +* Builds the Nerve of a point cloud sampled on a 3D human shape (human.off). +The cover C comes from the preimages of intervals (10 intervals with gain 0.3) covering the height function (coordinate 2). + +`python KeplerMapperVisuFromTxtFile.py -f ../../data/points/human.off_sc.txt` + +* Constructs `human.off_sc.html` file. You can now use your favorite web browser to visualize it. + +## VoronoiGIC ## + +This util builds the Graph Induced Complex (GIC) of a point cloud. +It subsamples *N* points in the point cloud, which act as seeds of a geodesic Voronoï diagram. +Each cell of the diagram is then an element of C. + +The program also writes a file `*_sc.off`, that is an OFF file that can be visualized with GeomView. + +**Usage** + +`VoroniGIC <OFF input file> samples_number [-v]` + +where + +* `samples_number` is the number of samples to take from the point cloud +* `-v` is optional, it activates verbose mode. + +**Example** + +`VoroniGIC ../../data/points/human.off 700` + +* Builds the Voronoi Graph Induced Complex with 700 subsamples from `human.off` file. +`../../data/points/human_sc.off` can be visualized with GeomView. + diff --git a/src/Nerve_GIC/utilities/km.py b/src/Nerve_GIC/utilities/km.py new file mode 100755 index 00000000..53024aab --- /dev/null +++ b/src/Nerve_GIC/utilities/km.py @@ -0,0 +1,390 @@ +from __future__ import division +import numpy as np +from collections import defaultdict +import json +import itertools +from sklearn import cluster, preprocessing, manifold +from datetime import datetime +import sys + +class KeplerMapper(object): + # With this class you can build topological networks from (high-dimensional) data. + # + # 1) Fit a projection/lens/function to a dataset and transform it. + # For instance "mean_of_row(x) for x in X" + # 2) Map this projection with overlapping intervals/hypercubes. + # Cluster the points inside the interval + # (Note: we cluster on the inverse image/original data to lessen projection loss). + # If two clusters/nodes have the same members (due to the overlap), then: + # connect these with an edge. + # 3) Visualize the network using HTML and D3.js. + # + # functions + # --------- + # fit_transform: Create a projection (lens) from a dataset + # map: Apply Mapper algorithm on this projection and build a simplicial complex + # visualize: Turns the complex dictionary into a HTML/D3.js visualization + + def __init__(self, verbose=2): + self.verbose = verbose + + self.chunk_dist = [] + self.overlap_dist = [] + self.d = [] + self.nr_cubes = 0 + self.overlap_perc = 0 + self.clusterer = False + + def fit_transform(self, X, projection="sum", scaler=preprocessing.MinMaxScaler()): + # Creates the projection/lens from X. + # + # Input: X. Input features as a numpy array. + # Output: projected_X. original data transformed to a projection (lens). + # + # parameters + # ---------- + # projection: Projection parameter is either a string, + # a scikit class with fit_transform, like manifold.TSNE(), + # or a list of dimension indices. + # scaler: if None, do no scaling, else apply scaling to the projection + # Default: Min-Max scaling + + self.scaler = scaler + self.projection = str(projection) + + # Detect if projection is a class (for scikit-learn) + #if str(type(projection))[1:6] == "class": #TODO: de-ugly-fy + # reducer = projection + # if self.verbose > 0: + # try: + # projection.set_params(**{"verbose":self.verbose}) + # except: + # pass + # print("\n..Projecting data using: \n\t%s\n"%str(projection)) + # X = reducer.fit_transform(X) + + # Detect if projection is a string (for standard functions) + if isinstance(projection, str): + if self.verbose > 0: + print("\n..Projecting data using: %s"%(projection)) + # Stats lenses + if projection == "sum": # sum of row + X = np.sum(X, axis=1).reshape((X.shape[0],1)) + if projection == "mean": # mean of row + X = np.mean(X, axis=1).reshape((X.shape[0],1)) + if projection == "median": # mean of row + X = np.median(X, axis=1).reshape((X.shape[0],1)) + if projection == "max": # max of row + X = np.max(X, axis=1).reshape((X.shape[0],1)) + if projection == "min": # min of row + X = np.min(X, axis=1).reshape((X.shape[0],1)) + if projection == "std": # std of row + X = np.std(X, axis=1).reshape((X.shape[0],1)) + + if projection == "dist_mean": # Distance of x to mean of X + X_mean = np.mean(X, axis=0) + X = np.sum(np.sqrt((X - X_mean)**2), axis=1).reshape((X.shape[0],1)) + + # Detect if projection is a list (with dimension indices) + if isinstance(projection, list): + if self.verbose > 0: + print("\n..Projecting data using: %s"%(str(projection))) + X = X[:,np.array(projection)] + + # Scaling + if scaler is not None: + if self.verbose > 0: + print("\n..Scaling with: %s\n"%str(scaler)) + X = scaler.fit_transform(X) + + return X + + def map(self, projected_X, inverse_X=None, clusterer=cluster.DBSCAN(eps=0.5,min_samples=3), nr_cubes=10, overlap_perc=0.1): + # This maps the data to a simplicial complex. Returns a dictionary with nodes and links. + # + # Input: projected_X. A Numpy array with the projection/lens. + # Output: complex. A dictionary with "nodes", "links" and "meta information" + # + # parameters + # ---------- + # projected_X projected_X. A Numpy array with the projection/lens. Required. + # inverse_X Numpy array or None. If None then the projection itself is used for clustering. + # clusterer Scikit-learn API compatible clustering algorithm. Default: DBSCAN + # nr_cubes Int. The number of intervals/hypercubes to create. + # overlap_perc Float. The percentage of overlap "between" the intervals/hypercubes. + + start = datetime.now() + + # Helper function + def cube_coordinates_all(nr_cubes, nr_dimensions): + # Helper function to get origin coordinates for our intervals/hypercubes + # Useful for looping no matter the number of cubes or dimensions + # Example: if there are 4 cubes per dimension and 3 dimensions + # return the bottom left (origin) coordinates of 64 hypercubes, + # as a sorted list of Numpy arrays + # TODO: elegance-ify... + l = [] + for x in range(nr_cubes): + l += [x] * nr_dimensions + return [np.array(list(f)) for f in sorted(set(itertools.permutations(l,nr_dimensions)))] + + nodes = defaultdict(list) + links = defaultdict(list) + complex = {} + self.nr_cubes = nr_cubes + self.clusterer = clusterer + self.overlap_perc = overlap_perc + + if self.verbose > 0: + print("Mapping on data shaped %s using dimensions\n"%(str(projected_X.shape))) + + # If inverse image is not provided, we use the projection as the inverse image (suffer projection loss) + if inverse_X is None: + inverse_X = projected_X + + # We chop up the min-max column ranges into 'nr_cubes' parts + self.chunk_dist = (np.max(projected_X, axis=0) - np.min(projected_X, axis=0))/nr_cubes + + # We calculate the overlapping windows distance + self.overlap_dist = self.overlap_perc * self.chunk_dist + + # We find our starting point + self.d = np.min(projected_X, axis=0) + + # Use a dimension index array on the projected X + # (For now this uses the entire dimensionality, but we keep for experimentation) + di = np.array([x for x in range(projected_X.shape[1])]) + + # Prefix'ing the data with ID's + ids = np.array([x for x in range(projected_X.shape[0])]) + projected_X = np.c_[ids,projected_X] + inverse_X = np.c_[ids,inverse_X] + + # Subdivide the projected data X in intervals/hypercubes with overlap + if self.verbose > 0: + total_cubes = len(cube_coordinates_all(nr_cubes,projected_X.shape[1])) + print("Creating %s hypercubes."%total_cubes) + + for i, coor in enumerate(cube_coordinates_all(nr_cubes,di.shape[0])): + # Slice the hypercube + hypercube = projected_X[ np.invert(np.any((projected_X[:,di+1] >= self.d[di] + (coor * self.chunk_dist[di])) & + (projected_X[:,di+1] < self.d[di] + (coor * self.chunk_dist[di]) + self.chunk_dist[di] + self.overlap_dist[di]) == False, axis=1 )) ] + + if self.verbose > 1: + print("There are %s points in cube_%s / %s with starting range %s"% + (hypercube.shape[0],i,total_cubes,self.d[di] + (coor * self.chunk_dist[di]))) + + # If at least one sample inside the hypercube + if hypercube.shape[0] > 0: + # Cluster the data point(s) in the cube, skipping the id-column + # Note that we apply clustering on the inverse image (original data samples) that fall inside the cube. + inverse_x = inverse_X[[int(nn) for nn in hypercube[:,0]]] + + clusterer.fit(inverse_x[:,1:]) + + if self.verbose > 1: + print("Found %s clusters in cube_%s\n"%(np.unique(clusterer.labels_[clusterer.labels_ > -1]).shape[0],i)) + + #Now for every (sample id in cube, predicted cluster label) + for a in np.c_[hypercube[:,0],clusterer.labels_]: + if a[1] != -1: #if not predicted as noise + cluster_id = str(coor[0])+"_"+str(i)+"_"+str(a[1])+"_"+str(coor)+"_"+str(self.d[di] + (coor * self.chunk_dist[di])) # TODO: de-rudimentary-ify + nodes[cluster_id].append( int(a[0]) ) # Append the member id's as integers + else: + if self.verbose > 1: + print("Cube_%s is empty.\n"%(i)) + + # Create links when clusters from different hypercubes have members with the same sample id. + candidates = itertools.combinations(nodes.keys(),2) + for candidate in candidates: + # if there are non-unique members in the union + if len(nodes[candidate[0]]+nodes[candidate[1]]) != len(set(nodes[candidate[0]]+nodes[candidate[1]])): + links[candidate[0]].append( candidate[1] ) + + # Reporting + if self.verbose > 0: + nr_links = 0 + for k in links: + nr_links += len(links[k]) + print("\ncreated %s edges and %s nodes in %s."%(nr_links,len(nodes),str(datetime.now()-start))) + + complex["nodes"] = nodes + complex["links"] = links + complex["meta"] = self.projection + + return complex + + def visualize(self, complex, color_function="", path_html="mapper_visualization_output.html", title="My Data", + graph_link_distance=30, graph_gravity=0.1, graph_charge=-120, custom_tooltips=None, width_html=0, + height_html=0, show_tooltips=True, show_title=True, show_meta=True, res=0,gain=0,minimum=0,maximum=0): + # Turns the dictionary 'complex' in a html file with d3.js + # + # Input: complex. Dictionary (output from calling .map()) + # Output: a HTML page saved as a file in 'path_html'. + # + # parameters + # ---------- + # color_function string. Not fully implemented. Default: "" (distance to origin) + # path_html file path as string. Where to save the HTML page. + # title string. HTML page document title and first heading. + # graph_link_distance int. Edge length. + # graph_gravity float. "Gravity" to center of layout. + # graph_charge int. charge between nodes. + # custom_tooltips None or Numpy Array. You could use "y"-label array for this. + # width_html int. Width of canvas. Default: 0 (full width) + # height_html int. Height of canvas. Default: 0 (full height) + # show_tooltips bool. default:True + # show_title bool. default:True + # show_meta bool. default:True + + # Format JSON for D3 graph + json_s = {} + json_s["nodes"] = [] + json_s["links"] = [] + k2e = {} # a key to incremental int dict, used for id's when linking + + for e, k in enumerate(complex["nodes"]): + # Tooltip and node color formatting, TODO: de-mess-ify + if custom_tooltips is not None: + tooltip_s = "<h2>Cluster %s</h2>"%k + " ".join(str(custom_tooltips[complex["nodes"][k][0]]).split(" ")) + if maximum == minimum: + tooltip_i = 0 + else: + tooltip_i = int(30*(custom_tooltips[complex["nodes"][k][0]]-minimum)/(maximum-minimum)) + json_s["nodes"].append({"name": str(k), "tooltip": tooltip_s, "group": 2 * int(np.log(complex["nodes"][k][2])), "color": tooltip_i}) + else: + tooltip_s = "<h2>Cluster %s</h2>Contains %s members."%(k,len(complex["nodes"][k])) + json_s["nodes"].append({"name": str(k), "tooltip": tooltip_s, "group": 2 * int(np.log(len(complex["nodes"][k]))), "color": str(k.split("_")[0])}) + k2e[k] = e + for k in complex["links"]: + for link in complex["links"][k]: + json_s["links"].append({"source": k2e[k], "target":k2e[link],"value":1}) + + # Width and height of graph in HTML output + if width_html == 0: + width_css = "100%" + width_js = 'document.getElementById("holder").offsetWidth-20' + else: + width_css = "%spx" % width_html + width_js = "%s" % width_html + if height_html == 0: + height_css = "100%" + height_js = 'document.getElementById("holder").offsetHeight-20' + else: + height_css = "%spx" % height_html + height_js = "%s" % height_html + + # Whether to show certain UI elements or not + if show_tooltips == False: + tooltips_display = "display: none;" + else: + tooltips_display = "" + + if show_meta == False: + meta_display = "display: none;" + else: + meta_display = "" + + if show_title == False: + title_display = "display: none;" + else: + title_display = "" + + with open(path_html,"wb") as outfile: + html = """<!DOCTYPE html> + <meta charset="utf-8"> + <meta name="generator" content="KeplerMapper"> + <title>%s | KeplerMapper</title> + <link href='https://fonts.googleapis.com/css?family=Roboto:700,300' rel='stylesheet' type='text/css'> + <style> + * {margin: 0; padding: 0;} + html { height: 100%%;} + body {background: #111; height: 100%%; font: 100 16px Roboto, Sans-serif;} + .link { stroke: #999; stroke-opacity: .333; } + .divs div { border-radius: 50%%; background: red; position: absolute; } + .divs { position: absolute; top: 0; left: 0; } + #holder { position: relative; width: %s; height: %s; background: #111; display: block;} + h1 { %s padding: 20px; color: #fafafa; text-shadow: 0px 1px #000,0px -1px #000; position: absolute; font: 300 30px Roboto, Sans-serif;} + h2 { text-shadow: 0px 1px #000,0px -1px #000; font: 700 16px Roboto, Sans-serif;} + .meta { position: absolute; opacity: 0.9; width: 220px; top: 80px; left: 20px; display: block; %s background: #000; line-height: 25px; color: #fafafa; border: 20px solid #000; font: 100 16px Roboto, Sans-serif;} + div.tooltip { position: absolute; width: 380px; display: block; %s padding: 20px; background: #000; border: 0px; border-radius: 3px; pointer-events: none; z-index: 999; color: #FAFAFA;} + } + </style> + <body> + <div id="holder"> + <h1>%s</h1> + <p class="meta"> + <b>Lens</b><br>%s<br><br> + <b>Length of intervals</b><br>%s<br><br> + <b>Overlap percentage</b><br>%s%%<br><br> + <b>Color Function</b><br>%s + </p> + </div> + <script src="https://cdnjs.cloudflare.com/ajax/libs/d3/3.5.5/d3.min.js"></script> + <script> + var width = %s, + height = %s; + var color = d3.scale.ordinal() + .domain(["0","1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13","14","15","16","17","18","19","20","21","22","23","24","25","26","27","28","29","30"]) + .range(["#FF0000","#FF1400","#FF2800","#FF3c00","#FF5000","#FF6400","#FF7800","#FF8c00","#FFa000","#FFb400","#FFc800","#FFdc00","#FFf000","#fdff00","#b0ff00","#65ff00","#17ff00","#00ff36","#00ff83","#00ffd0","#00e4ff","#00c4ff","#00a4ff","#00a4ff","#0084ff","#0064ff","#0044ff","#0022ff","#0002ff","#0100ff","#0300ff","#0500ff"]); + var force = d3.layout.force() + .charge(%s) + .linkDistance(%s) + .gravity(%s) + .size([width, height]); + var svg = d3.select("#holder").append("svg") + .attr("width", width) + .attr("height", height); + + var div = d3.select("#holder").append("div") + .attr("class", "tooltip") + .style("opacity", 0.0); + + var divs = d3.select('#holder').append('div') + .attr('class', 'divs') + .attr('style', function(d) { return 'overflow: hidden; width: ' + width + 'px; height: ' + height + 'px;'; }); + + graph = %s; + force + .nodes(graph.nodes) + .links(graph.links) + .start(); + var link = svg.selectAll(".link") + .data(graph.links) + .enter().append("line") + .attr("class", "link") + .style("stroke-width", function(d) { return Math.sqrt(d.value); }); + var node = divs.selectAll('div') + .data(graph.nodes) + .enter().append('div') + .on("mouseover", function(d) { + div.transition() + .duration(200) + .style("opacity", .9); + div .html(d.tooltip + "<br/>") + .style("left", (d3.event.pageX + 100) + "px") + .style("top", (d3.event.pageY - 28) + "px"); + }) + .on("mouseout", function(d) { + div.transition() + .duration(500) + .style("opacity", 0); + }) + .call(force.drag); + + node.append("title") + .text(function(d) { return d.name; }); + force.on("tick", function() { + link.attr("x1", function(d) { return d.source.x; }) + .attr("y1", function(d) { return d.source.y; }) + .attr("x2", function(d) { return d.target.x; }) + .attr("y2", function(d) { return d.target.y; }); + node.attr("cx", function(d) { return d.x; }) + .attr("cy", function(d) { return d.y; }) + .attr('style', function(d) { return 'width: ' + (d.group * 2) + 'px; height: ' + (d.group * 2) + 'px; ' + 'left: '+(d.x-(d.group))+'px; ' + 'top: '+(d.y-(d.group))+'px; background: '+color(d.color)+'; box-shadow: 0px 0px 3px #111; box-shadow: 0px 0px 33px '+color(d.color)+', inset 0px 0px 5px rgba(0, 0, 0, 0.2);'}) + ; + }); + </script>"""%(title,width_css, height_css, title_display, meta_display, tooltips_display, title,complex["meta"],res,gain*100,color_function,width_js,height_js,graph_charge,graph_link_distance,graph_gravity,json.dumps(json_s)) + outfile.write(html.encode("utf-8")) + if self.verbose > 0: + print("\nWrote d3.js graph to '%s'"%path_html) diff --git a/src/Nerve_GIC/utilities/km.py.COPYRIGHT b/src/Nerve_GIC/utilities/km.py.COPYRIGHT new file mode 100644 index 00000000..bef7b121 --- /dev/null +++ b/src/Nerve_GIC/utilities/km.py.COPYRIGHT @@ -0,0 +1,26 @@ +km.py is a fork of https://github.com/MLWave/kepler-mapper. +Only the visualization part has been kept (Mapper part has been removed). + +This file has te following Copyright : + +The MIT License (MIT) + +Copyright (c) 2015 Triskelion - HJ van Veen - info@mlwave.com + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. |
