Add parallelism using OpenMP

Closes #10

.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

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(

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());
         }

geodesic_library/geodesic_utils.h

@@ -1,7 +1,15 @@
 #include <vector>
+#include <omp.h>
 
 #include "geodesic_algorithm_exact.h"
 
+class SparseMatrix {
+public:
+    std::vector<unsigned> rows, columns;
+    std::vector<double> data;
+};
+
+
 std::vector<double> compute_gdist_impl(
     std::vector<double> points,
     std::vector<unsigned> faces,
@@ -35,3 +43,48 @@ std::vector<double> compute_gdist_impl(
     }
     return distances;
 }
+
+SparseMatrix local_gdist_matrix_impl(
+    std::vector<double> points,
+    std::vector<unsigned> 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<unsigned> rows, columns;
+    std::vector<double> data;
+    #pragma omp parallel
+    {
+        geodesic::GeodesicAlgorithmExact algorithm(&mesh);
+        std::vector<unsigned> rows_private, columns_private;
+        std::vector<double> data_private;
+        #pragma omp for nowait
+        for (int i = 0; i < (int)mesh.vertices().size(); ++i) {
+            std::vector<geodesic::SurfacePoint> 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;
+}

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"],
     )
 ]