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.pyimport numpy as npdef 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)
Calculating the volume of a tetrahedron
python;validation;mathematics;numpy
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 dovertices = 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 likeInvalid ... in tetrahedron_volume()to error messages: the traceback tells you where it's raised.Instead ofvertex1, vertex2 = vertex_pair_indexes[:,0], vertex_pair_indexes[:,1]You can just dovertex1, vertex2 = vertex_pair_indexes.TEven better would bevertex1, vertex2 = np.triu_indices(4, k=1)To deal with the second usage, you can doM[1:, 1:][np.triu_indices(4, k=1)] = sides_squaredFinally, it might be a little nicer to start with a triangular, rather than an empty, matrix:M = np.tri(5, k=-1).TI wuld also make it take positional-only arguments if Python 2 support isn't required:def tetrahedron_volume(*, vertices=None, sides=None):Larger simplificationsAfter considering WolframMathWorld's explanation, you can just dodistances = scipy.spatial.distance.pdist(points, metric='sqeuclidean')to get pairwise distances, and make the matrix withdistances_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) ** 2and 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 factorialimport numpy as npfrom scipy.spatial import distancedef 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)