Calculating the volume of a tetrahedron 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 &lt;= 0:
        raise ValueError('Provided vertices do not form a tetrahedron')
    return np.sqrt(det / 288)
"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.