Add faster c implementations for some functions in triangulation.py

[philippeitis]

I implemented some of the functions in triangulation.py in C for a considerable speed boost in several situations. They may not support all list types that are used, but I tried to implement functionality for the ones I saw.

The timing differences are as follows, for 100000 iterations of arbitrarily chosen function arguments:

C fast_norm: 0.025422s
Python fast_norm: 0.074730s
C fast_2d_circumcircle: 0.035352s
Python fast_2d_circumcircle: 1.077767s
C fast_3d_circumcircle: 0.068189s
Python fast_3d_circumcircle: 1.729581s
C fast_2d_point_in_simplex: 0.037834s
Python fast_2d_point_in_simplex: 0.170977s

The difference is considerably large for fast_2d_circumcircle and fast_3d_circumcircle, as these use numpy array operations that massively slow the program down - getting rid of

    points = np.array(points)
    pts = points[1:] - points[0]

in both of these functions would provide a considerable speedup - converting points to an array creates a lot of unnecessary overhead.

[basnijholt]

This is awesome! I'll check it out in more detail when I'm near my computer.

Do you know if this speeds up the LearnerND significantly in running ~10.000 points?

[basnijholt]

@philippeitis, I hope you don't mind, but I solved the style issues in setup.py.

Also, when trying it out, I seem to get the following error SystemError: <built-in function fast_3d_circumcircle> returned NULL without setting an error

When running

import adaptive
import numpy as np

def sphere(xyz):
    x, y, z = xyz
    a = 0.4
    return x + z**2 + np.exp(-(x**2 + y**2 + z**2 - 0.75**2)**2/a**4)

learner = adaptive.LearnerND(sphere, bounds=[(-1, 1), (-1, 1), (-1, 1)])
runner = adaptive.runner.simple(learner, goal=lambda l: l.npoints > 2000)

To use your module in the LearnerND I've applied the following patch:

diff --git a/adaptive/learner/triangulation.py b/adaptive/learner/triangulation.py
index 0cc0bde..4e29af1 100644
--- a/adaptive/learner/triangulation.py
+++ b/adaptive/learner/triangulation.py
@@ -7,28 +7,29 @@ from math import factorial
 import numpy as np
 import scipy.spatial
 
+from adaptive.triangulation import fast_norm, fast_2d_circumcircle, fast_3d_circumcircle, fast_2d_point_in_simplex
 
-def fast_norm(v):
-    # notice this method can be even more optimised
-    if len(v) == 2:
-        return math.sqrt(v[0] * v[0] + v[1] * v[1])
-    if len(v) == 3:
-        return math.sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2])
-    return math.sqrt(np.dot(v, v))
+# def fast_norm(v):
+#     # notice this method can be even more optimised
+#     if len(v) == 2:
+#         return math.sqrt(v[0] * v[0] + v[1] * v[1])
+#     if len(v) == 3:
+#         return math.sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2])
+#     return math.sqrt(np.dot(v, v))
 
 
-def fast_2d_point_in_simplex(point, simplex, eps=1e-8):
-    (p0x, p0y), (p1x, p1y), (p2x, p2y) = simplex
-    px, py = point
+# def fast_2d_point_in_simplex(point, simplex, eps=1e-8):
+#     (p0x, p0y), (p1x, p1y), (p2x, p2y) = simplex
+#     px, py = point
 
-    area = 0.5 * (-p1y * p2x + p0y * (p2x - p1x) + p1x * p2y + p0x * (p1y - p2y))
+#     area = 0.5 * (-p1y * p2x + p0y * (p2x - p1x) + p1x * p2y + p0x * (p1y - p2y))
 
