Calculating the volume of a tetrahedron

Question

Here's my function for calculating the volume of a tetrahedron. I've tried to comment well and perform a few checks on the types of the objects passed to it. Can the comments be improved?

I'd also like an opinion about my method for checking the object types. Is it better to use assert for this?

# eg9-tetrahedron.py
import numpy as np

def tetrahedron_volume(vertices=None, sides=None):
    """
    Return the volume of the tetrahedron with given vertices or sides. If
    vertices are given they must be in a NumPy array with shape (4,3): the
    position vectors of the 4 vertices in 3 dimensions; if the six sides are
    given, they must be an array of length 6. If both are given, the sides
    will be used in the calculation.

    Raises a ValueError if the vertices do not form a tetrahedron (for example,
    because they are coplanar, colinear or coincident). This method implements
    Tartaglia's formula using the Cayley-Menger determinant:
              |0   1    1    1    1  |
              |1   0   s1^2 s2^2 s3^2|
    288 V^2 = |1  s1^2  0   s4^2 s5^2|
              |1  s2^2 s4^2  0   s6^2|
              |1  s3^2 s5^2 s6^2  0  |
    where s1, s2, ..., s6 are the tetrahedron side lengths.

    Warning: this algorithm has not been tested for numerical stability.

    """

    # The indexes of rows in the vertices array corresponding to all
    # possible pairs of vertices
    vertex_pair_indexes = np.array(((0, 1), (0, 2), (0, 3),
                                    (1, 2), (1, 3), (2, 3)))
    if sides is None:
        # If no sides were provided, work them out from the vertices
        if type(vertices) != np.ndarray or vertices.shape != (4,3):
            raise TypeError('Invalid vertex array in tetrahedron_volume():'
                             ' vertices must be a numpy array with shape (4,3)')
        # Get all the squares of all side lengths from the differences between
        # the 6 different pairs of vertex positions
        vertex1, vertex2 = vertex_pair_indexes[:,0], vertex_pair_indexes[:,1]
        sides_squared = np.sum((vertices[vertex1] - vertices[vertex2])**2,
                               axis=-1)
    else:
        # Check that sides has been provided as a valid array and square it
        if type(sides) != np.ndarray or sides.shape != (6,):
            raise TypeError('Invalid argument to tetrahedron_volume():'
                             ' sides must be a numpy array with shape (6,)')
        sides_squared = sides**2

    # Set up the Cayley-Menger determinant
    M = np.zeros((5,5))
    # Fill in the upper triangle of the matrix
    M[0,1:] = 1
    # The squared-side length elements can be indexed using the vertex
    # pair indices (compare with the determinant illustrated above)
    M[tuple(zip(*(vertex_pair_indexes + 1)))] = sides_squared

    # The matrix is symmetric, so we can fill in the lower triangle by
    # adding the transpose
    M = M + M.T

    # Calculate the determinant and check it is positive (negative or zero
    # values indicate the vertices to not form a tetrahedron).
    det = np.linalg.det(M)
    if det <= 0:
        raise ValueError('Provided vertices do not form a tetrahedron')
    return np.sqrt(det / 288)

Answer 1

Your docstring includes implementation details; typically this should be in comments instead since it's not part of the API. Doing so gives implementation flexibility, too.

Instead of checking if type(vertices) != np.ndarray, you should do

vertices = np.asarray(vertices)

which wraps or converts non-np.ndarray types and does nothing to already-correct types. This is a form of duck-typing. ALternatively, you could not check nor call asarray on the type, which would be friendlier in some use-cases (fake ndarray-like types) and less friendly to others. The convention in Numpy seems to be to use asarray, although both are reasonable choices.

You should check that only one of vertices or sides is passed in.

I would be hesitant adding details like

Invalid ... in tetrahedron_volume()

to error messages: the traceback tells you where it's raised.

Instead of

vertex1, vertex2 = vertex_pair_indexes[:,0], vertex_pair_indexes[:,1]

You can just do

vertex1, vertex2 = vertex_pair_indexes.T

Even better would be

vertex1, vertex2 = np.triu_indices(4, k=1)

To deal with the second usage, you can do

M[1:, 1:][np.triu_indices(4, k=1)] = sides_squared

Finally, it might be a little nicer to start with a triangular, rather than an empty, matrix:

M = np.tri(5, k=-1).T

I wuld also make it take positional-only arguments if Python 2 support isn't required:

def tetrahedron_volume(*, vertices=None, sides=None):

---

Larger simplifications

After considering WolframMathWorld's explanation, you can just do

distances = scipy.spatial.distance.pdist(points, metric='sqeuclidean')

to get pairwise distances, and make the matrix with

distances_square = scipy.spatial.distance.squareform(distances)

plus a bit of concatenation.

This gives the opportunity to implement the N-dimensional variant by using the algorithm as given since the prior is so simple and supports it trivially. This requires addition of:

num_verts = distance.num_obs_y(sq_dists)
coeff = - (-2) ** num_verts * factorial(num_verts) ** 2

and of doing the division before checking if vol_square <= 0 instead of hardcoding a coefficient of +288.

If this means we get the same amount of code doing more general stuff with better error checking, great!

from math import factorial

import numpy as np
from scipy.spatial import distance

def simplex_volume(*, vertices=None, sides=None) -> float:
    """
    Return the volume of the simplex with given vertices or sides.

    If vertices are given they must be in a NumPy array with shape (N+1, N):
    the position vectors of the N+1 vertices in N dimensions. If the sides
    are given, they must be the compressed pairwise distance matrix as
    returned from scipy.spatial.distance.pdist.

    Raises a ValueError if the vertices do not form a simplex (for example,
    because they are coplanar, colinear or coincident).

    Warning: this algorithm has not been tested for numerical stability.
    """

    # Implements http://mathworld.wolfram.com/Cayley-MengerDeterminant.html

    if (vertices is None) == (sides is None):
        raise ValueError("Exactly one of vertices and sides must be given")

    # β_ij = |v_i - v_k|²
    if sides is None:
        vertices = np.asarray(vertices, dtype=float)
        sq_dists = distance.pdist(vertices, metric='sqeuclidean')

    else:
        sides = np.asarray(sides, dtype=float)
        if not distance.is_valid_y(sides):
            raise ValueError("Invalid number or type of side lengths")

        sq_dists = sides ** 2

    # Add border while compressed
    num_verts = distance.num_obs_y(sq_dists)
    bordered = np.concatenate((np.ones(num_verts), sq_dists)) 

    # Make matrix and find volume
    sq_dists_mat = distance.squareform(bordered)

    coeff = - (-2) ** (num_verts-1) * factorial(num_verts-1) ** 2
    vol_square = np.linalg.det(sq_dists_mat) / coeff

    if vol_square <= 0:
        raise ValueError('Provided vertices do not form a tetrahedron')

    return np.sqrt(vol_square)

Answer 2

"The most Pythonic way to check the type of an object is... not to check it." The pythonic way to do this is to assume that the calling code/user is sensible enough to pass an object which can be handled by the callee. This allows for subclass items or simply any object with the correct properties to be used in the code.

For maintainability you should:

• use self explanatory names such as numpy, matrix, determinant etc.,
• use descriptively named constants for any magic numbers such as 288,
• split the function into one for sides and one for vertices (one may end up calling the other or they may both call a third, internal function for the final calculation),
• run the code through pep8 and possibly a cyclomatic complexity checker like radon, and
• in general, go through all comments, and for every comment see if you can refactor the code to make the comment unnecessary. For example, instead of "Set up the Cayley-Menger determinant" create a function get_cayley_menger_determinant.