Copied Clément's example for 2 functions

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/subsampling_and_spatialsearching@1338 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: 42f29038057980965bfead2c73cbc013448493fb

5 files changed, 266 insertions, 26 deletions

diff --git a/src/Subsampling/example/CMakeLists.txt b/src/Subsampling/example/CMakeLists.txt
index e7a8a9f7..e1e7cc71 100644
--- a/src/Subsampling/example/CMakeLists.txt
+++ b/src/Subsampling/example/CMakeLists.txt
@@ -11,9 +11,9 @@ if(CGAL_FOUND)
       include( ${EIGEN3_USE_FILE} )
       include_directories (BEFORE "../../include")
 
-      #add_executable( Subsampling_example_pick_random_points example_pick_random_points.cpp )
+      add_executable( Subsampling_example_pick_random_points example_pick_random_points.cpp )
 
-      #add_executable( Subsampling_example_choose_farthest_point example_choose_farthest_point.cpp )
+      add_executable( Subsampling_example_choose_by_farthest_point example_choose_by_farthest_point.cpp )
       #target_link_libraries(Subsampling_example_choose_farthest_point ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
       
       add_executable(Subsampling_example_sparsify_point_set example_sparsify_point_set.cpp)
diff --git a/src/Subsampling/example/example_choose_by_farthest_point.cpp b/src/Subsampling/example/example_choose_by_farthest_point.cpp
new file mode 100644
index 00000000..5b81bc9d
--- /dev/null
+++ b/src/Subsampling/example/example_choose_by_farthest_point.cpp
@@ -0,0 +1,52 @@
+/*    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) 2016 INRIA
+ *
+ *    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 <http://www.gnu.org/licenses/>.
+ */
+
+
+#include <gudhi/choose_by_farthest_point.h>
+
+#include <CGAL/Epick_d.h>
+#include <CGAL/Random.h>
+
+#include <array>
+#include <vector>
+#include <iterator>
+
+int main (void)
+{
+  typedef CGAL::Epick_d<CGAL::Dimension_tag<4> >   K;
+  typedef typename K::FT                           FT;
+  typedef typename K::Point_d                      Point_d;
+  
+  CGAL::Random rd;
+
+  std::vector<Point_d> points;
+  for (int i = 0 ; i < 500 ; ++i)
+    points.push_back(Point_d(std::array<FT,4>({rd.get_double(-1.,1),rd.get_double(-1.,1),rd.get_double(-1.,1),rd.get_double(-1.,1)})));
+
+  K k;
+  std::vector<Point_d> results;
+  Gudhi::subsampling::choose_by_farthest_point(k, points, 100, std::back_inserter(results));
+  std::cout << "Before sparsification: " << points.size() << " points.\n";
+  std::cout << "After  sparsification: " << results.size() << " points.\n";
+
+  return 0;
+}
diff --git a/src/Subsampling/example/example_pick_random_points.cpp b/src/Subsampling/example/example_pick_random_points.cpp
new file mode 100644
index 00000000..49a027a4
--- /dev/null
+++ b/src/Subsampling/example/example_pick_random_points.cpp
@@ -0,0 +1,52 @@
+/*    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) 2016 INRIA
+ *
+ *    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 <http://www.gnu.org/licenses/>.
+ */
+
+
+#include <gudhi/pick_random_points.h>
+
+#include <CGAL/Epick_d.h>
+#include <CGAL/Random.h>
+
+#include <array>
+#include <vector>
+#include <iterator>
+
+int main (void)
+{
+  typedef CGAL::Epick_d<CGAL::Dimension_tag<4> >   K;
+  typedef typename K::FT                           FT;
+  typedef typename K::Point_d                      Point_d;
+  
+  CGAL::Random rd;
+
+  std::vector<Point_d> points;
+  for (int i = 0 ; i < 500 ; ++i)
+    points.push_back(Point_d(std::array<FT,4>({rd.get_double(-1.,1),rd.get_double(-1.,1),rd.get_double(-1.,1),rd.get_double(-1.,1)})));
+
+  K k;
+  std::vector<Point_d> results;
+  Gudhi::subsampling::pick_random_points(points, 100, std::back_inserter(results));
+  std::cout << "Before sparsification: " << points.size() << " points.\n";
+  std::cout << "After  sparsification: " << results.size() << " points.\n";
+
+  return 0;
+}
diff --git a/src/Subsampling/include/gudhi/choose_by_farthest_point.h b/src/Subsampling/include/gudhi/choose_by_farthest_point.h
index 52647c16..d1db0d1a 100644
--- a/src/Subsampling/include/gudhi/choose_by_farthest_point.h
+++ b/src/Subsampling/include/gudhi/choose_by_farthest_point.h
@@ -25,6 +25,14 @@
 
 #include <boost/range.hpp>
 