-    s = 1 / (2 * area) * (+p0y * p2x + (p2y - p0y) * px - p0x * p2y + (p0x - p2x) * py)
-    if s < -eps or s > 1 + eps:
-        return False
-    t = 1 / (2 * area) * (+p0x * p1y + (p0y - p1y) * px - p0y * p1x + (p1x - p0x) * py)
+#     s = 1 / (2 * area) * (+p0y * p2x + (p2y - p0y) * px - p0x * p2y + (p0x - p2x) * py)
+#     if s < -eps or s > 1 + eps:
+#         return False
+#     t = 1 / (2 * area) * (+p0x * p1y + (p0y - p1y) * px - p0y * p1x + (p1x - p0x) * py)
 
-    return (t >= -eps) and (s + t <= 1 + eps)
+#     return (t >= -eps) and (s + t <= 1 + eps)
 
 
 def point_in_simplex(point, simplex, eps=1e-8):
@@ -42,80 +43,80 @@ def point_in_simplex(point, simplex, eps=1e-8):
     return all(alpha > -eps) and sum(alpha) < 1 + eps
 
 
-def fast_2d_circumcircle(points):
-    """Compute the center and radius of the circumscribed circle of a triangle
-
-    Parameters
-    ----------
-    points: 2D array-like
-        the points of the triangle to investigate
-
-    Returns
-    -------
-    tuple
-        (center point : tuple(int), radius: int)
-    """
-    points = np.array(points)
-    # transform to relative coordinates
-    pts = points[1:] - points[0]
-
-    (x1, y1), (x2, y2) = pts
-    # compute the length squared
-    l1 = x1 * x1 + y1 * y1
-    l2 = x2 * x2 + y2 * y2
-
-    # compute some determinants
-    dx = +l1 * y2 - l2 * y1
-    dy = -l1 * x2 + l2 * x1
-    aa = +x1 * y2 - x2 * y1
-    a = 2 * aa
-
-    # compute center
-    x = dx / a
-    y = dy / a
-    radius = math.sqrt(x * x + y * y)  # radius = norm([x, y])
-
-    return (x + points[0][0], y + points[0][1]), radius
-
-
-def fast_3d_circumcircle(points):
-    """Compute the center and radius of the circumscribed shpere of a simplex.
-
-    Parameters
-    ----------
-    points: 2D array-like
-        the points of the triangle to investigate
-
-    Returns
-    -------
-    tuple
-        (center point : tuple(int), radius: int)
-    """
-    points = np.array(points)
-    pts = points[1:] - points[0]
-
-    (x1, y1, z1), (x2, y2, z2), (x3, y3, z3) = pts
-
-    l1 = x1 * x1 + y1 * y1 + z1 * z1
-    l2 = x2 * x2 + y2 * y2 + z2 * z2
-    l3 = x3 * x3 + y3 * y3 + z3 * z3
-
-    # Compute some determinants:
-    dx = +l1 * (y2 * z3 - z2 * y3) - l2 * (y1 * z3 - z1 * y3) + l3 * (y1 * z2 - z1 * y2)
-    dy = +l1 * (x2 * z3 - z2 * x3) - l2 * (x1 * z3 - z1 * x3) + l3 * (x1 * z2 - z1 * x2)
-    dz = +l1 * (x2 * y3 - y2 * x3) - l2 * (x1 * y3 - y1 * x3) + l3 * (x1 * y2 - y1 * x2)
-    aa = +x1 * (y2 * z3 - z2 * y3) - x2 * (y1 * z3 - z1 * y3) + x3 * (y1 * z2 - z1 * y2)
-    a = 2 * aa
-
-    center = (dx / a, -dy / a, dz / a)
-    radius = fast_norm(center)
-    center = (
-        center[0] + points[0][0],
-        center[1] + points[0][1],
-        center[2] + points[0][2],
-    )
-
-    return center, radius
+# def fast_2d_circumcircle(points):
+#     """Compute the center and radius of the circumscribed circle of a triangle
+
+#     Parameters
+#     ----------
+#     points: 2D array-like
+#         the points of the triangle to investigate
+
+#     Returns
+#     -------
+#     tuple
+#         (center point : tuple(int), radius: int)
+#     """
+#     points = np.array(points)
+#     # transform to relative coordinates
+#     pts = points[1:] - points[0]
+
+#     (x1, y1), (x2, y2) = pts
+#     # compute the length squared
+#     l1 = x1 * x1 + y1 * y1
+#     l2 = x2 * x2 + y2 * y2
+
+#     # compute some determinants
+#     dx = +l1 * y2 - l2 * y1
+#     dy = -l1 * x2 + l2 * x1
+#     aa = +x1 * y2 - x2 * y1
+#     a = 2 * aa
+
+#     # compute center
+#     x = dx / a
+#     y = dy / a
+#     radius = math.sqrt(x * x + y * y)  # radius = norm([x, y])
+
+#     return (x + points[0][0], y + points[0][1]), radius
+
+
+# def fast_3d_circumcircle(points):
+#     """Compute the center and radius of the circumscribed shpere of a simplex.
+
+#     Parameters
+#     ----------
+#     points: 2D array-like
+#         the points of the triangle to investigate
+
+#     Returns
+#     -------
+#     tuple
+#         (center point : tuple(int), radius: int)
+#     """
+#     points = np.array(points)
+#     pts = points[1:] - points[0]
+
+#     (x1, y1, z1), (x2, y2, z2), (x3, y3, z3) = pts
+
+#     l1 = x1 * x1 + y1 * y1 + z1 * z1
+#     l2 = x2 * x2 + y2 * y2 + z2 * z2
+#     l3 = x3 * x3 + y3 * y3 + z3 * z3
+
+#     # Compute some determinants:
+#     dx = +l1 * (y2 * z3 - z2 * y3) - l2 * (y1 * z3 - z1 * y3) + l3 * (y1 * z2 - z1 * y2)
+#     dy = +l1 * (x2 * z3 - z2 * x3) - l2 * (x1 * z3 - z1 * x3) + l3 * (x1 * z2 - z1 * x2)
+#     dz = +l1 * (x2 * y3 - y2 * x3) - l2 * (x1 * y3 - y1 * x3) + l3 * (x1 * y2 - y1 * x2)
+#     aa = +x1 * (y2 * z3 - z2 * y3) - x2 * (y1 * z3 - z1 * y3) + x3 * (y1 * z2 - z1 * y2)
+#     a = 2 * aa
+
+#     center = (dx / a, -dy / a, dz / a)
+#     radius = fast_norm(center)
+#     center = (
+#         center[0] + points[0][0],
+#         center[1] + points[0][1],
+#         center[2] + points[0][2],
+#     )
+
+#     return center, radius
 
 
 def fast_det(matrix):

[philippeitis]

I've updated the functions and their docstrings, so calling help(adaptive.triangulation.{}) will yield useful information, and I added helpful error messages where necessary.

fast_norm can handle one dimensional lists, tuples and numpy arrays. I've decided not to support arrays with more dimensions, as numpy's linear algebra library (eg. np.linalg.norm) is likely to be faster and more robust for these cases.

fast_2d_circumcircle can handle lists of tuples, tuples of tuples, and numpy arrays that have at least 3 rows and exactly 2 columns - numpy.array([(1, 0), (0, 1), (-1, 0)]) is automatically converted to such an array, but this may not be the case if you're appending tuples.
fast_3d_circumcircle can handle lists of tuples, tuples of tuples, and numpy arrays that have at least 4 rows and exactly 3 columns (so numpy.array([(0, 0, 0), (0, 0, 0), (0, 0, 0), (0, 0, 0)])). For these functions, I can also add 2x3 and 3x4 support, though it should be noted that these cases are highly unlikely (and should never occur if the function is called from circumsphere).

fast_2d_point_in_simplex takes a tuple as its first argument, a list of at least three tuples as its second argument, and a double as its third. I can add support for numpy arrays and lists here if necessary.