Skillversitas

Testing and Debugging

Testing NumPy code differs from testing ordinary Python in two ways: floating-point values need tolerance-based comparison rather than ==, and the dominant bug class is shape errors, not logic errors. The np.testing module is built for exactly this: assertions that compare arrays elementwise with tolerances and report a useful diff on failure. Combined with pytest for structure and a fixed seed for reproducibility, it is the complete testing workflow for numerical code.

The np.testing Assertions

import numpy as np

a = np.array([1.0, 2.0, 3.0])
b = np.array([1.0, 2.0, 3.0 + 1e-15])

# exact equality - for integers and provably identical values
np.testing.assert_array_equal(a, a)

# tolerance-based - the correct float comparison
np.testing.assert_allclose(a, b)     # passes: within 1e-15

# explicit tolerance control
np.testing.assert_allclose(a, a + 1e-12, rtol=1e-10, atol=0)

# a deliberate failure - the message shows where
try:
    np.testing.assert_allclose(a, a + 1.0)
except AssertionError as e:
    print(str(e)[:60])

The three assertions cover the practical cases: assert_array_equal for exact comparison, assert_allclose for floating-point with relative and absolute tolerances, and assert_almost_equal for a decimal-places specification. On failure they print a per-element diff - that diff is the debugging signal, showing which element differed and by how much rather than merely that something was wrong.

Structuring Tests with pytest

import numpy as np
import pytest

def normalize_rows(x):
    s = x.sum(axis=1, keepdims=True)
    return x / s

def test_normalize_rows_sum_to_one():
    x = np.array([[1.0, 2.0], [3.0, 0.0]])
    out = normalize_rows(x)
    np.testing.assert_allclose(out.sum(axis=1), [1.0, 1.0])

def test_normalize_shape_preserved():
    x = np.random.default_rng(0).normal(size=(7, 3))
    assert normalize_rows(x).shape == x.shape

@pytest.mark.parametrize("n", [1, 4, 100])
def test_diagonal_identity(n):
    out = np.eye(n) @ np.eye(n)
    np.testing.assert_array_equal(out, np.eye(n))

The structure is pytest's: one small test per behavior, named for the property it guarantees. Every numerical test asserts two things - the shape is preserved (assert out.shape == x.shape) and the values are correct within tolerance (assert_allclose). parametrize sweeps the size - a small case, a typical case, and an edge case - because shape bugs are scale-dependent: code that passes at n=4 often breaks at n=1 or n=100.

Reproducibility: The Fixed Seed

import numpy as np

# a test on random data needs a fixed seed to be reproducible
def test_gram_matrix_is_spd():
    rng = np.random.default_rng(123)           # fixed: failure reproducible
    X = rng.normal(size=(50, 5))
    G = X.T @ X + 1e-6 * np.eye(5)

    assert np.allclose(G, G.T)                  # symmetric
    w = np.linalg.eigvalsh(G)
    assert (w > 0).all()                         # positive definite

# seed inside the test, not globally: tests stay independent,
# and a failure points at one deterministic input

A test that uses random data without a seed is a flaky test - it fails on one run and passes on the next, and the failure cannot be reproduced. The fix is a default_rng(seed) local to the test: the input is random in form but fixed in value, so any failure is deterministic and debuggable. The alternative, np.random.seed at import time, couples the tests to global state and breaks parallel execution.

Debugging: Shapes First

import numpy as np

# the first reflex on any array bug: print the shapes
def pipeline(x):
    y = x @ x.T
    print(x.shape, y.shape)        # the drift is visible the moment it happens
    return y

# control the noise: limit how much NumPy prints
np.printoptions(precision=4, threshold=10, linewidth=80)
big = np.random.default_rng(0).normal(size=(5, 5))
print(big)                         # 4 decimals, compact output

# quantify the divergence when two results "should" match
a = np.linspace(0, 1, 50)
b = np.linspace(0, 1, 50, endpoint=False)
print(np.abs(a - b).max())         # one number, not an eyeball estimate

The debugging order for numerical code: (1) print the shapes of every intermediate - most failures are a shape that drifted somewhere in the pipeline; (2) check the dtype when values look "off" (an integer division truncating, a float32 losing precision); (3) quantify the actual divergence with np.abs(a - b).max() instead of trusting your eye. printoptions keeps large-array output readable while you hunt. And because the comparison is tolerance-based, "close enough" is a number you state - the max divergence is the metric, not a feeling.

The Pitfalls That Break Tests

  • Comparing floats with == - a == b on a float array is an exact elementwise comparison and almost never what you want; use allclose or assert_allclose.
  • Forgetting the shape assertion - a test that checks values but not the shape passes on the wrong axis and fails at production scale.
  • Seeding globally - np.random.seed at import couples tests to order and breaks parallelism; use a local default_rng.
  • Testing one size only - shape bugs are scale-dependent; parametrize the size across small, typical, and edge cases.

Summary

Test NumPy code with np.testing: assert_array_equal for exact values, assert_allclose for floats, and read the per-element diff on failure as the debugging signal. Structure it in pytest - one behavior per test, always asserting both shape and value, parametrized over sizes - and make every random test reproducible with a local fixed default_rng seed. When debugging, print shapes first, check the dtype second, and quantify divergence with abs(a - b).max(). The two habits that prevent most numerical bugs: never compare floats with ==, and never let a shape be implicit.