+#include <gudhi/Spatial_tree_data_structure.h>
+
+#include <gudhi/Clock.h>
+
+#include <CGAL/Search_traits.h>
+#include <CGAL/Search_traits_adapter.h>
+#include <CGAL/Fuzzy_sphere.h>
+
 #include <iterator>
 #include <algorithm>  // for sort
 #include <vector>
@@ -45,10 +53,11 @@ namespace subsampling {
   template < typename Kernel,
              typename Point_container,
              typename OutputIterator>
-  void choose_by_farthest_point( Kernel& k,
-                                 Point_container const &points,
-                                 int final_size,
-                                 OutputIterator output_it)
+  void choose_by_farthest_point_old( Kernel& k,
+                                     Point_container const &points,
+                                     int final_size,
+                                     int starting_point,
+                                     OutputIterator output_it)
   {
     typename Kernel::Squared_distance_d sqdist = k.squared_distance_d_object();
     
@@ -60,11 +69,7 @@ namespace subsampling {
     const double infty = std::numeric_limits<double>::infinity();  // infinity (see next entry)
     std::vector< double > dist_to_L(nb_points, infty);  // vector of current distances to L from points
 
-    // Choose randomly the first landmark 
-    std::random_device rd;
-    std::mt19937 gen(rd());
-    std::uniform_int_distribution<> dis(1, 6);
-    int curr_max_w = dis(gen);
+    int curr_max_w = starting_point;
 
     for (current_number_of_landmarks = 0; current_number_of_landmarks != final_size; current_number_of_landmarks++) {
       // curr_max_w at this point is the next landmark
@@ -86,7 +91,103 @@ namespace subsampling {
         }
     }
   }
+  
+  template < typename Kernel,
+             typename Point_container,
+             typename OutputIterator>
+  void choose_by_farthest_point_old( Kernel& k,
+                                     Point_container const &points,
+                                     int final_size,
+                                     OutputIterator output_it)
+  {
+    // Choose randomly the first landmark 
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_int_distribution<> dis(1, 6);
+    int starting_point = dis(gen);
+    choose_by_farthest_point_old(k, points, final_size, starting_point, output_it);
+  }
+
+  template < typename Kernel,
+             typename Point_container,
+             typename OutputIterator>
+  void choose_by_farthest_point( Kernel& k,
+                                 Point_container const &points,
+                                 int final_size,
+                                 int starting_point,
+                                 OutputIterator output_it)
+  {
+    // typedef typename Kernel::Point_d Point_d;
+    // typedef typename Kernel::FT FT;
+    // typedef CGAL::Search_traits<
+    //   FT, Point_d,
+    //   typename Kernel::Cartesian_const_iterator_d,
+    //   typename Kernel::Construct_cartesian_const_iterator_d> Traits_base;
 
+    // typedef CGAL::Search_traits_adapter< std::ptrdiff_t, Point_d*, Traits_base > STraits;
+    // typedef CGAL::Fuzzy_sphere< STraits > Fuzzy_sphere;
+    
+    typename Kernel::Squared_distance_d sqdist = k.squared_distance_d_object();
+    
+    int nb_points = boost::size(points);
+    assert(nb_points >= final_size);
+
+    Clock t;
+    Gudhi::spatial_searching::Spatial_tree_data_structure< Kernel, Point_container> tree(points);
+    t.end();
+    //std::cout << "Constructed the Kd tree: " << t.num_seconds()  << " s." << std::endl;
+  
+    //CGAL::Fuzzy_sphere< CGAL::Search_trai>
+    
+    int current_number_of_landmarks = 0;  // counter for landmarks
+    const double infty = std::numeric_limits<double>::infinity();  // infinity (see next entry)
+    double curr_max_dist = infty;  // used for defining the furhest point from L
+    std::vector< double > dist_to_L(nb_points, infty);  // vector of current distances to L from points
+    
+    // Choose randomly the first landmark 
+    int curr_max_w = starting_point;
+
+    for (current_number_of_landmarks = 0; current_number_of_landmarks != final_size; current_number_of_landmarks++) {
+      // curr_max_w at this point is the next landmark
+      *output_it++ = points[curr_max_w];
+      // std::cout << curr_max_w << "\n";
+      //for (auto& p : points) {
+      auto search = tree.query_incremental_ANN(points[curr_max_w]);
+      auto search_it = search.begin();
+      while (search_it != search.end() && search_it->second <= curr_max_dist ) {
+        //std::cout << search_it->second << " " << curr_max_dist << "\n";
+        if (dist_to_L[search_it->first] > search_it->second)
+          dist_to_L[search_it->first] = search_it->second;
+        search_it++;
+      }
+      // choose the next curr_max_w
+      curr_max_dist = 0;
+      for (unsigned i = 0; i < dist_to_L.size(); i++)
+        if (dist_to_L[i] > curr_max_dist) {
+          curr_max_dist = dist_to_L[i];
+          curr_max_w = i;
+        }
+    }
+  }
+  
+  template < typename Kernel,
+             typename Point_container,
+             typename OutputIterator>
+  void choose_by_farthest_point( Kernel& k,
+                                 Point_container const &points,
+                                 int final_size,
+                                 OutputIterator output_it)
+  {    
+    // Choose randomly the first landmark 
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_int_distribution<> dis(1, 6);
+    int starting_point = dis(gen);
+
+    choose_by_farthest_point_old(k, points, final_size, starting_point, output_it);
+  }
+
+  
 } // namespace subsampling
   
 }  // namespace Gudhi
