From bc9d9bdaec5e4ca1604cca2a5251031f25c64b0f Mon Sep 17 00:00:00 2001 From: skachano Date: Wed, 20 May 2015 15:16:13 +0000 Subject: Added perturbation witness complex construction git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/witness@590 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 5c478d127e0fefa8696be5ca72ab89efa2e2dfc6 --- src/Witness_complex/example/CMakeLists.txt | 29 +- .../example/witness_complex_perturbations.cpp | 377 +++++++++++++++++++++ 2 files changed, 404 insertions(+), 2 deletions(-) create mode 100644 src/Witness_complex/example/witness_complex_perturbations.cpp (limited to 'src/Witness_complex/example') diff --git a/src/Witness_complex/example/CMakeLists.txt b/src/Witness_complex/example/CMakeLists.txt index 315ffc36..a3e7b3f5 100644 --- a/src/Witness_complex/example/CMakeLists.txt +++ b/src/Witness_complex/example/CMakeLists.txt @@ -56,12 +56,37 @@ if(CGAL_FOUND) add_executable ( witness_complex_knn_landmarks witness_complex_knn_landmarks.cpp ) target_link_libraries(witness_complex_knn_landmarks ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) add_test(witness_complex_knn_landmarks ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_knn_landmarks ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) - + #add_executable ( witness_complex_perturbations witness_complex_perturbations.cpp ) + #target_link_libraries(witness_complex_perturbations ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + #add_test(witness_complex_perturbations ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_perturbations ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) else() - message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.") + message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.") endif() else() message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile Alpha shapes feature. Version 4.6.0 is required.") endif () endif() +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.6.0) + message(STATUS "CGAL version: ${CGAL_VERSION}.") + include( ${CGAL_USE_FILE} ) + + find_package(Eigen3 3.1.0) + if (EIGEN3_FOUND) + message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.") + include( ${EIGEN3_USE_FILE} ) + include_directories (BEFORE "../../include") + add_executable ( witness_complex_perturbations witness_complex_perturbations.cpp ) + target_link_libraries(witness_complex_perturbations ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + add_test(witness_complex_perturbations ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_perturbations ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) + add_executable ( range_queries range_queries.cpp ) + target_link_libraries(range_queries ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + add_test(range_queries ${CMAKE_CURRENT_BINARY_DIR}/range_queries ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) + else() + message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.") + endif() + else() + message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile Alpha shapes feature. Version 4.6.0 is required.") + endif () +endif() diff --git a/src/Witness_complex/example/witness_complex_perturbations.cpp b/src/Witness_complex/example/witness_complex_perturbations.cpp new file mode 100644 index 00000000..afb42807 --- /dev/null +++ b/src/Witness_complex/example/witness_complex_perturbations.cpp @@ -0,0 +1,377 @@ +/* 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 . + */ + +#include +#include +#include +#include +#include +#include + +#include +#include +//#include + +//#include "gudhi/graph_simplicial_complex.h" +#include "gudhi/Witness_complex.h" +#include "gudhi/reader_utils.h" +//#include + +//#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +using namespace Gudhi; +//using namespace boost::filesystem; + +typedef std::vector< Vertex_handle > typeVectorVertex; + +//typedef std::pair typeSimplex; +//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; + +typedef CGAL::Epick_d K; +typedef K::FT FT; +typedef K::Point_d Point_d; +typedef CGAL::Search_traits< + FT, Point_d, + typename K::Cartesian_const_iterator_d, + typename K::Construct_cartesian_const_iterator_d> Traits_base; +typedef CGAL::Search_traits_adapter< + std::ptrdiff_t, Point_d*, Traits_base> STraits; +//typedef K TreeTraits; +typedef CGAL::Orthogonal_k_neighbor_search K_neighbor_search; +typedef K_neighbor_search::Tree Tree; +typedef K_neighbor_search::Distance Distance; +typedef K_neighbor_search::iterator KNS_iterator; +typedef K_neighbor_search::iterator KNS_range; +typedef boost::container::flat_map Point_etiquette_map; +typedef CGAL::Kd_tree Tree2; + +typedef CGAL::Fuzzy_sphere Fuzzy_sphere; + +typedef std::vector Point_Vector; + +typedef CGAL::Euclidean_distance Euclidean_distance; +//typedef K::Equal_d Equal_d; +typedef CGAL::Random_points_in_ball_d Random_point_iterator; +/** + * \brief Customized version of read_points + * which takes into account a possible nbP first line + * + */ +inline void +read_points_cust ( std::string file_name , Point_Vector & points) +{ + std::ifstream in_file (file_name.c_str(),std::ios::in); + if(!in_file.is_open()) + { + std::cerr << "Unable to open file " << file_name << std::endl; + return; + } + std::string line; + double x; + while( getline ( in_file , line ) ) + { + std::vector< double > point; + std::istringstream iss( line ); + while(iss >> x) { point.push_back(x); } + Point_d p(point.begin(), point.end()); + if (point.size() != 1) + points.push_back(p); + } + in_file.close(); +} + +/* +void read_points_to_tree (std::string file_name, Tree& tree) +{ + //I assume here that tree is empty + std::ifstream in_file (file_name.c_str(),std::ios::in); + if(!in_file.is_open()) + { + std::cerr << "Unable to open file " << file_name << std::endl; + return; + } + std::string line; + double x; + while( getline ( in_file , line ) ) + { + std::vector coords; + std::istringstream iss( line ); + while(iss >> x) { coords.push_back(x); } + if (coords.size() != 1) + { + Point_d point(coords.begin(), coords.end()); + tree.insert(point); + } + } + in_file.close(); +} +*/ + +void write_wl( std::string file_name, std::vector< std::vector > & WL) +{ + std::ofstream ofs (file_name, std::ofstream::out); + for (auto w : WL) + { + for (auto l: w) + ofs << l << " "; + ofs << "\n"; + } + ofs.close(); +} + + + +/** Function that chooses landmarks from W and place it in the kd-tree L. + * Note: nbL hould be removed if the code moves to Witness_complex + */ +void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, std::vector& landmarks_ind) +{ + std::cout << "Enter landmark choice to kd tree\n"; + //std::vector landmarks; + int chosen_landmark; + //std::pair res = std::make_pair(L_i.begin(),false); + Point_d* p; + srand(24660); + for (int i = 0; i < nbL; i++) + { + // while (!res.second) + // { + chosen_landmark = rand()%nbP; + p = &W[chosen_landmark]; + //L_i.emplace(chosen_landmark,i); + // } + landmarks.push_back(*p); + landmarks_ind.push_back(chosen_landmark); + //std::cout << "Added landmark " << chosen_landmark << std::endl; + } + } + + +int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector& landmarks_ind) +{ + //******************** Constructing a WL matrix + int nbP = W.size(); + int nbL = landmarks.size(); + //Point_Vector landmarks_ = landmarks; + Euclidean_distance ed; + //Equal_d ed; + FT lambda = ed.transformed_distance(landmarks[0],landmarks[1]); + //FT lambda = 0.1;//Euclidean_distance(); + std::vector< std::vector > WL(nbP); + Tree L(boost::counting_iterator(0), + boost::counting_iterator(nbL), + typename Tree::Splitter(), + STraits(&(landmarks[0]))); + /*Tree2 L2(boost::counting_iterator(0), + boost::counting_iterator(nbL), + typename Tree::Splitter(), + STraits(&(landmarks[0]))); + */ + std::cout << "Enter (D+1) nearest landmarks\n"; + //std::cout << "Size of the tree is " << L.size() << std::endl; + int D = W[0].size(); + for (int i = 0; i < nbP; i++) + { + //std::cout << "Entered witness number " << i << std::endl; + Point_d& w = W[i]; + //std::cout << "Safely constructed a point\n"; + ////Search D+1 nearest neighbours from the tree of landmarks L + K_neighbor_search search(L, w, D+1, FT(0), true, + CGAL::Distance_adapter>(&(landmarks[0])) ); + //std::cout << "Safely found nearest landmarks\n"; + for(K_neighbor_search::iterator it = search.begin(); it != search.end(); ++it) + { + //std::cout << "Entered KNN_it with point at distance " << it->second << "\n"; + //Point_etiquette_map::iterator itm = L_i.find(it->first); + //assert(itm != L_i.end()); + //std::cout << "Entered KNN_it with point at distance " << it->second << "\n"; + WL[i].push_back(it->first); + //std::cout << "ITFIRST " << it->first << std::endl; + //std::cout << i << " " << it->first << ": " << it->second << std::endl; + } + if (i == landmarks_ind[WL[i][0]]) + { + //std::cout << "'"; + FT dist = ed.transformed_distance(W[i], landmarks[WL[i][1]]); + if (dist < lambda) + lambda = dist; + } + } + //std::cout << "\n"; + /* + std::string out_file = "wl_result"; + write_wl(out_file,WL); + */ + //******************** Constructng a witness complex + std::cout << "Entered witness complex construction\n"; + Witness_complex<> witnessComplex; + witnessComplex.setNbL(nbL); + witnessComplex.witness_complex(WL); + //******************** Making a set of bad link landmarks + std::cout << "Entered bad links\n"; + std::set< int > perturbL; + int count_badlinks = 0; + std::cout << "Bad links around "; + for (auto u: witnessComplex.complex_vertex_range()) + if (!witnessComplex.has_good_link(u)) + { + //std::cout << "Landmark " << u << " start!" << std::endl; + //perturbL.insert(u); + count_badlinks++; + std::cout << u << " "; + Point_d& l = landmarks[u]; + Fuzzy_sphere fs(l, sqrt(lambda)*2, 0, STraits(&(landmarks[0]))); + L.search(std::insert_iterator>(perturbL,perturbL.begin()),fs); + //L.search(std::inserter(perturbL,perturbL.begin()),fs); + //L.search(std::ostream_iterator(std::cout,"\n"),fs); + //std::cout << "PerturbL size is " << perturbL.size() << std::endl; + } + std::cout << "\nBad links total: " << count_badlinks << " Points to perturb: " << perturbL.size() << std::endl; + //std::cout << "landmark[0][0] before" << landmarks[0][0] << std::endl; + //*********************** Perturb bad link landmarks + + for (auto u: perturbL) + { + Random_point_iterator rp(D,sqrt(lambda)/2); + //std::cout << landmarks[u] << std::endl; + + std::vector point; + for (int i = 0; i < D; i++) + { + point.push_back(W[landmarks_ind[u]][i] + (*rp)[i]); + } + landmarks[u] = Point_d(point); + //std::cout << landmarks[u] << std::endl; + } + + //std::cout << "landmark[0][0] after" << landmarks[0][0] << std::endl; + std::cout << "lambda=" << lambda << std::endl; + // Write the WL matrix in a file + /* + mkdir("output", S_IRWXU); + const size_t last_slash_idx = file_name.find_last_of("/"); + if (std::string::npos != last_slash_idx) + { + file_name.erase(0, last_slash_idx + 1); + } + + char buffer[100]; + int i = sprintf(buffer,"stree_result.txt"); + + if (i >= 0) + { + std::string out_file = (std::string)buffer; + std::ofstream ofs (out_file, std::ofstream::out); + witnessComplex.st_to_file(ofs); + ofs.close(); + } + */ + return count_badlinks; +} + + +int main (int argc, char * const argv[]) +{ + if (argc != 3) + { + std::cerr << "Usage: " << argv[0] + << " path_to_point_file nbL \n"; + return 0; + } + /* + boost::filesystem::path p; + + for (; argc > 2; --argc, ++argv) + p /= argv[1]; + */ + std::string file_name = argv[1]; + int nbL = atoi(argv[2]); + + //clock_t start, end; + //Construct the Simplex Tree + //Witness_complex<> witnessComplex; + + std::cout << "Let the carnage begin!\n"; + Point_Vector point_vector; + read_points_cust(file_name, point_vector); + //std::cout << "Successfully read the points\n"; + //witnessComplex.setNbL(nbL); + // witnessComplex.witness_complex_from_points(point_vector); + int nbP = point_vector.size(); + //std::vector > WL(nbP); + //std::set L; + Point_Vector L; + std::vector chosen_landmarks; + //Point_etiquette_map L_i; + //start = clock(); + //witnessComplex.landmark_choice_by_furthest_points(point_vector, point_vector.size(), WL); + landmark_choice(point_vector, nbP, nbL, L, chosen_landmarks); + int bl = 1; + for (int i = 0; bl != 0; i++) + { + std::cout << "========== Start iteration " << i << " ========\n"; + bl = landmark_perturbation(point_vector, L, chosen_landmarks); + } + //end = clock(); + + /* + std::cout << "Landmark choice took " + << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; + start = clock(); + witnessComplex.witness_complex(WL); + // + end = clock(); + std::cout << "Howdy world! The process took " + << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; + */ + + /* + out_file = "output/"+file_name+"_"+argv[2]+".stree"; + std::ofstream ofs (out_file, std::ofstream::out); + witnessComplex.st_to_file(ofs); + ofs.close(); + + out_file = "output/"+file_name+"_"+argv[2]+".badlinks"; + std::ofstream ofs2(out_file, std::ofstream::out); + witnessComplex.write_bad_links(ofs2); + ofs2.close(); + */ +} -- cgit v1.2.3