Calculating the distance squared between all vertex-pairs of a number of 2D polygons I have implemented my code using Cython. It is the current bottleneck in my computations.

There are two non-numpy functions involved:


calculate_2D_dist_squared which calculates the distance squared between two points
calculate_2D_dist_squared_matrix which generates the distances squared between every two distinct vertices. calculate_2D_dist_squared_matrix organizes results so that dist_squared_matrix[1, 2, 3, 4] = "the distance between polygon 1, vertex 2 and polygon 3, vertex 4". All indexing starts from 0.




cdef double calculate_2D_dist_squared(self, np.ndarray[np.float64_t, ndim=1] p1, np.ndarray[np.float64_t, ndim=1] p2):
    cdef np.ndarray[np.float64_t, ndim=1] relative_vector = p1 - p2

    return relative_vector[0]**2 + relative_vector[1]**2




cdef np.ndarray[np.float64_t, ndim=4] calculate_2D_dist_squared_matrix(self, np.ndarray[np.float64_t, ndim=3] polygons_vertex_coords, int num_polygons, int num_vertices):
    cdef:
        int pi_focus
        int vi_focus
        int pi
        int vi
        # at initialization, set all dist_squared values to be -1,
        # indicating that they have been initialized, but not set properly
        # since by definition a dist_squared value has to be &gt;= 0
        np.ndarray[np.float64_t, ndim=4] result  = -1*np.ones((num_polygons, num_vertices, num_polygons, num_vertices), dtype=np.float64)

    for pi_focus in range(num_polygons):
        for vi_focus in range(num_vertices):
            for pi in range(num_polygons):
                for vi in range(num_vertices):
                    # if a dist_squared &lt; 0, then it means that it
                    # it has not been changed since initialization, and 
                    # needs to be updated, this way I avoid repeating work
                    if result[pi_focus, vi_focus, pi, vi] &lt; 0:
                        dist_squared = self.calculate_2D_dist_squared(polygons_vertex_coords[pi_focus, vi_focus], polygons_vertex_coords[pi, vi])
                        result[pi_focus, vi_focus, pi, vi] = dist_squared
                        result[pi, vi, pi_focus, vi_focus] = dist_squared

    return result


What are some things I could think about in order to increase the performance of my code?

For the time being, I got a significant improvement in performance by only re-calculating updates to the dist_squared_matrix, rather than always re-calculating the dist_squared_matrix entirely every step.
It would be helpful to have a simple test harness for the code.

In reviewing your code I can see some possible bottlenecks, in Cython there can be quite a large overhead in passing arrays (views) to functions, also the array operations may be convenient but when performing them on very small arrays of only a few values performance will tend to suffer considerably compared with doing the arithmetic directly — Cython can make some optimizations when it comes to calling functions into python-space and a few numpy array operations bypass python-space — but arithmetic on basic types in Cython is guaranteed to be at the speed of c as it translates directly to c code, as does array indexing.

Hence you would probably gain some performance benefits by passing floats directly instead of arrays.

cdef double calculate_2D_dist_squared(self, np.float64_t x1, np.float64_t y1, 
                                      np.float64_t x2, np.float64_t y2):
    return (x1 - x2) ** 2 + (y1 - y2) ** 2


As noted previously, indexing is very very fast in Cython, but intermediate views can result in performance losses, hence my instinct would be to call the function in this way:

dist_squared = self.calculate_2D_dist_squared(
    polygons_vertex_coords[pi_focus, vi_focus, 0], 
    polygons_vertex_coords[pi_focus, vi_focus, 1], 
    polygons_vertex_coords[pi, vi, 0],
    polygons_vertex_coords[pi, vi, 1])


Instead of the other possible way which would involve creating a variable to hold polygons_vertex_coords[pi_focus, vi_focus] (etc).