diff --git a/src/Subsampling/test/test_choose_farthest_point.cpp b/src/Subsampling/test/test_choose_farthest_point.cpp
index 87c2c38d..dff2cd4e 100644
--- a/src/Subsampling/test/test_choose_farthest_point.cpp
+++ b/src/Subsampling/test/test_choose_farthest_point.cpp
@@ -24,14 +24,15 @@
 // # define TBB_USE_THREADING_TOOL
 // #endif
 
-#define BOOST_TEST_DYN_LINK
-#define BOOST_TEST_MODULE "witness_complex_points"
-#include <boost/test/unit_test.hpp>
-#include <boost/mpl/list.hpp>
+// #define BOOST_TEST_DYN_LINK
+// #define BOOST_TEST_MODULE "test_choose_farthest_point"
+//#include <boost/test/unit_test.hpp>
+//#include <boost/mpl/list.hpp>
 
 #include <gudhi/choose_by_farthest_point.h>
 #include <vector>
 #include <iterator>
+#include <gudhi/Clock.h>
 
 #include <CGAL/Epick_d.h>
 
@@ -40,18 +41,52 @@ typedef typename K::FT                                            FT;
 typedef typename K::Point_d                                       Point_d;
 
 
-BOOST_AUTO_TEST_CASE(test_choose_farthest_point) {
-  std::vector< Point_d > points, results;
-  // Add grid points (625 points)
-  for (FT i = 0; i < 5; i += 1.0)
-    for (FT j = 0; j < 5; j += 1.0)
-      for (FT k = 0; k < 5; k += 1.0)
-        for (FT l = 0; l < 5; l += 1.0)
+//BOOST_AUTO_TEST_CASE(test_choose_farthest_point)
+int main() {
+  std::vector< Point_d > points, results, results2;
+  K k;
+  Clock t;
+  // Add grid points (810000 points)
+  for (FT i = 0; i < 30; i += 1.0)
+    for (FT j = 0; j < 30; j += 1.0)
+      for (FT k = 0; k < 30; k += 1.0)
+        for (FT l = 0; l < 30; l += 1.0)
           points.push_back(Point_d(std::vector<FT>({i, j, k, l})));
 
-  results.clear();
-  K k;
-  Gudhi::subsampling::choose_by_farthest_point(k, points, 100, std::back_inserter(results));
+  unsigned final_size = 100, numeral = 1;
+  std::cout << "Test   New     Old\n";
+  while (final_size < 100001) {
+    std::cout << final_size << ": ";
+    results.clear();
+    t.begin();
+    Gudhi::subsampling::choose_by_farthest_point(k, points, final_size, 0, std::back_inserter(results));
+    t.end();
+    std::cout << t.num_seconds()  << " s, ";
   
-  assert(results.size() == 100);
+  // std::cout << "New algorithm result:\n";
+  // for (auto p: results)
+  //   std::cout << p << std::endl;
+  
+    results2.clear();
+    t.begin();
+    Gudhi::subsampling::choose_by_farthest_point_old(k, points, final_size, 0, std::back_inserter(results2));
+    t.end();
+    std::cout << t.num_seconds()  << " s" << std::endl;
+
+  
+  // std::cout << "Old algorithm result:\n";
+  // for (auto p: results2)
+  //   std::cout << p << std::endl;
+  
+    assert(results.size() == final_size);
+    assert(results2.size() == final_size);
+    assert(results == results2);
+
+    switch (numeral) {
+    case 1: numeral = 2; final_size *= 2; break;
+    case 2: numeral = 5; final_size = final_size/2*5; break;
+    case 5: numeral = 1; final_size *= 2; break;
+    default: assert(false);
+    }
+  }
 }