From bf3a53d0d73f57fb801875a5dde11b9ede10b0aa Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Fri, 17 Mar 2017 13:11:22 +0000 Subject: 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 --- src/cython/CMakeLists.txt | 2 +- .../cython/euclidean_strong_witness_complex.pyx | 97 ++++++++++++++++++++++ src/cython/cython/euclidean_witness_complex.pyx | 13 +++ ...ex_diagram_persistence_from_off_file_example.py | 78 +++++++++++++++++ src/cython/include/Alpha_complex_interface.h | 2 - .../Euclidean_strong_witness_complex_interface.h | 87 +++++++++++++++++++ .../include/Euclidean_witness_complex_interface.h | 46 +++------- src/cython/test/test_euclidean_witness_complex.py | 34 +++++++- 8 files changed, 319 insertions(+), 40 deletions(-) create mode 100644 src/cython/cython/euclidean_strong_witness_complex.pyx create mode 100755 src/cython/example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py create mode 100644 src/cython/include/Euclidean_strong_witness_complex_interface.h 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 . +""" + +__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 . +""" + +__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>&points) { alpha_complex_ = new Alpha_complex(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 . + */ + +#ifndef EUCLIDEAN_STRONG_WITNESS_COMPLEX_INTERFACE_H +#define EUCLIDEAN_STRONG_WITNESS_COMPLEX_INTERFACE_H + +#include +#include + +#include "Simplex_tree_interface.h" + +#include + +#include +#include // std::pair +#include +#include + +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>&landmarks, std::vector>&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 get_point(unsigned vh) { + std::vector 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 landmarks_; + std::vector witnesses_; + Euclidean_strong_witness_complex 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&landmarks, std::vector&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 landmarks_; - std::vector witnesses_; - Euclidean_witness_complex witness_complex_; -}; - -/*template> -class Euclidean_witness_complex_interface : public Euclidean_witness_complex { -//class Euclidean_witness_complex_interface { - using Point_d = Kernel::Point_d; - - public: - Euclidean_witness_complex_interface(const std::vector>& landmarks, const std::vector>& witnesses) { - : landmarks_(std::begin(landmarks), std::end(landmarks)), - witnesses_(std::begin(witnesses), std::end(witnesses)), - Gudhi::witness_complex::Euclidean_witness_complex(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 get_point(unsigned vh) { + std::vector 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* witness_complex_; std::vector landmarks_; std::vector witnesses_; - -};*/ + Euclidean_witness_complex 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)] + -- cgit v1.2.3