From a06a649e591bc0e82b87a916e6eea16653608f43 Mon Sep 17 00:00:00 2001 From: ayan-b Date: Wed, 29 Jul 2020 22:01:23 +0530 Subject: [PATCH] Add parallelism using OpenMP Closes #10 --- .travis.yml | 14 ++- gdist.pyx | 105 +++++--------------- geodesic_library/geodesic_algorithm_exact.h | 2 +- geodesic_library/geodesic_utils.h | 53 ++++++++++ setup.py | 21 +++- 5 files changed, 107 insertions(+), 88 deletions(-) diff --git a/.travis.yml b/.travis.yml index 4fe9937..95414b8 100644 --- a/.travis.yml +++ b/.travis.yml @@ -3,7 +3,6 @@ notifications: install: - pip3 install . - # - python setup.py build_ext --inplace script: - pip3 install pytest~=3.6.1 - pytest @@ -25,6 +24,19 @@ jobs: os: osx osx_image: xcode11.2 language: shell + before_cache: + - brew cleanup; + # Credit https://discourse.brew.sh/t/best-practice-for-homebrew-on-travis-brew-update-is-5min-to-build-time/5215/9 + # Cache only .git files under "/usr/local/Homebrew" so "brew update" does not take 5min every build + - find /usr/local/Homebrew \! -regex ".+\.git.+" -delete; + cache: + directories: + - $HOME/Library/Caches/Homebrew + - /usr/local/Homebrew + before_install: + - brew pin gdal; brew pin postgis; brew pin cairo; brew pin glib; brew pin gnutls; brew pin p11-kit; brew pin gnupg; brew pin poppler; + - brew install llvm + - echo 'export PATH="/usr/local/opt/llvm/bin:$PATH"' >> ~/.bash_profile - stage: lint language: python python: 3.8 diff --git a/gdist.pyx b/gdist.pyx index dadc28e..271accf 100644 --- a/gdist.pyx +++ b/gdist.pyx @@ -62,61 +62,19 @@ from libcpp.vector cimport vector ################################################################################ ############ Defininitions to access the C++ geodesic distance library ######### ################################################################################ -cdef extern from "geodesic_mesh_elements.h" namespace "geodesic": - cdef cppclass Vertex: - Vertex() - -cdef extern from "geodesic_mesh_elements.h" namespace "geodesic": - cdef cppclass SurfacePoint: - SurfacePoint() - SurfacePoint(Vertex*) - double& x() - double& y() - double& z() - -cdef extern from "geodesic_mesh.h" namespace "geodesic": - cdef cppclass Mesh: - Mesh() - void initialize_mesh_data(vector[double]&, vector[unsigned]&, bool) - vector[Vertex]& vertices() - cdef extern from "geodesic_utils.h": + cdef cppclass SparseMatrix: + vector[unsigned] rows + vector[unsigned] columns + vector[double] data vector[double] compute_gdist_impl(vector[double], vector[unsigned], vector[unsigned], vector[unsigned], double, bool, bool) - -cdef extern from "geodesic_algorithm_exact.h" namespace "geodesic": - cdef cppclass GeodesicAlgorithmExact: - GeodesicAlgorithmExact(Mesh*) - void propagate(vector[SurfacePoint]&, double, vector[SurfacePoint]*) - unsigned best_source(SurfacePoint&, double&) + SparseMatrix local_gdist_matrix_impl(vector[double], vector[unsigned], double, bool) cdef extern from "geodesic_constants_and_simple_functions.h" namespace "geodesic": double GEODESIC_INF ################################################################################ -cdef get_mesh( - numpy.ndarray[numpy.float64_t, ndim=2] vertices, - numpy.ndarray[numpy.int32_t, ndim=2] triangles, - Mesh &amesh, - bool is_one_indexed -): - # Define C++ vectors to contain the mesh surface components. - cdef vector[double] points - cdef vector[unsigned] faces - - # Map numpy array of mesh "vertices" to C++ vector of mesh "points" - cdef numpy.float64_t coord - for coord in vertices.flatten(): - points.push_back(coord) - - # Map numpy array of mesh "triangles" to C++ vector of mesh "faces" - cdef numpy.int32_t indx - for indx in triangles.flatten(): - faces.push_back(indx) - - amesh.initialize_mesh_data(points, faces, is_one_indexed) - - def compute_gdist(numpy.ndarray[numpy.float64_t, ndim=2] vertices, numpy.ndarray[numpy.int32_t, ndim=2] triangles, numpy.ndarray[numpy.int32_t, ndim=1] source_indices = None, @@ -246,43 +204,26 @@ def local_gdist_matrix(numpy.ndarray[numpy.float64_t, ndim=2] vertices, from the propgate step... """ - cdef Mesh amesh - get_mesh(vertices, triangles, amesh, is_one_indexed) - - # Define and create object for exact algorithm on that mesh: - cdef GeodesicAlgorithmExact *algorithm = new GeodesicAlgorithmExact(&amesh) - - cdef vector[SurfacePoint] source, targets cdef Py_ssize_t N = vertices.shape[0] - cdef Py_ssize_t k - cdef Py_ssize_t kk - cdef numpy.float64_t distance = 0 - - # Add all vertices as targets - for k in range(N): - targets.push_back(SurfacePoint(&amesh.vertices()[k])) - - rows = [] - columns = [] - data = [] - for k in range(N): - source.push_back(SurfacePoint(&amesh.vertices()[k])) - algorithm.propagate(source, max_distance, NULL) - source.pop_back() - - for kk in range(N): #TODO: Reduce to vertices reached during propagate. - algorithm.best_source(targets[kk], distance) - - if ( - distance is not GEODESIC_INF - and distance is not 0 - and distance <= max_distance - ): - rows.append(k) - columns.append(kk) - data.append(distance) + + cdef vector[double] points + cdef vector[unsigned] faces + + for k in vertices.flatten(): + points.push_back(k) + for k in triangles.flatten(): + faces.push_back(k) + + cdef SparseMatrix distances = local_gdist_matrix_impl( + points, + faces, + max_distance, + is_one_indexed, + ) - return scipy.sparse.csc_matrix((data, (rows, columns)), shape=(N, N)) + return scipy.sparse.csc_matrix( + (distances.data, (distances.rows, distances.columns)), shape=(N, N) + ) def distance_matrix_of_selected_points( diff --git a/geodesic_library/geodesic_algorithm_exact.h b/geodesic_library/geodesic_algorithm_exact.h index fad8c41..2f2027b 100644 --- a/geodesic_library/geodesic_algorithm_exact.h +++ b/geodesic_library/geodesic_algorithm_exact.h @@ -1163,7 +1163,7 @@ inline void GeodesicAlgorithmExact::construct_propagated_intervals(bool invert, p->min() = 0; - assert(p->start() < p->stop()); + // assert(p->start() < p->stop()); assert(p->start() >= 0.0); assert(p->stop() <= edge->length()); } diff --git a/geodesic_library/geodesic_utils.h b/geodesic_library/geodesic_utils.h index c2dedc6..67c21ed 100644 --- a/geodesic_library/geodesic_utils.h +++ b/geodesic_library/geodesic_utils.h @@ -1,7 +1,15 @@ #include +#include #include "geodesic_algorithm_exact.h" +class SparseMatrix { +public: + std::vector rows, columns; + std::vector data; +}; + + std::vector compute_gdist_impl( std::vector points, std::vector faces, @@ -35,3 +43,48 @@ std::vector compute_gdist_impl( } return distances; } + +SparseMatrix local_gdist_matrix_impl( + std::vector points, + std::vector faces, + double max_distance, + bool is_one_indexed +) { + geodesic::Mesh mesh; + mesh.initialize_mesh_data(points, faces, is_one_indexed=is_one_indexed); // create internal mesh data structure including edges + + SparseMatrix distances; + + std::vector rows, columns; + std::vector data; + #pragma omp parallel + { + geodesic::GeodesicAlgorithmExact algorithm(&mesh); + std::vector rows_private, columns_private; + std::vector data_private; + #pragma omp for nowait + for (int i = 0; i < (int)mesh.vertices().size(); ++i) { + std::vector sources {&mesh.vertices()[i]}; + algorithm.propagate(sources, geodesic::GEODESIC_INF, NULL); // cover the whole mesh + for(int j = 0; j < (int)mesh.vertices().size(); ++j) { + geodesic::SurfacePoint p(&mesh.vertices()[j]); + + double distance; + unsigned best_source = algorithm.best_source(p, distance); // for a given surface point, find closest source and distance to this source + + if (distance != 0 && distance <= geodesic::GEODESIC_INF && distance <= max_distance) { + rows_private.push_back(i); + columns_private.push_back(j); + data_private.push_back(distance); + } + } + } + rows.insert(rows.end(), rows_private.begin(), rows_private.end()); + columns.insert(columns.end(), columns_private.begin(), columns_private.end()); + data.insert(data.end(), data_private.begin(), data_private.end()); + } + distances.rows = rows; + distances.columns = columns; + distances.data = data; + return distances; +} diff --git a/setup.py b/setup.py index 102a2f1..ce0f44b 100644 --- a/setup.py +++ b/setup.py @@ -41,11 +41,18 @@ """ import os -import numpy +import sys import shutil import setuptools + +import numpy from Cython.Distutils import build_ext + +if sys.platform == "darwin": + os.environ["CC"] = "/usr/local/opt/llvm/bin/clang -L/usr/local/opt/llvm/lib -Wl,-rpath,/usr/local/opt/llvm/lib" + os.environ["CXX"] = "/usr/local/opt/llvm/bin/clang++ -L/usr/local/opt/llvm/lib -Wl,-rpath,/usr/local/opt/llvm/lib" + GEODESIC_NAME = "gdist" GEODESIC_MODULE = [ @@ -54,9 +61,15 @@ sources=["gdist.pyx"], # Filename of Cython source language="c++", # Cython create C++ source # Disable assertions; one is failing geodesic_mesh.h:405 - define_macros=[('NDEBUG', 1)], - extra_compile_args=['--std=c++14'], - extra_link_args=['--std=c++14'], + # define_macros=[('NDEBUG', 1)], + extra_compile_args=[ + '--std=c++14', + '-fopenmp' if sys.platform != 'win32' else '/openmp', + ], + extra_link_args=[ + '--std=c++14', + '-fopenmp' if sys.platform != 'win32' else '/openmp', + ], include_dirs=[numpy.get_include(), "geodesic_library"], ) ]