triangle_vol.py

import numpy as np

import perfplot


def dot(a, b):
    a_shape = a.shape
    b_shape = b.shape
    b = b.reshape(b.shape[0], -1)
    return np.dot(a, b).reshape(a_shape[:-1] + b_shape[1:])


def setup(n):
    pts = np.random.rand(3, n, 2)
    e = np.array(
        [
            pts[2] - pts[1],
            pts[0] - pts[2],
            pts[1] - pts[0],
        ]
    )
    ei_dot_ei = np.einsum("...k, ...k->...", e, e)
    ei_dot_ej = ei_dot_ei - np.sum(ei_dot_ei, axis=0) / 2
    s1 = ei_dot_ei[0] + ei_dot_ei[1] + ei_dot_ei[2]
    s2 = np.sum(ei_dot_ei ** 2, axis=0)
    return pts, ei_dot_ei, ei_dot_ej, s1, s2


def test1(data):
    _, _, ei_dot_ej, _, _ = data
    # ei_dot_ej = ei_dot_ei - np.sum(ei_dot_ei, axis=0) / 2
    out = 0.25 * (
        ei_dot_ej[2] * ei_dot_ej[0]
        + ei_dot_ej[0] * ei_dot_ej[1]
        + ei_dot_ej[1] * ei_dot_ej[2]
    )
    return out


def test2(data):
    _, ei_dot_ei, _, _, _ = data
    return (0.5 * np.sum(ei_dot_ei, axis=0) ** 2 - np.sum(ei_dot_ei ** 2, axis=0)) / 8


def test3(data):
    _, ei_dot_ei, ei_dot_ej, _, _ = data
    return -np.einsum("ij,ij->j", ei_dot_ei, ei_dot_ej) / 8


def test4(data):
    pts, ei_dot_ei, ei_dot_ej, _, _ = data

    A = pts[1:] - pts[0]
    ATA = np.einsum("i...j,k...j->...ik", A, A)
    e = np.ones(ATA.shape[:-1])
    sol = np.sum(np.linalg.solve(ATA, e), axis=1)
    out = 0.25 / sol * ei_dot_ei[0]
    return out


# same as test4, but with explicit solve()
def test5(data):
    pts, ei_dot_ei, ei_dot_ej, _, _ = data

    A = pts[1:] - pts[0]
    ATA = np.einsum("i...j,k...j->ik...", A, A)
    det = ATA[0][0] * ATA[1][1] - ATA[1][0] * ATA[0][1]
    sol = (ATA[0][0] + ATA[1][1] - ATA[1][0] - ATA[0][1]) / det
    out = 0.25 / sol * ei_dot_ei[0]
    return out


# cayley-menger
def test6(data):
    pts, ei_dot_ei, ei_dot_ej, _, _ = data
    z = np.zeros(ei_dot_ei.shape[1])
    e = np.ones(ei_dot_ei.shape[1])
    cm = np.array(
        [
            [z, e, e, e],
            [e, z, ei_dot_ei[0], ei_dot_ei[1]],
            [e, ei_dot_ei[0], z, ei_dot_ei[2]],
            [e, ei_dot_ei[1], ei_dot_ei[2], z],
        ]
    )
    cm = np.moveaxis(cm, -1, 0)
    out = -np.linalg.det(cm) / 16
    return out


perfplot.save(
    "vol.png",
    setup=setup,
    kernels=[test1, test2, test3, test4, test5, test6],
    n_range=[2 ** k for k in range(5, 23)],
)