4
\$\begingroup\$

I have implemented my code using Cython. It is the current bottleneck in my computations.

There are two non-numpy functions involved:

  1. calculate_2D_dist_squared which calculates the distance squared between two points
  2. 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 >= 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 < 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] < 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.

\$\endgroup\$
1
\$\begingroup\$

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).

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.