diff options
Diffstat (limited to 'src/Witness_complex/include/gudhi/Witness_complex.h')
-rw-r--r-- | src/Witness_complex/include/gudhi/Witness_complex.h | 312 |
1 files changed, 44 insertions, 268 deletions
diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h index 8316fe3e..b218611b 100644 --- a/src/Witness_complex/include/gudhi/Witness_complex.h +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -27,9 +27,7 @@ #include <boost/iterator/transform_iterator.hpp> #include <algorithm> #include <utility> -#include "gudhi/reader_utils.h" #include "gudhi/distance_functions.h" -#include "gudhi/Simplex_tree.h" #include <vector> #include <list> #include <set> @@ -47,47 +45,32 @@ namespace Gudhi { - /** Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: - * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ - * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that - * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. + /** + \class Witness_complex + \brief Constructs the witness complex for the given set of witnesses and landmarks. + \ingroup witness_complex */ - template<typename FiltrationValue = double, - typename SimplexKey = int, - typename VertexHandle = int> - class Witness_complex: public Simplex_tree<> { - + template< class Simplicial_complex> + class Witness_complex { + private: struct Active_witness { int witness_id; int landmark_id; - Simplex_handle simplex_handle; - Active_witness(int witness_id_, int landmark_id_, Simplex_handle simplex_handle_) + Active_witness(int witness_id_, int landmark_id_) : witness_id(witness_id_), - landmark_id(landmark_id_), - simplex_handle(simplex_handle_) + landmark_id(landmark_id_) {} }; - - - public: - - - /** \brief Type for the vertex handle. - * - * Must be a signed integer type. It admits a total order <. */ - typedef VertexHandle Vertex_handle; - - /* Type of node in the simplex tree. */ - typedef Simplex_tree_node_explicit_storage<Simplex_tree> Node; - /* Type of dictionary Vertex_handle -> Node for traversing the simplex tree. */ - typedef typename boost::container::flat_map<Vertex_handle, Node> Dictionary; - typedef typename Dictionary::iterator Simplex_handle; + private: + typedef typename Simplicial_complex::Simplex_handle Simplex_handle; + typedef typename Simplicial_complex::Vertex_handle Vertex_handle; + typedef std::vector< double > Point_t; typedef std::vector< Point_t > Point_Vector; @@ -102,33 +85,29 @@ namespace Gudhi { private: int nbL; // Number of landmarks double density; // Desired density - + Simplicial_complex& sc; // Simplicial complex + public: - /** \brief Set number of landmarks to nbL_ + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** @name Constructor */ - void setNbL(int nbL_) - { - nbL = nbL_; - } - /** \brief Set density to density_ - */ - void setDensity(double density_) - { - density = density_; - } + //@{ /** - * /brief Iterative construction of the witness complex basing on a matrix of k nearest neighbours of the form {witnesses}x{landmarks}. - * Landmarks are supposed to be in [0,nbL-1] - */ - + * \brief Iterative construction of the witness complex. + * \details The witness complex is written in sc_ basing on a matrix knn of + * nearest neighbours of the form {witnesses}x{landmarks}. + * Parameter dim serves as the limit for the number of closest landmarks to consider. + * Landmarks are supposed to be in [0,nbL_-1] + */ template< typename KNearestNeighbours > - void witness_complex(KNearestNeighbours & knn) - //void witness_complex(std::vector< std::vector< Vertex_handle > > & knn) + Witness_complex(KNearestNeighbours & knn, + Simplicial_complex & sc_, + int nbL_, + int dim ): nbL(nbL_), sc(sc_) { - std::cout << "**Start the procedure witness_complex" << std::endl; //Construction of the active witness list int nbW = knn.size(); typeVectorVertex vv; @@ -144,16 +123,14 @@ namespace Gudhi { // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore counter++; vv = {i}; - returnValue = insert_simplex(vv, Filtration_value(0.0)); + returnValue = sc.insert_simplex(vv, Filtration_value(0.0)); /* TODO Error if not inserted : normally no need here though*/ } int k=1; /* current dimension in iterative construction */ for (int i=0; i != nbW; ++i) active_w.push_back(i); - std::cout << "k=0, active witnesses: " << active_w.size() << std::endl; //std::cout << "Successfully added edges" << std::endl; - int D = knn[0].size(); - while (!active_w.empty() && k < D ) + while (!active_w.empty() && k < dim ) { //std::cout << "Started the step k=" << k << std::endl; typename ActiveWitnessList::iterator it = active_w.begin(); @@ -166,84 +143,20 @@ namespace Gudhi { { for (int i = 0; i != k+1; ++i) simplex_vector.push_back(knn[*it][i]); - returnValue = insert_simplex(simplex_vector,0.0); + returnValue = sc.insert_simplex(simplex_vector,0.0); it++; } else active_w.erase(it++); //First increase the iterator and then erase the previous element } - std::cout << "k=" << k << ", active witnesses: " << active_w.size() << std::endl; k++; } } - -private: - /** \brief Print functions - */ - void print_sc(Siblings * sibl) - { - if (sibl == NULL) - std::cout << "&"; - else - print_children(sibl->members_); - } + //@} - void print_children(Dictionary map) - { - std::cout << "("; - if (!map.empty()) - { - std::cout << map.begin()->first; - if (has_children(map.begin())) - print_sc(map.begin()->second.children()); - typename Dictionary::iterator it; - for (it = map.begin()+1; it != map.end(); ++it) - { - std::cout << "," << it->first; - if (has_children(it)) - print_sc(it->second.children()); - } - } - std::cout << ")"; - } - - public: - /** \brief Print functions - */ - - void st_to_file(std::ofstream& out_file) - { - sc_to_file(out_file, root()); - } - private: - void sc_to_file(std::ofstream& out_file, Siblings * sibl) - { - assert(sibl); - children_to_file(out_file, sibl->members_); - } - void children_to_file(std::ofstream& out_file, Dictionary& map) - { - out_file << "(" << std::flush; - if (!map.empty()) - { - out_file << map.begin()->first << std::flush; - if (has_children(map.begin())) - sc_to_file(out_file, map.begin()->second.children()); - typename Dictionary::iterator it; - for (it = map.begin()+1; it != map.end(); ++it) - { - out_file << "," << it->first << std::flush; - if (has_children(it)) - sc_to_file(out_file, it->second.children()); - } - } - out_file << ")" << std::flush; - } - - /** \brief Check if the facets of the k-dimensional simplex witnessed * by witness witness_id are already in the complex. * inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id @@ -252,8 +165,8 @@ private: bool all_faces_in(KNearestNeighbours &knn, int witness_id, int k) { //std::cout << "All face in with the landmark " << inserted_vertex << std::endl; - std::vector< VertexHandle > facet; - //VertexHandle curr_vh = curr_sh->first; + std::vector< Vertex_handle > facet; + //Vertex_handle curr_vh = curr_sh->first; // CHECK ALL THE FACETS for (int i = 0; i != k+1; ++i) { @@ -265,14 +178,14 @@ private: facet.push_back(knn[witness_id][j]); } }//endfor - if (find(facet) == null_simplex()) + if (sc.find(facet) == sc.null_simplex()) return false; //std::cout << "++++ finished loop safely\n"; } //endfor return true; } - - template <typename T> + + template <typename T> void print_vector(std::vector<T> v) { std::cout << "["; @@ -288,162 +201,25 @@ private: std::cout << "]"; } - template <typename T> - void print_vvector(std::vector< std::vector <T> > vv) - { - std::cout << "["; - if (!vv.empty()) - { - print_vector(*(vv.begin())); - for (auto it = vv.begin()+1; it != vv.end(); ++it) - { - std::cout << ","; - print_vector(*it); - } - } - std::cout << "]\n"; - } - public: -/** - * \brief Landmark choice strategy by iteratively adding the landmark the furthest from the - * current landmark set - * \arg W is the vector of points which will be the witnesses - * \arg nbP is the number of witnesses - * \arg nbL is the number of landmarks - * \arg WL is the matrix of the nearest landmarks with respect to witnesses (output) - */ - - template <typename KNearestNeighbours> - void landmark_choice_by_furthest_points(Point_Vector &W, int nbP, KNearestNeighbours &WL) - { - Point_Vector wit_land_dist(nbP,std::vector<double>()); // distance matrix witness x landmarks - typeVectorVertex chosen_landmarks; // landmark list - - WL = KNearestNeighbours(nbP,std::vector<int>()); - int current_number_of_landmarks=0; // counter for landmarks - double curr_max_dist = 0; // used for defining the furhest point from L - double curr_dist; // used to stock the distance from the current point to L - double infty = std::numeric_limits<double>::infinity(); // infinity (see next entry) - std::vector< double > dist_to_L(nbP,infty); // vector of current distances to L from points - int curr_max_w=0; // the point currently furthest from L - int j; - int temp_swap_int; - double temp_swap_double; - - //CHOICE OF THE FIRST LANDMARK - std::cout << "Enter the first landmark stage\n"; - srand(354698); - int rand_int = rand()% nbP; - curr_max_w = rand_int; //For testing purposes a pseudo-random number is used here - - for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) - { - //curr_max_w at this point is the next landmark - chosen_landmarks.push_back(curr_max_w); - for (auto v: WL) - v.push_back(current_number_of_landmarks); - for (int i = 0; i < nbP; ++i) - { - curr_dist = euclidean_distance(W[i],W[chosen_landmarks[current_number_of_landmarks]]); - wit_land_dist[i].push_back(curr_dist); - WL[i].push_back(current_number_of_landmarks); - if (curr_dist < dist_to_L[i]) - dist_to_L[i] = curr_dist; - j = current_number_of_landmarks; - while (j > 0 && wit_land_dist[i][j-1] > wit_land_dist[i][j]) - { - temp_swap_int = WL[i][j]; - WL[i][j] = WL[i][j-1]; - WL[i][j-1] = temp_swap_int; - temp_swap_double = wit_land_dist[i][j]; - wit_land_dist[i][j] = wit_land_dist[i][j-1]; - wit_land_dist[i][j-1] = temp_swap_double; - --j; - } - } - curr_max_dist = 0; - for (int i = 0; i < nbP; ++i) { - if (dist_to_L[i] > curr_max_dist) - { - curr_max_dist = dist_to_L[i]; - curr_max_w = i; - } - } - } - } - - /** \brief Landmark choice strategy by taking random vertices for landmarks. - * - */ - - // template <typename KNearestNeighbours> - void landmark_choice_by_random_points(Point_Vector &W, int nbP, std::set<int> &L) - { - int current_number_of_landmarks=0; // counter for landmarks - - srand(24660); - int chosen_landmark = rand()%nbP; - for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) - { - while (L.find(chosen_landmark) != L.end()) - { - srand((int)clock()); - chosen_landmark = rand()% nbP; - } - L.insert(chosen_landmark); - } - } - - - /** \brief Construct the matrix |W|x(D+1) of D+1 closest landmarks - * where W is the set of witnesses and D is the ambient dimension - */ - template <typename KNearestNeighbours> - void nearest_landmarks(Point_Vector &W, std::set<int> &L, KNearestNeighbours &WL) - { - int D = W[0].size(); - int nbP = W.size(); - WL = KNearestNeighbours(nbP,std::vector<int>()); - typedef std::pair<double,int> dist_i; - typedef bool (*comp)(dist_i,dist_i); - for (int W_i = 0; W_i < nbP; W_i++) - { - std::priority_queue<dist_i, std::vector<dist_i>, comp> l_heap([&](dist_i j1, dist_i j2){return j1.first > j2.first;}); - std::set<int>::iterator L_it; - int L_i; - for (L_it = L.begin(), L_i=0; L_it != L.end(); L_it++, L_i++) - { - dist_i dist = std::make_pair(euclidean_distance(W[W_i],W[*L_it]), L_i); - l_heap.push(dist); - } - for (int i = 0; i < D+1; i++) - { - dist_i dist = l_heap.top(); - WL[W_i].push_back(dist.second); - l_heap.pop(); - } - } - } - - - public: - /** \brief Verification if every simplex in the complex is witnessed + /** + * \brief Verification if every simplex in the complex is witnessed. + * \remark Added for debugging purposes. */ template< class KNearestNeighbors > - bool is_witness_complex(KNearestNeighbors WL) + bool is_witness_complex(KNearestNeighbors & knn) { //bool final_result = true; - for (Simplex_handle sh: complex_simplex_range()) + for (Simplex_handle sh: sc.complex_simplex_range()) { bool is_witnessed = false; typeVectorVertex simplex; int nbV = 0; //number of verticed in the simplex - for (int v: simplex_vertex_range(sh)) + for (int v: sc.simplex_vertex_range(sh)) simplex.push_back(v); nbV = simplex.size(); - for (typeVectorVertex w: WL) + for (typeVectorVertex w: knn) { bool has_vertices = true; for (int v: simplex) |