Euclidean strong witness complex cytonization

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

Former-commit-id: de15e40b8f0efcc59b2f07ffb5f5d59ccab8d51c

8 files changed, 319 insertions, 40 deletions

diff --git a/src/cython/CMakeLists.txt b/src/cython/CMakeLists.txt
index efe2005c..25cf952e 100644
--- a/src/cython/CMakeLists.txt
+++ b/src/cython/CMakeLists.txt
@@ -128,7 +128,7 @@ if(PYTHON_PATH AND CYTHON_PATH)
     if (NOT CGAL_VERSION VERSION_LESS 4.6.0)
       # If CGAL_VERSION >= 4.6.0, include euclidean versions of witness complex
       set(GUDHI_CYTHON_EUCLIDEAN_WITNESS_COMPLEX
-          "include 'cython/euclidean_witness_complex.pyx'\n")
+          "include 'cython/euclidean_witness_complex.pyx'\ninclude 'cython/euclidean_strong_witness_complex.pyx'\n")
     else (NOT CGAL_VERSION VERSION_LESS 4.6.0)
       # Remove alpha complex unitary tests
       file(REMOVE ${CMAKE_CURRENT_BINARY_DIR}/test/test_euclidean_witness_complex.py)
diff --git a/src/cython/cython/euclidean_strong_witness_complex.pyx b/src/cython/cython/euclidean_strong_witness_complex.pyx
new file mode 100644
index 00000000..73d31708
--- /dev/null
+++ b/src/cython/cython/euclidean_strong_witness_complex.pyx
@@ -0,0 +1,97 @@
+from cython cimport numeric
+from libcpp.vector cimport vector
+from libcpp.utility cimport pair
+
+"""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):       Vincent Rouvreau
+
+   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/>.
+"""
+
+__author__ = "Vincent Rouvreau"
+__copyright__ = "Copyright (C) 2016 INRIA"
+__license__ = "GPL v3"
+
+cdef extern from "Euclidean_strong_witness_complex_interface.h" namespace "Gudhi":
+    cdef cppclass Euclidean_strong_witness_complex_interface "Gudhi::witness_complex::Euclidean_strong_witness_complex_interface":
+        Euclidean_strong_witness_complex_interface(vector[vector[double]] landmarks, vector[vector[double]] witnesses)
+        void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square)
+        void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square,
+            unsigned limit_dimension)
+        vector[double] get_point(unsigned vertex)
+
+# EuclideanStrongWitnessComplex python interface
+cdef class EuclideanStrongWitnessComplex:
+    """EuclideanStrongWitnessComplex is a simplicial complex constructed from ...
+
+    """
+
+    cdef Euclidean_strong_witness_complex_interface * thisptr
+
+    # Fake constructor that does nothing but documenting the constructor
+    def __init__(self, landmarks=None, witnesses=None):
+        """WitnessComplex constructor.
+
+        :param landmarks: A list of landmarks (in the point cloud).
+        :type landmarks: list of list of double
+
+        :param witnesses: The point cloud.
+        :type witnesses: list of list of double
+        """
+
+    # The real cython constructor
+    def __cinit__(self, landmarks=None, witnesses=None):
+        if landmarks is not None and witnesses is not None:
+            self.thisptr = new Euclidean_strong_witness_complex_interface(landmarks, witnesses)
+
+    def __dealloc__(self):
+        if self.thisptr != NULL:
+            del self.thisptr
+
+    def __is_defined(self):
+        """Returns true if WitnessComplex pointer is not NULL.
+         """
+        return self.thisptr != NULL
+
+    def create_simplex_tree(self, max_alpha_square, limit_dimension = -1):
+        """
+        :param max_alpha_square: The maximum alpha square threshold the
+            simplices shall not exceed. Default is set to infinity.
+        :type max_alpha_square: float
+        :returns: A simplex tree created from the Delaunay Triangulation.
+        :rtype: SimplexTree
+        """
+        simplex_tree = SimplexTree()
+        if limit_dimension is not -1:
+            self.thisptr.create_simplex_tree(simplex_tree.thisptr, max_alpha_square, limit_dimension)
+        else:
+            self.thisptr.create_simplex_tree(simplex_tree.thisptr, max_alpha_square)
+        return simplex_tree
+
+    def get_point(self, vertex):
+        """This function returns the point corresponding to a given vertex.
+
+        :param vertex: The vertex.
+        :type vertex: int.
+        :returns:  The point.
+        :rtype: list of float
+        """
+        cdef vector[double] point = self.thisptr.get_point(vertex)
+        return point
+
diff --git a/src/cython/cython/euclidean_witness_complex.pyx b/src/cython/cython/euclidean_witness_complex.pyx
index efbbf06d..6197c322 100644
--- a/src/cython/cython/euclidean_witness_complex.pyx
+++ b/src/cython/cython/euclidean_witness_complex.pyx
@@ -34,6 +34,7 @@ cdef extern from "Euclidean_witness_complex_interface.h" namespace "Gudhi":
         void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square)
         void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square,
             unsigned limit_dimension)
