From 8f54c437e0b895368c6151584811ce7df1575ea0 Mon Sep 17 00:00:00 2001 From: skachano Date: Wed, 9 Dec 2015 13:03:03 +0000 Subject: Modified Witness_complex.h + more doc git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/witness@936 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 771e0a9cbfd0035738aee1a07b5f7c1f56bf13fc --- .../example/simple_witness_complex.cpp | 38 +- .../gudhi/Landmark_choice_by_furthest_point.h | 98 +++ .../gudhi/Landmark_choice_by_random_point.h | 82 ++ .../include/gudhi/Relaxed_witness_complex.h | 886 --------------------- .../include/gudhi/Witness_complex.h | 312 +------- .../include/gudhi/Witness_complex_doc.h | 37 + 6 files changed, 283 insertions(+), 1170 deletions(-) create mode 100644 src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h create mode 100644 src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h delete mode 100644 src/Witness_complex/include/gudhi/Relaxed_witness_complex.h create mode 100644 src/Witness_complex/include/gudhi/Witness_complex_doc.h (limited to 'src/Witness_complex') diff --git a/src/Witness_complex/example/simple_witness_complex.cpp b/src/Witness_complex/example/simple_witness_complex.cpp index e95f67a8..6731f135 100644 --- a/src/Witness_complex/example/simple_witness_complex.cpp +++ b/src/Witness_complex/example/simple_witness_complex.cpp @@ -23,30 +23,36 @@ #include #include //#include "gudhi/graph_simplicial_complex.h" +#include "gudhi/Simplex_tree.h" #include "gudhi/Witness_complex.h" using namespace Gudhi; typedef std::vector< Vertex_handle > typeVectorVertex; +typedef Witness_complex> WitnessComplex; //typedef std::pair typeSimplex; //typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; int main (int argc, char * const argv[]) { - Witness_complex<> witnessComplex = Witness_complex<>(); - std::vector< typeVectorVertex > KNN; - typeVectorVertex witness0 = {1,0,5,2,6,3,4}; KNN.push_back(witness0 ); - typeVectorVertex witness1 = {2,6,4,5,0,1,3}; KNN.push_back(witness1 ); - typeVectorVertex witness2 = {3,4,2,1,5,6,0}; KNN.push_back(witness2 ); - typeVectorVertex witness3 = {4,2,1,3,5,6,0}; KNN.push_back(witness3 ); - typeVectorVertex witness4 = {5,1,6,0,2,3,4}; KNN.push_back(witness4 ); - typeVectorVertex witness5 = {6,0,5,2,1,3,4}; KNN.push_back(witness5 ); - typeVectorVertex witness6 = {0,5,6,1,2,3,4}; KNN.push_back(witness6 ); - typeVectorVertex witness7 = {2,6,4,5,3,1,0}; KNN.push_back(witness7 ); - typeVectorVertex witness8 = {1,2,5,4,3,6,0}; KNN.push_back(witness8 ); - typeVectorVertex witness9 = {3,4,0,6,5,1,2}; KNN.push_back(witness9 ); - typeVectorVertex witness10 = {5,0,1,3,6,2,4}; KNN.push_back(witness10); - typeVectorVertex witness11 = {5,6,1,0,2,3,4}; KNN.push_back(witness11); - typeVectorVertex witness12 = {1,6,0,5,2,3,4}; KNN.push_back(witness12); - witnessComplex.witness_complex(KNN); + Simplex_tree<> complex; + std::vector< typeVectorVertex > knn; + typeVectorVertex witness0 = {1,0,5,2,6,3,4}; knn.push_back(witness0 ); + typeVectorVertex witness1 = {2,6,4,5,0,1,3}; knn.push_back(witness1 ); + typeVectorVertex witness2 = {3,4,2,1,5,6,0}; knn.push_back(witness2 ); + typeVectorVertex witness3 = {4,2,1,3,5,6,0}; knn.push_back(witness3 ); + typeVectorVertex witness4 = {5,1,6,0,2,3,4}; knn.push_back(witness4 ); + typeVectorVertex witness5 = {6,0,5,2,1,3,4}; knn.push_back(witness5 ); + typeVectorVertex witness6 = {0,5,6,1,2,3,4}; knn.push_back(witness6 ); + typeVectorVertex witness7 = {2,6,4,5,3,1,0}; knn.push_back(witness7 ); + typeVectorVertex witness8 = {1,2,5,4,3,6,0}; knn.push_back(witness8 ); + typeVectorVertex witness9 = {3,4,0,6,5,1,2}; knn.push_back(witness9 ); + typeVectorVertex witness10 = {5,0,1,3,6,2,4}; knn.push_back(witness10); + typeVectorVertex witness11 = {5,6,1,0,2,3,4}; knn.push_back(witness11); + typeVectorVertex witness12 = {1,6,0,5,2,3,4}; knn.push_back(witness12); + WitnessComplex witnessComplex(knn, complex, 7); + if (witnessComplex.is_witness_complex(knn)) + std::cout << "Witness complex is good\n"; + else + std::cout << "Witness complex is bad\n"; } diff --git a/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h b/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h new file mode 100644 index 00000000..44bf9dbb --- /dev/null +++ b/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h @@ -0,0 +1,98 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Siargey Kachanovich + * + * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef GUDHI_LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ +#define GUDHI_LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ + +/** + * \class Landmark_choice_by_furthest_point + * \brief The class `Landmark_choice_by_furthest_point` allows to construct the matrix + * of closest landmarks per witness by iteratively choosing the furthest witness + * from the set of already chosen landmarks as the new landmark. + * \ingroup witness_complex + */ + +class Landmark_choice_by_random_point { + +/** + * \brief Landmark choice strategy by iteratively adding the furthest witness from the + * current landmark set as the new landmark. It takes a random access range `points` and + * writes {witness}*{closest landmarks} matrix in `knn`. + */ + + template + Landmark_choice_by_furthest_points(Point_random_access_range &points, + KNearestNeighbours &knn) + { + int nb_points = points.end() - points.begin(); + std::vector> wit_land_dist(nb_points, std::vector()); // distance matrix witness x landmarks + typeVectorVertex chosen_landmarks; // landmark list + + knn = KNearestNeighbours(nb_points, std::vector()); + int current_number_of_landmarks=0; // counter for landmarks + double curr_max_dist = 0; // used for defining the furhest point from L + const double infty = std::numeric_limits::infinity(); // infinity (see next entry) + std::vector< double > dist_to_L(nb_points,infty); // vector of current distances to L from points + + //CHOICE OF THE FIRST LANDMARK + int rand_int = rand() % nb_points; + int 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: knn) + v.push_back(current_number_of_landmarks); + int i = 0; + for (const auto& p: points) + { + // used to stock the distance from the current point to L + double curr_dist = euclidean_distance(p, points.begin() + chosen_landmarks[current_number_of_landmarks]); + wit_land_dist[i].push_back(curr_dist); + knn[i].push_back(current_number_of_landmarks); + if (curr_dist < dist_to_L[i]) + dist_to_L[i] = curr_dist; + int j = current_number_of_landmarks; + while (j > 0 && wit_land_dist[i][j-1] > wit_land_dist[i][j]) + { + std::swap(knn[i][j], knn[i][j-1]); + std::swap(wit_land_dist[i][j-1], wit_land_dist[i][j-1]); + --j; + } + ++i; + } + curr_max_dist = 0; + for (auto dist: dist_to_L) { + if (dist > curr_max_dist) + { + curr_max_dist = dist; + curr_max_w = i; + } + } + } + } + +}; + +#endif diff --git a/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h b/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h new file mode 100644 index 00000000..bc3e72d9 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h @@ -0,0 +1,82 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Siargey Kachanovich + * + * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef GUDHI_LANDMARK_CHOICE_BY_RANDOM_POINT_H_ +#define GUDHI_LANDMARK_CHOICE_BY_RANDOM_POINT_H_ + +/** + * \class Landmark_choice_by_random_point + * \brief The class `Landmark_choice_by_random_point` allows to construct the matrix + * of closest landmarks per witness by iteratively choosing a random non-chosen witness + * as a new landmark. + * \ingroup witness_complex + */ + +class Landmark_choice_by_random_point { + + + /** \brief Landmark choice strategy by taking random vertices for landmarks. + * It takes a random access range points and outputs a matrix {witness}*{closest landmarks} + * in knn. + */ + + template + void landmark_choice_by_random_points(Point_random_access_range &points, KNearestNeighbours &knn) + { + int nbP = points.end() - points.begin(); + std::set &landmarks; + int current_number_of_landmarks=0; // counter for landmarks + + int chosen_landmark = rand()%nbP; + for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) + { + while (landmarks.find(chosen_landmark) != landmarks.end()) + chosen_landmark = rand()% nbP; + landmarks.insert(chosen_landmark); + } + + int D = points.begin().size(); + typedef std::pair dist_i; + typedef bool (*comp)(dist_i,dist_i); + for (int points_i = 0; points_i < nbP; points_i++) + { + std::priority_queue, comp> l_heap([&](dist_i j1, dist_i j2){return j1.first > j2.first;}); + std::set::iterator landmarks_it; + int landmarks_i = 0; + for (landmarks_it = landmarks.begin(), landmarks_i=0; landmarks_it != landmarks.end(); landmarks_it++, landmarks_i++) + { + dist_i dist = std::make_pair(euclidean_distance(points[points_i],points[*landmarks_it]), landmarks_i); + l_heap.push(dist); + } + for (int i = 0; i < D+1; i++) + { + dist_i dist = l_heap.top(); + knn[points_i].push_back(dist.second); + l_heap.pop(); + } + } + } + +}; + +#endif diff --git a/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h b/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h deleted file mode 100644 index c869628f..00000000 --- a/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h +++ /dev/null @@ -1,886 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -#ifndef GUDHI_RELAXED_WITNESS_COMPLEX_H_ -#define GUDHI_RELAXED_WITNESS_COMPLEX_H_ - -#include -#include -#include -#include -#include "gudhi/reader_utils.h" -#include "gudhi/distance_functions.h" -#include "gudhi/Simplex_tree.h" -#include -#include -#include -#include -#include -#include -#include -#include - -// Needed for nearest neighbours -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -// Needed for the adjacency graph in bad link search -#include -#include -#include - -namespace Gudhi { - - - /** \addtogroup simplex_tree - * 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. - */ - template - class Witness_complex: public Simplex_tree<> { - - 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_) - : witness_id(witness_id_), - landmark_id(landmark_id_), - simplex_handle(simplex_handle_) - {} - }; - - - - - 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 Node; - /* Type of dictionary Vertex_handle -> Node for traversing the simplex tree. */ - typedef typename boost::container::flat_map Dictionary; - typedef typename Dictionary::iterator Simplex_handle; - - typedef std::vector< double > Point_t; - typedef std::vector< Point_t > Point_Vector; - - typedef std::vector< Vertex_handle > typeVectorVertex; - typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; - typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; - - typedef int Witness_id; - typedef int Landmark_id; - typedef std::list< Vertex_handle > ActiveWitnessList; - - private: - /** Number of landmarks - */ - int nbL; - /** Desired density - */ - double density; - - public: - - /** \brief Set number of landmarks to nbL_ - */ - void setNbL(int nbL_) - { - nbL = nbL_; - } - - /** \brief Set density to density_ - */ - void setDensity(double density_) - { - density = density_; - } - - /** - * /brief Iterative construction of the relaxed witness complex basing on a matrix of k nearest neighbours of the form {witnesses}x{landmarks} and (1+epsilon)-limit table {witnesses}*{landmarks} consisting of iterators of k nearest neighbor matrix. - * The line lengths can differ, however both matrices have the same corresponding line lengths. - */ - - template< typename KNearestNeighbours, typename OPELimits > - void relaxed_witness_complex(KNearestNeighbours & knn, OPELimits & rl) - //void witness_complex(std::vector< std::vector< Vertex_handle > > & knn) - { - std::cout << "**Start the procedure witness_complex" << std::endl; - //Construction of the active witness list - int nbW = knn.size(); - //int nbL = knn.at(0).size(); - typeVectorVertex vv; - //typeSimplex simplex; - //typePairSimplexBool returnValue; - //int counter = 0; - /* The list of still useful witnesses - * it will diminuish in the course of iterations - */ - ActiveWitnessList active_w;// = new ActiveWitnessList(); - for (int i=0; i != nbL; ++i) { - // initial fill of 0-dimensional simplices - // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore - //counter++; - vv = {i}; - 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 */ - //std::cout << "Successfully added landmarks" << std::endl; - // PRINT2 - //print_sc(root()); std::cout << std::endl; - for (int i=0; i != nbW; ++i) - active_w.push_back(i); - /* - int u,v; // two extremities of an edge - if (nbL > 1) // if the supposed dimension of the complex is >0 - { - for (int i=0; i != nbW; ++i) - { - // initial fill of active witnesses list - u = knn[i][0]; - v = knn[i][1]; - vv = {u,v}; - this->insert_simplex(vv,Filtration_value(0.0)); - //print_sc(root()); std::cout << std::endl; - //std::cout << "Added edges" << std::endl; - } - //print_sc(root()); - - } - */ - std::cout << "k=0, active witnesses: " << active_w.size() << std::endl; - //std::cout << "Successfully added edges" << std::endl; - //count_good = {0,0}; - //count_bad = {0,0}; - while (!active_w.empty() && k < nbL ) - { - //count_good.push_back(0); - //count_bad.push_back(0); - //std::cout << "Started the step k=" << k << std::endl; - typename ActiveWitnessList::iterator aw_it = active_w.begin(); - while (aw_it != active_w.end()) - { - std::vector simplex; - bool ok = add_all_faces_of_dimension(k, knn[*aw_it].begin(), rl[*aw_it].begin(), simplex, knn[*aw_it].end(), knn[*aw_it].end()); - if (!ok) - active_w.erase(aw_it++); //First increase the iterator and then erase the previous element - else - aw_it++; - } - std::cout << "k=" << k << ", active witnesses: " << active_w.size() << std::endl; - k++; - } - //print_sc(root()); std::cout << std::endl; - } - - /* \brief Adds recursively all the faces of a certain dimension dim witnessed by the same witness - * Iterator is needed to know until how far we can take landmarks to form simplexes - * simplex is the prefix of the simplexes to insert - * The output value indicates if the witness rests active or not - */ - bool add_all_faces_of_dimension(int dim, std::vector::iterator curr_l, typename std::vector< std::vector::iterator >::iterator curr_until, std::vector& simplex, std::vector::iterator until, std::vector::iterator end) - { - /* - std::ofstream ofs ("stree_result.txt", std::ofstream::out); - st_to_file(ofs); - ofs.close(); - */ - //print_sc(root()); - bool will_be_active = false; - if (dim > 0) - for (std::vector::iterator it = curr_l; it != until && it != end; ++it, ++curr_until) - { - simplex.push_back(*it); - if (find(simplex) != null_simplex()) - will_be_active = will_be_active || add_all_faces_of_dimension(dim-1, it+1, curr_until+1, simplex, until, end); - simplex.pop_back(); - if (until == end) - until = *curr_until; - } - else if (dim == 0) - for (std::vector::iterator it = curr_l; it != until && it != end; ++it, ++curr_until) - { - simplex.push_back(*it); - if (all_faces_in(simplex)) - { - will_be_active = true; - insert_simplex(simplex, 0.0); - } - simplex.pop_back(); - if (until == end) - until = *curr_until; - } - return will_be_active; - } - - /** \brief Construction of witness complex from points given explicitly - * nbL must be set to the right value of landmarks for strategies - * FURTHEST_POINT_STRATEGY and RANDOM_POINT_STRATEGY and - * density must be set to the right value for DENSITY_STRATEGY - */ - // void witness_complex_from_points(Point_Vector point_vector) - // { - // std::vector > WL; - // landmark_choice_by_random_points(point_vector, point_vector.size(), WL); - // witness_complex(WL); - // } - -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 - */ - bool all_faces_in(std::vector& simplex) - { - //std::cout << "All face in with the landmark " << inserted_vertex << std::endl; - std::vector< VertexHandle > facet; - //VertexHandle curr_vh = curr_sh->first; - // CHECK ALL THE FACETS - for (std::vector::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) - { - facet.clear(); - //facet = {}; - for (std::vector::iterator it = simplex.begin(); it != simplex.end(); ++it) - if (it != not_it) - facet.push_back(*it); - if (find(facet) == null_simplex()) - return false; - } //endfor - return true; - } - - template - void print_vector(std::vector v) - { - std::cout << "["; - if (!v.empty()) - { - std::cout << *(v.begin()); - for (auto it = v.begin()+1; it != v.end(); ++it) - { - std::cout << ","; - std::cout << *it; - } - } - std::cout << "]"; - } - - template - void print_vvector(std::vector< std::vector > 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 - void landmark_choice_by_furthest_points(Point_Vector &W, int nbP, KNearestNeighbours &WL) - { - //std::cout << "Enter landmark_choice_by_furthest_points "<< std::endl; - //std::cout << "W="; print_vvector(W); - //double density = 5.; - Point_Vector wit_land_dist(nbP,std::vector()); // distance matrix witness x landmarks - typeVectorVertex chosen_landmarks; // landmark list - - WL = KNearestNeighbours(nbP,std::vector()); - 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::infinity(); // infinity (see next entry) - std::vector< double > dist_to_L(nbP,infty); // vector of current distances to L from points - // double mindist = infty; - 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); - //std::cout << "**********Entered loop with current number of landmarks = " << current_number_of_landmarks << std::endl; - //std::cout << "WL="; print_vvector(WL); - //std::cout << "WLD="; print_vvector(wit_land_dist); - //std::cout << "landmarks="; print_vector(chosen_landmarks); std::cout << std::endl; - for (auto v: WL) - v.push_back(current_number_of_landmarks); - for (int i = 0; i < nbP; ++i) - { - // iteration on points in W. update of distance vectors - - //std::cout << "In the loop with i=" << i << " and landmark=" << chosen_landmarks[current_number_of_landmarks] << std::endl; - //std::cout << "W[i]="; print_vector(W[i]); std::cout << " W[landmark]="; print_vector(W[chosen_landmarks[current_number_of_landmarks]]); std::cout << std::endl; - curr_dist = euclidean_distance(W[i],W[chosen_landmarks[current_number_of_landmarks]]); - //std::cout << "The problem is not in distance function\n"; - wit_land_dist[i].push_back(curr_dist); - WL[i].push_back(current_number_of_landmarks); - //std::cout << "Push't back\n"; - if (curr_dist < dist_to_L[i]) - dist_to_L[i] = curr_dist; - j = current_number_of_landmarks; - //std::cout << "First half complete\n"; - while (j > 0 && wit_land_dist[i][j-1] > wit_land_dist[i][j]) - { - // sort the closest landmark vector for every witness - 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; - } - //std::cout << "result WL="; print_vvector(WL); - //std::cout << "result WLD="; print_vvector(wit_land_dist); - //std::cout << "result distL="; print_vector(dist_to_L); std::cout << std::endl; - //std::cout << "End loop\n"; - } - //std::cout << "Distance to landmarks="; print_vector(dist_to_L); std::cout << std::endl; - 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; - } - } - //std::cout << "Chose " << curr_max_w << " as new landmark\n"; - } - //std::cout << endl; - } - - /** \brief Landmark choice strategy by taking random vertices for landmarks. - * - */ - - // template - // void landmark_choice_by_random_points(Point_Vector &W, int nbP, KNearestNeighbours &WL) - // { - // std::cout << "Enter landmark_choice_by_random_points "<< std::endl; - // //std::cout << "W="; print_vvector(W); - // std::unordered_set< int > chosen_landmarks; // landmark set - - // Point_Vector wit_land_dist(nbP,std::vector()); // distance matrix witness x landmarks - - // WL = KNearestNeighbours(nbP,std::vector()); - // int current_number_of_landmarks=0; // counter for landmarks - - // srand(24660); - // int chosen_landmark = rand()%nbP; - // double curr_dist; - - // //int j; - // //int temp_swap_int; - // //double temp_swap_double; - - - // for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) - // { - // while (chosen_landmarks.find(chosen_landmark) != chosen_landmarks.end()) - // { - // srand((int)clock()); - // chosen_landmark = rand()% nbP; - // //std::cout << chosen_landmark << "\n"; - // } - // chosen_landmarks.insert(chosen_landmark); - // //std::cout << "**********Entered loop with current number of landmarks = " << current_number_of_landmarks << std::endl; - // //std::cout << "WL="; print_vvector(WL); - // //std::cout << "WLD="; print_vvector(wit_land_dist); - // //std::cout << "landmarks="; print_vector(chosen_landmarks); std::cout << std::endl; - // for (auto v: WL) - // v.push_back(current_number_of_landmarks); - // for (int i = 0; i < nbP; ++i) - // { - // // iteration on points in W. update of distance vectors - - // //std::cout << "In the loop with i=" << i << " and landmark=" << chosen_landmarks[current_number_of_landmarks] << std::endl; - // //std::cout << "W[i]="; print_vector(W[i]); std::cout << " W[landmark]="; print_vector(W[chosen_landmarks[current_number_of_landmarks]]); std::cout << std::endl; - // curr_dist = euclidean_distance(W[i],W[chosen_landmark]); - // //std::cout << "The problem is not in distance function\n"; - // wit_land_dist[i].push_back(curr_dist); - // WL[i].push_back(current_number_of_landmarks); - // //std::cout << "Push't back\n"; - // //j = current_number_of_landmarks; - // //std::cout << "First half complete\n"; - // //std::cout << "result WL="; print_vvector(WL); - // //std::cout << "result WLD="; print_vvector(wit_land_dist); - // //std::cout << "End loop\n"; - // } - // } - // for (int i = 0; i < nbP; i++) - // { - // sort(WL[i].begin(), WL[i].end(), [&](int j1, int j2){return wit_land_dist[i][j1] < wit_land_dist[i][j2];}); - // } - // //std::cout << endl; - // } - - /** \brief Landmark choice strategy by taking random vertices for landmarks. - * - */ - - // template - void landmark_choice_by_random_points(Point_Vector &W, int nbP, std::set &L) - { - std::cout << "Enter landmark_choice_by_random_points "<< std::endl; - //std::cout << "W="; print_vvector(W); - //std::unordered_set< int > chosen_landmarks; // landmark set - - //Point_Vector wit_land_dist(nbP,std::vector()); // distance matrix witness x landmarks - - //WL = KNearestNeighbours(nbP,std::vector()); - int current_number_of_landmarks=0; // counter for landmarks - - srand(24660); - int chosen_landmark = rand()%nbP; - //double curr_dist; - //int j; - //int temp_swap_int; - //double temp_swap_double; - 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; - //std::cout << chosen_landmark << "\n"; - } - L.insert(chosen_landmark); - //std::cout << "**********Entered loop with current number of landmarks = " << current_number_of_landmarks << std::endl; - //std::cout << "WL="; print_vvector(WL); - //std::cout << "WLD="; print_vvector(wit_land_dist); - //std::cout << "landmarks="; print_vector(chosen_landmarks); std::cout << std::endl; - // for (auto v: WL) - // v.push_back(current_number_of_landmarks); - // for (int i = 0; i < nbP; ++i) - // { - // // iteration on points in W. update of distance vectors - - // //std::cout << "In the loop with i=" << i << " and landmark=" << chosen_landmarks[current_number_of_landmarks] << std::endl; - // //std::cout << "W[i]="; print_vector(W[i]); std::cout << " W[landmark]="; print_vector(W[chosen_landmarks[current_number_of_landmarks]]); std::cout << std::endl; - // curr_dist = euclidean_distance(W[i],W[chosen_landmark]); - // //std::cout << "The problem is not in distance function\n"; - // wit_land_dist[i].push_back(curr_dist); - // WL[i].push_back(current_number_of_landmarks); - // //std::cout << "Push't back\n"; - // //j = current_number_of_landmarks; - // //std::cout << "First half complete\n"; - // //std::cout << "result WL="; print_vvector(WL); - // //std::cout << "result WLD="; print_vvector(wit_land_dist); - // //std::cout << "End loop\n"; - // } - } - // for (int i = 0; i < nbP; i++) - // { - // sort(WL[i].begin(), WL[i].end(), [&](int j1, int j2){return wit_land_dist[i][j1] < wit_land_dist[i][j2];}); - // } - //std::cout << endl; - } - - - /** \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 - void nearest_landmarks(Point_Vector &W, std::set &L, KNearestNeighbours &WL) - { - int D = W[0].size(); - int nbP = W.size(); - WL = KNearestNeighbours(nbP,std::vector()); - typedef std::pair dist_i; - typedef bool (*comp)(dist_i,dist_i); - for (int W_i = 0; W_i < nbP; W_i++) - { - //std::cout << "<<<<<<<<<<<<<<" << W_i <<"\n"; - std::priority_queue, comp> l_heap([&](dist_i j1, dist_i j2){return j1.first > j2.first;}); - std::set::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); - //WL[W_i].insert(WL[W_i].begin(),dist.second); - //std::cout << dist.first << " " << dist.second << std::endl; - l_heap.pop(); - } - } - } - - /** \brief Search and output links around vertices that are not pseudomanifolds - * - */ - void write_bad_links(std::ofstream& out_file) - { - out_file << "Bad links list\n"; - std::cout << "Entered write_bad_links\n"; - //typeVectorVertex testv = {9,15,17}; - //int count = 0; - for (auto v: complex_vertex_range()) - { - //std::cout << "Vertex " << v << ":\n"; - std::vector< Vertex_handle > link_vertices; - // Fill link_vertices - for (auto u: complex_vertex_range()) - { - typeVectorVertex edge = {u,v}; - if (u != v && find(edge) != null_simplex()) - link_vertices.push_back(u); - } - /* - print_vector(link_vertices); - std::cout << "\n"; - */ - // Find the dimension - typeVectorVertex empty_simplex = {}; - int d = link_dim(link_vertices, link_vertices.begin(),-1, empty_simplex); - //std::cout << " dim " << d << "\n"; - //Siblings* curr_sibl = root(); - if (link_is_pseudomanifold(link_vertices,d)) - count_good[d]++; - //out_file << "Bad link at " << v << "\n"; - } - //out_file << "Number of bad links: " << count << "/" << root()->size(); - //std::cout << "Number of bad links: " << count << "/" << root()->size() << std::endl; - nc = nbL; - for (unsigned int i = 0; i != count_good.size(); i++) - { - out_file << "count_good[" << i << "] = " << count_good[i] << std::endl; - nc -= count_good[i]; - if (count_good[i] != 0) - std::cout << "count_good[" << i << "] = " << count_good[i] << std::endl; - } - for (unsigned int i = 0; i != count_bad.size(); i++) - { - out_file << "count_bad[" << i << "] = " << count_bad[i] << std::endl; - nc -= count_bad[i]; - if (count_bad[i] != 0) - std::cout << "count_bad[" << i << "] = " << count_bad[i] << std::endl; - } - std::cout << "not_connected = " << nc << std::endl; - } - - private: - - std::vector count_good; - std::vector count_bad; - int nc; - - int link_dim(std::vector< Vertex_handle >& link_vertices, - typename std::vector< Vertex_handle >::iterator curr_v, - int curr_d, - typeVectorVertex& curr_simplex) - { - //std::cout << "Entered link_dim for " << *(curr_v-1) << "\n"; - Simplex_handle sh; - int final_d = curr_d; - typename std::vector< Vertex_handle >::iterator it; - for (it = curr_v; it != link_vertices.end(); ++it) - { - curr_simplex.push_back(*it); - /* - std::cout << "Searching for "; - print_vector(curr_simplex); - std::cout << " curr_dim " << curr_d << " final_dim " << final_d; - */ - sh = find(curr_simplex); - if (sh != null_simplex()) - { - //std::cout << " -> " << *it << "\n"; - int d = link_dim(link_vertices, it+1, curr_d+1, curr_simplex); - if (d > final_d) - final_d = d; - } - /* - else - std::cout << "\n"; - */ - curr_simplex.pop_back(); - } - return final_d; - } - - // color is false is a (d-1)-dim face, true is a d-dim face - //typedef bool Color; - // graph is an adjacency list - typedef typename boost::adjacency_list Adj_graph; - // map that gives to a certain simplex its node in graph and its dimension - //typedef std::pair Reference; - typedef boost::graph_traits::vertex_descriptor Vertex_t; - typedef boost::graph_traits::edge_descriptor Edge_t; - - typedef boost::container::flat_map Graph_map; - - /* \brief Verifies if the simplices formed by vertices given by link_vertices - * form a pseudomanifold. - * The idea is to make a bipartite graph, where vertices are the d- and (d-1)-dimensional - * faces and edges represent adjacency between them. - */ - bool link_is_pseudomanifold(std::vector< Vertex_handle >& link_vertices, - int dimension) - { - Adj_graph adj_graph; - Graph_map d_map, f_map; // d_map = map for d-dimensional simplices - // f_map = map for its facets - typeVectorVertex empty_vector = {}; - add_vertices(link_vertices, - link_vertices.begin(), - adj_graph, - d_map, - f_map, - empty_vector, - 0, dimension); - //std::cout << "DMAP_SIZE: " << d_map.size() << "\n"; - //std::cout << "FMAP_SIZE: " << f_map.size() << "\n"; - add_edges(adj_graph, d_map, f_map); - for (auto f_map_it : f_map) - { - //std::cout << "Degree of " << f_map_it.first->first << " is " << boost::out_degree(f_map_it.second, adj_graph) << "\n"; - if (boost::out_degree(f_map_it.second, adj_graph) != 2) - { - count_bad[dimension]++; - return false; - } - } - // At this point I know that all (d-1)-simplices are adjacent to exactly 2 d-simplices - // What is left is to check the connexity - std::vector components(boost::num_vertices(adj_graph)); - return (boost::connected_components(adj_graph, &components[0]) == 1); - } - - void add_vertices(typeVectorVertex& link_vertices, - typename typeVectorVertex::iterator curr_v, - Adj_graph& adj_graph, - Graph_map& d_map, - Graph_map& f_map, - typeVectorVertex& curr_simplex, - int curr_d, - int dimension) - { - Simplex_handle sh; - Vertex_t vert; - typename typeVectorVertex::iterator it; - std::pair resPair; - //typename Graph_map::iterator resPair; - //Add vertices - //std::cout << "Entered add vertices\n"; - for (it = curr_v; it != link_vertices.end(); ++it) - { - curr_simplex.push_back(*it); - /* - std::cout << "Searching for "; - print_vector(curr_simplex); - std::cout << " curr_dim " << curr_d << " d " << dimension << ""; - */ - sh = find(curr_simplex); - if (sh != null_simplex()) - { - //std::cout << " added\n"; - if (curr_d == dimension) - { - vert = boost::add_vertex(adj_graph); - resPair = d_map.emplace(sh,vert); - } - else - { - if (curr_d == dimension-1) - { - vert = boost::add_vertex(adj_graph); - resPair = f_map.emplace(sh,vert); - } - add_vertices(link_vertices, - it+1, - adj_graph, - d_map, - f_map, - curr_simplex, - curr_d+1, dimension); - } - } - /* - else - std::cout << "\n"; - */ - curr_simplex.pop_back(); - } - } - - void add_edges(Adj_graph& adj_graph, - Graph_map& d_map, - Graph_map& f_map) - { - Simplex_handle sh; - // Add edges - //std::cout << "Entered add edges:\n"; - typename Graph_map::iterator map_it; - for (auto d_map_pair : d_map) - { - //std::cout << "*"; - sh = d_map_pair.first; - Vertex_t d_vert = d_map_pair.second; - for (auto facet_sh : boundary_simplex_range(sh)) - //for (auto f_map_it : f_map) - { - //std::cout << "'"; - map_it = f_map.find(facet_sh); - //We must have all the facets in the graph at this point - assert(map_it != f_map.end()); - Vertex_t f_vert = map_it->second; - //std::cout << "Added edge " << sh->first << "-" << map_it->first->first << "\n"; - boost::add_edge(d_vert,f_vert,adj_graph); - } - } - } - - ////////////////////////////////////////////////////////////////////////////////////////////////// - //***********COLLAPSES**************************************************************************// - ////////////////////////////////////////////////////////////////////////////////////////////////// - - - - - - - -}; //class Witness_complex - - - -} // namespace Guhdi - -#endif 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 #include #include -#include "gudhi/reader_utils.h" #include "gudhi/distance_functions.h" -#include "gudhi/Simplex_tree.h" #include #include #include @@ -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 - 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 Node; - /* Type of dictionary Vertex_handle -> Node for traversing the simplex tree. */ - typedef typename boost::container::flat_map 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 + + template void print_vector(std::vector v) { std::cout << "["; @@ -288,162 +201,25 @@ private: std::cout << "]"; } - template - void print_vvector(std::vector< std::vector > 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 - void landmark_choice_by_furthest_points(Point_Vector &W, int nbP, KNearestNeighbours &WL) - { - Point_Vector wit_land_dist(nbP,std::vector()); // distance matrix witness x landmarks - typeVectorVertex chosen_landmarks; // landmark list - - WL = KNearestNeighbours(nbP,std::vector()); - 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::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 - void landmark_choice_by_random_points(Point_Vector &W, int nbP, std::set &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 - void nearest_landmarks(Point_Vector &W, std::set &L, KNearestNeighbours &WL) - { - int D = W[0].size(); - int nbP = W.size(); - WL = KNearestNeighbours(nbP,std::vector()); - typedef std::pair dist_i; - typedef bool (*comp)(dist_i,dist_i); - for (int W_i = 0; W_i < nbP; W_i++) - { - std::priority_queue, comp> l_heap([&](dist_i j1, dist_i j2){return j1.first > j2.first;}); - std::set::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) diff --git a/src/Witness_complex/include/gudhi/Witness_complex_doc.h b/src/Witness_complex/include/gudhi/Witness_complex_doc.h new file mode 100644 index 00000000..22cfe992 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Witness_complex_doc.h @@ -0,0 +1,37 @@ +#ifndef WITNESS_COMPLEX_DOC_ +#define WITNESS_COMPLEX_DOC_ + +/** + \defgroup witness_complex Witness complex + + \author Siargey Kachanovich + + \section Definitions + + Witness complex \f$ Wit(W,L) \f$ is a simplicial complex defined on two sets of points in \f$\mathbb{R}^D\f$: + + \li \f$W\f$ set of **witnesses** and + \li \f$L \subseteq W\f$ set of **landmarks**. + + The simplices are based on landmarks + and a simplex belongs to the witness complex if and only if it is witnessed, that is: + + \f$ \sigma \subset L \f$ is witnessed if there exists a point \f$w \in W\f$ such that + w is closer to the vertices of \f$ \sigma \f$ than other points in \f$ L \f$ and all of its faces are witnessed as well. + + \section Implementation + + Two classes are implemented in this module: Gudhi::Witness_complex and Gudhi::Relaxed_witness_complex. + + While Gudhi::Witness_complex represents the classical witness complex, Gudhi::Relaxed_witness_complex takes an additional positive real parameter \f$ \alpha \f$ and constructs simplices \f$ \sigma \f$, for which + there exists \f$ w \in W \f$, such that \f$ d(p,w) < d(q,w) + \alpha \f$ for all \f$ p \in \sigma, q \in L\setminus \sigma \f$. + + In both cases, the constructors take a {witness}x{closest_landmarks} table, + which can be constructed by two additional classes Landmark_choice_by_furthest_point and Landmark_choice_by_random_point also included in the module. + + \copyright GNU General Public License v3. + + + */ + +#endif -- cgit v1.2.3