Skip to content

Identities, round trips and sibling code

Question: does the code agree with itself and with textbook identities, where no exact solution is needed?

Most of the bugs found in widely used libraries during the corpus study were of three kinds, and none of them needed an order study:

  • Two code paths that must agree, but do not. In one acoustics library, Eyring's reverberation formula had the air-absorption term with the wrong sign, while Sabine's formula right next to it was correct. Two related functions differed in the same way elsewhere: a max_moment that ignored a flag its min_moment sibling passed on, and a zone plate that ignored its design wavelength while the binary version used it.
  • An operation that its inverse does not undo. A solar-image rotation by zero time moved the map by half a pixel, because the reference pixel was written as 1 + n/2 instead of (n + 1)/2.
  • A textbook identity that fails. Poisson's equation for a potential and its density, unitarity of a lossless scattering matrix, E² = p² + m², a sliding average that returns to the mean on stationary noise.

physverify.identities has a small check for each. All of them return a Result; wrap one in Check("limit", label, lambda: ...) to run it with the other checks.

"""Checks with no exact solution: sibling agreement, a round trip and an analytic derivative."""

import math

import numpy as np

from physverify.identities import check_agree, check_derivative, check_roundtrip

C = 343.0  # speed of sound, m/s


def rt60_sabine(S, V, a, m):
    return 24 * math.log(10) / C * V / (a * S + 4 * m * V)


def rt60_eyring(S, V, a, m):  # bug: the air term has the wrong sign
    return -(24 * math.log(10) / C) * V / (S * math.log(1 - a) + 4 * m * V)


# Eyring's formula reduces to Sabine's for small wall absorption, with or without air absorption.
cases = [(94.0, 60.0, 1e-4, m) for m in (0.0, 0.005, 0.05)]
print(check_agree(rt60_eyring, rt60_sabine, cases, rtol=1e-3, names=("eyring", "sabine")).status, "sibling agreement")


# A reference pixel written as 1 + n/2 instead of (n + 1)/2: pixel -> sky -> pixel is off by half a pixel.
n = 20
to_sky = lambda p: (p + 1 - (1 + n / 2)) * 2.0  # noqa: E731
to_pixel = lambda s: s / 2.0 + (n + 1) / 2 - 1  # noqa: E731
r = check_roundtrip(to_sky, to_pixel, [(n - 1) / 2])
print(r.status, "round trip:", r.message.split("; ")[1].split(" (")[0])


# An analytic derivative against finite differences of the function itself.
f = lambda x: np.exp(-x) * np.sin(3 * x)  # noqa: E731
df = lambda x: np.exp(-x) * (3 * np.cos(3 * x) - np.sin(3 * x))  # noqa: E731
print(check_derivative(f, df, [0.0, 0.4, 2.0]).status, "derivative")

Output:

FAIL sibling agreement
FAIL round trip: got 9.0, expected 9.5
PASS derivative
Function What it checks
check_agree(a, b, inputs, rtol, atol) a(*x) == b(*x) for each input: sibling paths, two formulas that coincide in a limit
check_roundtrip(forward, inverse, inputs) inverse(forward(x)) == x
check_derivative(f, df, points, order=1) an analytic first or second derivative against finite differences of f
check_identity([(label, lhs, rhs), ...]) identities you compute from the physics
check_unitary(M) M^H M = I, for one matrix or a stack
fd_derivative, laplacian_fd the finite-difference helpers, for building your own identities

Values may be real or complex, scalars or arrays; a non-finite value is a FAIL.

What these checks cannot tell you

A round trip or a sibling comparison uses two pieces of the code under test, so it passes when both share the same mistake. Use them to find inconsistencies, and keep at least one check against a reference written from the physics. A derivative check settles whether a derivative matches its function, not whether the function is right.

A derivative check also works the other way: it can show that a suspected bug is not one. A code reader claimed that an activity-coefficient derivative in a thermodynamics library dropped a term; comparing it with finite differences of the activity coefficients agreed to 1e-10, so no issue was filed.