+        vector[double] get_point(unsigned vertex)
 
 # EuclideanWitnessComplex python interface
 cdef class EuclideanWitnessComplex:
@@ -82,3 +83,15 @@ cdef class EuclideanWitnessComplex:
         else:
             self.thisptr.create_simplex_tree(simplex_tree.thisptr, max_alpha_square)
         return simplex_tree
+
+    def get_point(self, vertex):
+        """This function returns the point corresponding to a given vertex.
+
+        :param vertex: The vertex.
+        :type vertex: int.
+        :returns:  The point.
+        :rtype: list of float
+        """
+        cdef vector[double] point = self.thisptr.get_point(vertex)
+        return point
+
diff --git a/src/cython/example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py b/src/cython/example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py
new file mode 100755
index 00000000..6c3e3a8a
--- /dev/null
+++ b/src/cython/example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py
@@ -0,0 +1,78 @@
+#!/usr/bin/env python
+
+import gudhi
+import argparse
+
+"""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):       Vincent Rouvreau
+
+   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/>.
+"""
+
+__author__ = "Vincent Rouvreau"
+__copyright__ = "Copyright (C) 2016 INRIA"
+__license__ = "GPL v3"
+
+parser = argparse.ArgumentParser(description='EuclideanStrongWitnessComplex creation from '
+                                 'points read in a OFF file.',
+                                 epilog='Example: '
+                                 'example/witness_complex_diagram_persistence_from_off_file_example.py '
+                                 '-f ../data/points/tore3D_300.off -a 1.0 -n 20 -d 2'
+                                 '- Constructs a alpha complex with the '
+                                 'points from the given OFF file.')
+parser.add_argument("-f", "--file", type=str, required=True)
+parser.add_argument("-a", "--max_alpha_square", type=float, required=True)
+parser.add_argument("-n", "--number_of_landmarks", type=int, required=True)
+parser.add_argument("-d", "--limit_dimension", type=int, required=True)
+parser.add_argument('--no-diagram', default=False, action='store_true' , help='Flag for not to display the diagrams')
+
+args = parser.parse_args()
+
+with open(args.file, 'r') as f:
+    first_line = f.readline()
+    if (first_line == 'OFF\n') or (first_line == 'nOFF\n'):
+        print("#####################################################################")
+        print("EuclideanStrongWitnessComplex creation from points read in a OFF file")
+
+        witnesses = gudhi.read_off(off_file=args.file)
+        landmarks = gudhi.pick_n_random_points(points=witnesses, nb_points=args.number_of_landmarks)
+
+        message = "EuclideanStrongWitnessComplex with max_edge_length=" + repr(args.max_alpha_square) + \
+            " - Number of landmarks=" + repr(args.number_of_landmarks)
+        print(message)
+
+        witness_complex = gudhi.EuclideanStrongWitnessComplex(witnesses=witnesses, landmarks=landmarks)
+        simplex_tree = witness_complex.create_simplex_tree(max_alpha_square=args.max_alpha_square,
+            limit_dimension=args.limit_dimension)
+
+        message = "Number of simplices=" + repr(simplex_tree.num_simplices())
+        print(message)
+
+        diag = simplex_tree.persistence()
+
+        print("betti_numbers()=")
+        print(simplex_tree.betti_numbers())
+
+        if args.no_diagram == False:
+            gudhi.diagram_persistence(diag)
+
+    else:
+        print(args.file, "is not a valid OFF file")
+
+    f.close()
diff --git a/src/cython/include/Alpha_complex_interface.h b/src/cython/include/Alpha_complex_interface.h
index 9f308ae9..bcd2849b 100644
--- a/src/cython/include/Alpha_complex_interface.h
+++ b/src/cython/include/Alpha_complex_interface.h
@@ -41,8 +41,6 @@ class Alpha_complex_interface {
   using Dynamic_kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >;
   using Point_d = Dynamic_kernel::Point_d;
 
-  typedef typename Simplex_tree<>::Simplex_key Simplex_key;
-
  public:
   Alpha_complex_interface(std::vector<std::vector<double>>&points) {
     alpha_complex_ = new Alpha_complex<Dynamic_kernel>(points);
diff --git a/src/cython/include/Euclidean_strong_witness_complex_interface.h b/src/cython/include/Euclidean_strong_witness_complex_interface.h
new file mode 100644
index 00000000..739014b9
--- /dev/null
+++ b/src/cython/include/Euclidean_strong_witness_complex_interface.h
@@ -0,0 +1,87 @@
+/*    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):       Vincent Rouvreau
+ *
+ *    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/>.
+ */
+
+#ifndef EUCLIDEAN_STRONG_WITNESS_COMPLEX_INTERFACE_H
+#define	EUCLIDEAN_STRONG_WITNESS_COMPLEX_INTERFACE_H
+
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Euclidean_strong_witness_complex.h>
+
+#include "Simplex_tree_interface.h"
+
+#include <CGAL/Epick_d.h>
+
+#include <vector>
+#include <utility>  // std::pair
+#include <iostream>
+#include <cstddef>
+
+namespace Gudhi {
+
+namespace witness_complex {
+
+
+class Euclidean_strong_witness_complex_interface {
+  using Dynamic_kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >;
+  using Point_d = Dynamic_kernel::Point_d;
+
+  typedef typename Simplex_tree<>::Simplex_key Simplex_key;
+
+ public:
+  Euclidean_strong_witness_complex_interface(std::vector<std::vector<double>>&landmarks, std::vector<std::vector<double>>&witnesses)
+    : landmarks_(landmarks.begin(), landmarks.end()),
+      witnesses_(witnesses.begin(), witnesses.end()),
+      witness_complex_(landmarks_, witnesses_) {
+  }
+
+  void create_simplex_tree(Gudhi::Simplex_tree<>* simplex_tree, double max_alpha_square, std::size_t limit_dimension) {
+    witness_complex_.create_complex(*simplex_tree, max_alpha_square, limit_dimension);
+    simplex_tree->initialize_filtration();
+  }
+
+  void create_simplex_tree(Gudhi::Simplex_tree<>* simplex_tree, double max_alpha_square) {
+    witness_complex_.create_complex(*simplex_tree, max_alpha_square);
+    simplex_tree->initialize_filtration();
+  }
+
+  std::vector<double> get_point(unsigned vh) {
+    std::vector<double> vd;
+    if (vh < landmarks_.size()) {
+      Point_d ph = witness_complex_.get_point(vh);
+      for (auto coord = ph.cartesian_begin(); coord < ph.cartesian_end(); coord++)
+        vd.push_back(*coord);
+    }
+    return vd;
+  }
+
+ private:
+  std::vector<Point_d> landmarks_;
+  std::vector<Point_d> witnesses_;
+  Euclidean_strong_witness_complex<Dynamic_kernel> witness_complex_;
+};
+
+}  // namespace witness_complex
+
+} // namespace Gudhi
+
+#endif  // EUCLIDEAN_STRONG_WITNESS_COMPLEX_INTERFACE_H
+
diff --git a/src/cython/include/Euclidean_witness_complex_interface.h b/src/cython/include/Euclidean_witness_complex_interface.h
index f79d93ab..efa8885f 100644
--- a/src/cython/include/Euclidean_witness_complex_interface.h
+++ b/src/cython/include/Euclidean_witness_complex_interface.h
@@ -52,11 +52,6 @@ class Euclidean_witness_complex_interface {
       witnesses_(witnesses.begin(), witnesses.end()),
       witness_complex_(landmarks_, witnesses_) {
   }
-  Euclidean_witness_complex_interface(std::vector<Point_d>&landmarks, std::vector<Point_d>&witnesses)
-    : landmarks_(landmarks.begin(), landmarks.end()),
-      witnesses_(witnesses.begin(), witnesses.end()),
-      witness_complex_(landmarks_, witnesses_) {
-  }
 
   void create_simplex_tree(Gudhi::Simplex_tree<>* simplex_tree, double max_alpha_square, std::size_t limit_dimension) {
     witness_complex_.create_complex(*simplex_tree, max_alpha_square, limit_dimension);
@@ -68,42 +63,21 @@ class Euclidean_witness_complex_interface {
     simplex_tree->initialize_filtration();
   }
 
- private:
-  std::vector<Point_d> landmarks_;
-  std::vector<Point_d> witnesses_;
-  Euclidean_witness_complex<Dynamic_kernel> witness_complex_;
-};
-
-/*template<typename Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>>
-class Euclidean_witness_complex_interface : public Euclidean_witness_complex<Kernel> {
-//class Euclidean_witness_complex_interface {
-  using Point_d = Kernel::Point_d;
-
- public:
-  Euclidean_witness_complex_interface(const std::vector<std::vector<double>>& landmarks, const std::vector<std::vector<double>>& witnesses) {
-    : landmarks_(std::begin(landmarks), std::end(landmarks)),
-      witnesses_(std::begin(witnesses), std::end(witnesses)),
-      Gudhi::witness_complex::Euclidean_witness_complex<Kernel>(landmarks_, witnesses_) {
-  }
-
-  void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double  max_alpha_square,
-                           std::size_t limit_dimension) {
-    create_complex(*simplex_tree, max_alpha_square, limit_dimension);
-    simplex_tree->initialize_filtration();
-  }
-
-  void create_simplex_tree(Simplex_tree_interface<>* simplex_tree,
-                           double  max_alpha_square) {
-    create_complex(*simplex_tree, max_alpha_square);
-    simplex_tree->initialize_filtration();
+  std::vector<double> get_point(unsigned vh) {
+    std::vector<double> vd;
+    if (vh < landmarks_.size()) {
+      Point_d ph = witness_complex_.get_point(vh);
+      for (auto coord = ph.cartesian_begin(); coord < ph.cartesian_end(); coord++)
+        vd.push_back(*coord);
+    }
+    return vd;
   }
 
  private:
-//  Euclidean_witness_complex<Kernel>* witness_complex_;
   std::vector<Point_d> landmarks_;
   std::vector<Point_d> witnesses_;
-
-};*/
+  Euclidean_witness_complex<Dynamic_kernel> witness_complex_;
+};
 
 }  // namespace witness_complex
 
diff --git a/src/cython/test/test_euclidean_witness_complex.py b/src/cython/test/test_euclidean_witness_complex.py
index cb150b7e..0947cc09 100755
--- a/src/cython/test/test_euclidean_witness_complex.py
+++ b/src/cython/test/test_euclidean_witness_complex.py
@@ -34,6 +34,38 @@ def test_empty_euclidean_witness_complex():
 def test_witness_complex():
     point_cloud = [[1.0, 1.0], [7.0, 0.0], [4.0, 6.0], [9.0, 6.0],
         [0.0, 14.0], [2.0, 19.0], [9.0, 17.0]]
-    landmarks = gudhi.pick_n_random_points(points=point_cloud,nb_points=3)
+    landmarks = [[1.0, 1.0], [7.0, 0.0], [4.0, 6.0]]
     euclidean_witness_complex = gudhi.EuclideanWitnessComplex(landmarks=landmarks, witnesses = point_cloud)
     simplex_tree = euclidean_witness_complex.create_simplex_tree(max_alpha_square=4.1)
+
+    assert landmarks[0] == euclidean_witness_complex.get_point(0)
+    assert landmarks[1] == euclidean_witness_complex.get_point(1)
+    assert landmarks[2] == euclidean_witness_complex.get_point(2)
+
+    assert simplex_tree.get_filtered_tree() == [([0], 0.0), ([1], 0.0),
+        ([0, 1], 0.0), ([2], 0.0), ([0, 2], 0.0), ([1, 2], 0.0),
+        ([0, 1, 2], 0.0)]
+
+def test_empty_euclidean_strong_witness_complex():
+    euclidean_strong_witness = gudhi.EuclideanStrongWitnessComplex()
+    assert euclidean_strong_witness.__is_defined() == False
+
+def test_strong_witness_complex():
+    point_cloud = [[1.0, 1.0], [7.0, 0.0], [4.0, 6.0], [9.0, 6.0],
+        [0.0, 14.0], [2.0, 19.0], [9.0, 17.0]]
+    landmarks = [[1.0, 1.0], [7.0, 0.0], [4.0, 6.0]]
+    euclidean_strong_witness_complex = gudhi.EuclideanStrongWitnessComplex(landmarks=landmarks, witnesses = point_cloud)
+    simplex_tree = euclidean_strong_witness_complex.create_simplex_tree(max_alpha_square=14.9)
+
+    assert landmarks[0] == euclidean_strong_witness_complex.get_point(0)
+    assert landmarks[1] == euclidean_strong_witness_complex.get_point(1)
+    assert landmarks[2] == euclidean_strong_witness_complex.get_point(2)
+
+    assert simplex_tree.get_filtered_tree() == [([0], 0.0), ([1], 0.0), ([2], 0.0)]
+
+    simplex_tree = euclidean_strong_witness_complex.create_simplex_tree(max_alpha_square=100.0)
+
+    assert simplex_tree.get_filtered_tree() == [([0], 0.0), ([1], 0.0),
+    ([2], 0.0), ([1, 2], 15.0), ([0, 2], 34.0), ([0, 1], 37.0),
+    ([0, 1, 2], 37.0)]
+