Add parallelism using OpenMP (#49)

.travis.yml

@@ -56,6 +56,12 @@ jobs:
       os: osx
       osx_image: xcode11.2
       language: shell
+      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
+      env:
+        - CC="/usr/local/opt/llvm/bin/clang -L/usr/local/opt/llvm/lib -Wl,-rpath,/usr/local/opt/llvm/lib"
+        - CXX="/usr/local/opt/llvm/bin/clang++ -L/usr/local/opt/llvm/lib -Wl,-rpath,/usr/local/opt/llvm/lib"
     - stage: lint
       language: python
       install: pip install flake8

README.rst

@@ -182,3 +182,13 @@ Notes
 
 * In order for the algorithm to work the mesh must not be numbered incorrectly
   or disconnected or of somehow degenerate.
+
+* The c++ compiler must have OpenMP installed. This is generally not an issue
+  on g++ or Microsoft Visual Studio. However, in macOS one may need to install
+  llvm from brew in order to use OpenMP: ``brew install llvm``. Then, change the
+  ``CC`` and ``CXX`` environment variables to the following:
+
+  .. code-block:: txt
+
+     CC="/usr/local/opt/llvm/bin/clang -L/usr/local/opt/llvm/lib -Wl,-rpath,/usr/local/opt/llvm/lib"
+     CXX="/usr/local/opt/llvm/bin/clang++ -L/usr/local/opt/llvm/lib -Wl,-rpath,/usr/local/opt/llvm/lib"

gdist.pyx

@@ -61,61 +61,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,
@@ -245,43 +203,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)
+    cdef Py_ssize_t N = vertices.shape[0]
 
-    # Define and create object for exact algorithm on that mesh:
-    cdef GeodesicAlgorithmExact *algorithm = new GeodesicAlgorithmExact(&amesh)
+    cdef vector[double] points
+    cdef vector[unsigned] faces
 
-    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)
-
-    return scipy.sparse.csc_matrix((data, (rows, columns)), shape=(N, N))
+    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(
+        (distances.data, (distances.rows, distances.columns)), shape=(N, N)
+    )
 
 
 def distance_matrix_of_selected_points(

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,6 +41,7 @@
 """
 
 import os
+import sys
 import setuptools
 
 import numpy
@@ -66,9 +67,16 @@
         name=GEODESIC_NAME,  # Name of extension
         sources=["gdist.pyx"],  # Filename of Cython source
         language="c++",  # Cython create C++ source
+        # Disable assertions; one is failing geodesic_mesh.h:405
         define_macros=define_macros,
-        extra_compile_args=["--std=c++14"],
-        extra_link_args=["--std=c++14"],
+        extra_compile_args=[
+            '--std=c++14',
+            '/openmp' if sys.platform == 'win32' else '-fopenmp',
+        ],
+        extra_link_args=[
+            '--std=c++14',
+            '/openmp' if sys.platform == 'win32' else '-fopenmp',
+        ],
         include_dirs=[numpy.get_include(), "geodesic_library"],
     )
 ]