Testing NumPy Applications: From Unit Tests to Integration
NumPy is the backbone of scientific computing in Python, powering everything from data analysis pipelines to machine learning models. But numerical code has a nasty habit of failing silently — a wrong axis argument, an unintended broadcast, or a floating-point precision drift can corrupt results without raising a single exception. That's why a disciplined testing strategy is essential. This tutorial walks you through testing NumPy applications comprehensively, starting with focused unit tests and scaling up to integration tests that validate entire data pipelines.
Why Testing NumPy Code Is Different
Testing numerical code differs from testing typical business logic in several important ways. First, floating-point arithmetic is inherently imprecise, so exact equality checks are often inappropriate. Second, NumPy operations can silently broadcast arrays into unexpected shapes, producing results that "look right" but are mathematically wrong. Third, performance matters — a function that returns correct results but is 100x slower than expected may indicate a bug in how arrays are being processed.
A robust testing strategy for NumPy applications should address these concerns by combining:
- Array-aware assertions that handle floating-point tolerance and shape comparison.
- Property-based tests that generate random arrays to explore edge cases.
- Integration tests that validate multi-step pipelines end to end.
- Performance regression tests that catch accidental slowdowns.
Setting Up Your Testing Environment
Before writing tests, install pytest, numpy, and hypothesis (for property-based testing). The numpy.testing module ships with NumPy itself, so no extra dependency is needed there.
pip install numpy pytest hypothesis
A typical project layout looks like this:
my_numpy_app/
├── my_numpy_app/
│ ├── __init__.py
│ ├── stats.py
│ └── pipeline.py
├── tests/
│ ├── __init__.py
│ ├── test_stats.py
│ ├── test_stats_property.py
│ └── test_pipeline.py
├── pyproject.toml
└── pytest.ini
Configure pytest to run from the project root with a pytest.ini file:
[pytest]
testpaths = tests
addopts = -v --tb=short
Unit Testing with numpy.testing
The numpy.testing module provides a suite of assertion functions specifically designed for comparing arrays. These are far more informative than plain assert statements because they report mismatched shapes, indices of differing elements, and maximum deviations.
Key Assertion Functions
np.testing.assert_array_equal— exact equality, useful for integer arrays or boolean masks.np.testing.assert_allclose— compares with relative and absolute tolerance, the workhorse for floating-point arrays.np.testing.assert_array_almost_equal— legacy decimal-place comparison.np.testing.assert_raises— context manager for expected exceptions.np.testing.assert_array_less— element-wise less-than comparison.
A Sample Module to Test
Let's create a small statistics module that we will test throughout this tutorial:
# my_numpy_app/stats.py
import numpy as np
def zscore(data: np.ndarray, axis: int = 0) -> np.ndarray:
"""Compute z-scores along the given axis."""
mean = np.mean(data, axis=axis, keepdims=True)
std = np.std(data, axis=axis, keepdims=True)
# Guard against division by zero for constant columns
std = np.where(std == 0, 1.0, std)
return (data - mean) / std
def moving_average(data: np.ndarray, window: int) -> np.ndarray:
"""Compute a simple moving average with the given window size."""
if window <= 0:
raise ValueError("window must be a positive integer")
if window > data.shape[-1]:
raise ValueError("window cannot exceed data length")
weights = np.ones(window) / window
return np.convolve(data, weights, mode="valid")
def normalize_rows(matrix: np.ndarray) -> np.ndarray:
"""Scale each row so its values sum to 1."""
row_sums = matrix.sum(axis=1, keepdims=True)
row_sums = np.where(row_sums == 0, 1.0, row_sums)
return matrix / row_sums
Writing the Unit Tests
Now let's write unit tests for each function. Notice how we use assert_allclose for floating-point results and assert_array_equal where exact comparison is appropriate.
# tests/test_stats.py
import numpy as np
import pytest
from numpy.testing import assert_allclose, assert_array_equal, assert_raises
from my_numpy_app.stats import zscore, moving_average, normalize_rows
class TestZscore:
def test_basic_1d(self):
data = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
result = zscore(data)
assert_allclose(result.mean(), 0.0, atol=1e-10)
assert_allclose(result.std(), 1.0, atol=1e-10)
def test_axis_parameter(self):
data = np.array([[1.0, 2.0], [3.0, 4.0]])
result_axis0 = zscore(data, axis=0)
result_axis1 = zscore(data, axis=1)
# Each column should have zero mean along axis 0
assert_allclose(result_axis0.mean(axis=0), [0.0, 0.0], atol=1e-10)
# Each row should have zero mean along axis 1
assert_allclose(result_axis1.mean(axis=1), [0.0, 0.0], atol=1e-10)
def test_constant_column_does_not_explode(self):
data = np.array([[5.0, 1.0], [5.0, 2.0], [5.0, 3.0]])
result = zscore(data, axis=0)
# The constant column should produce zeros, not NaNs
assert not np.any(np.isnan(result))
assert_allclose(result[:, 0], [0.0, 0.0, 0.0], atol=1e-10)
def test_preserves_shape(self):
data = np.random.RandomState(42).randn(10, 5)
assert zscore(data).shape == data.shape
class TestMovingAverage:
def test_simple_case(self):
data = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
result = moving_average(data, window=3)
expected = np.array([2.0, 3.0, 4.0])
assert_allclose(result, expected)
def test_window_equals_length(self):
data = np.array([2.0, 4.0, 6.0])
result = moving_average(data, window=3)
assert_allclose(result, [4.0])
def test_invalid_window_raises(self):
data = np.array([1.0, 2.0, 3.0])
with assert_raises(ValueError):
moving_average(data, window=0)
with assert_raises(ValueError):
moving_average(data, window=5)
class TestNormalizeRows:
def test_rows_sum_to_one(self):
matrix = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
result = normalize_rows(matrix)
assert_allclose(result.sum(axis=1), [1.0, 1.0], atol=1e-10)
def test_zero_row_handled(self):
matrix = np.array([[0.0, 0.0], [1.0, 1.0]])
result = normalize_rows(matrix)
assert not np.any(np.isnan(result))
assert_allclose(result[0], [0.0, 0.0])
assert_allclose(result[1], [0.5, 0.5])
Run the tests with pytest tests/test_stats.py -v. Each test class groups related tests, and each test method focuses on a single behavior. This structure makes failures easy to localize.
Property-Based Testing with Hypothesis
Example-based testing is powerful, but it only checks the specific inputs you thought to write down. Property-based testing flips the model: you describe invariants that should always hold, and a library generates hundreds of random inputs to try to break them. For NumPy code, the hypothesis library with the hypothesis-extra numpy strategies is ideal.
Defining Array Strategies
Hypothesis provides arrays and array_shapes strategies that generate random NumPy arrays with controlled shapes, dtypes, and value ranges.
# tests/test_stats_property.py
import numpy as np
from hypothesis import given, strategies as st, settings, HealthCheck
from hypothesis.extra.numpy import arrays, array_shapes, floating_dtypes
from numpy.testing import assert_allclose
from my_numpy_app.stats import zscore, normalize_rows
@given(
arrays(
dtype=np.float64,
shape=array_shapes(min_dims=1, max_dims=2, min_side=2, max_side=10),
elements=st.floats(min_value=-1e6, max_value=1e6, allow_nan=False, allow_inf=False),
)
)
@settings(suppress_health_check=[HealthCheck.too_slow], max_examples=200)
def test_zscore_always_has_zero_mean_and_unit_std(arr):
result = zscore(arr, axis=-1)
# The z-score of any non-constant slice has mean ~0 and std ~1
means = result.mean(axis=-1)
stds = result.std(axis=-1)
# Allow tolerance for floating point
assert_allclose(means, np.zeros_like(means), atol=1e-8)
# Constant slices produce std 0 (guarded), so only check non-constant
non_constant = arr.std(axis=-1) > 1e-8
if np.any(non_constant):
assert_allclose(stds[non_constant], np.ones(np.sum(non_constant)), atol=1e-8)
@given(
arrays(
dtype=np.float64,
shape=st.tuples(st.integers(min_value=1, max_value=10), st.integers(min_value=2, max_value=10)),
elements=st.floats(min_value=0.1, max_value=100.0, allow_nan=False, allow_inf=False),
)
)
def test_normalize_rows_always_sums_to_one(matrix):
result = normalize_rows(matrix)
row_sums = result.sum(axis=1)
assert_allclose(row_sums, np.ones_like(row_sums), atol=1e-10)
@given(
arrays(
dtype=np.float64,
shape=array_shapes(min_dims=1, max_dims=1, min_side=3, max_side=50),
elements=st.floats(min_value=-100.0, max_value=100.0, allow_nan=False, allow_inf=False),
),
window=st.integers(min_value=1, max_value=10),
)
def test_moving_average_output_length(data, window):
from my_numpy_app.stats import moving_average
if window <= data.shape[-1]:
result = moving_average(data, window)
assert result.shape[-1] == data.shape[-1] - window + 1
Property-based tests are especially valuable for catching edge cases like empty arrays, single-element arrays, arrays with extreme values, and arrays that trigger broadcasting surprises. When Hypothesis finds a failing case, it shrinks the input to the smallest example that still fails, which is invaluable for debugging.
Testing Custom Vectorized Functions
When you write custom ufuncs or use np.vectorize, you should test both correctness and broadcasting behavior. A common bug is a function that works for 1-D input but breaks when given a 2-D array because of implicit axis assumptions.
# my_numpy_app/stats.py (continued)
def safe_log1p(x: np.ndarray) -> np.ndarray:
"""Numerically stable log(1 + x) that clips very negative inputs."""
x = np.asarray(x, dtype=np.float64)
clipped = np.clip(x, -1.0 + 1e-12, None)
return np.log1p(clipped)
# tests/test_stats.py (continued)
class TestSafeLog1p:
def test_scalar_input(self):
assert_allclose(safe_log1p(0.0), 0.0)
def test_array_input(self):
result = safe_log1p(np.array([0.0, 1.0, np.e - 1]))
assert_allclose(result, [0.0, np.log(2.0), 1.0])
def test_clips_negative_one(self):
# log(1 + (-1)) would be -inf; clipping prevents that
result = safe_log1p(np.array([-1.0, -2.0, -100.0]))
assert np.all(np.isfinite(result))
assert np.all(result <= 0.0)
def test_preserves_shape(self):
arr = np.random.RandomState(0).randn(4, 5)
assert safe_log1p(arr).shape == (4, 5)
def test_handles_integer_input(self):
result = safe_log1p(np.array([0, 1, 2]))
assert_allclose(result, [0.0, np.log(2.0), np.log(3.0)])
Integration Testing NumPy Pipelines
Unit tests verify individual functions in isolation. Integration tests verify that those functions compose correctly into a full pipeline. This is where you catch issues like dtype mismatches between stages, unexpected shape changes, and cumulative floating-point error.
A Sample Pipeline
# my_numpy_app/pipeline.py
import numpy as np
from my_numpy_app.stats import zscore, normalize_rows
def preprocess_pipeline(raw_data: np.ndarray) -> np.ndarray:
"""Full preprocessing pipeline:
1. Replace NaNs with column means.
2. Z-score normalize each column.
3. Normalize each row to sum to 1.
"""
data = raw_data.astype(np.float64, copy=True)
# Step 1: impute NaNs with column means
col_means = np.nanmean(data, axis=0)
nan_mask = np.isnan(data)
if np.any(nan_mask):
indices = np.where(nan_mask)
data[indices] = np.take(col_means, indices[1])
# Step 2: z-score columns
data = zscore(data, axis=0)
# Step 3: shift to non-negative and normalize rows
data = data - data.min(axis=1, keepdims=True)
data = normalize_rows(data)
return data
Writing the Integration Test
# tests/test_pipeline.py
import numpy as np
import pytest
from numpy.testing import assert_allclose
from my_numpy_app.pipeline import preprocess_pipeline
class TestPreprocessPipeline:
def test_end_to_end_with_clean_data(self):
raw = np.array([
[1.0, 10.0, 100.0],
[2.0, 20.0, 200.0],
[3.0, 30.0, 300.0],
])
result = preprocess_pipeline(raw)
# Output rows must sum to 1
assert_allclose(result.sum(axis=1), [1.0, 1.0, 1.0], atol=1e-10)
# Output must be non-negative
assert np.all(result >= 0.0)
# Output must have same shape
assert result.shape == raw.shape
def test_handles_nan_imputation(self):
raw = np.array([
[1.0, np.nan, 3.0],
[4.0, 5.0, 6.0],
[7.0, 8.0, 9.0],
])
result = preprocess_pipeline(raw)
assert not np.any(np.isnan(result))
assert_allclose(result.sum(axis=1), [1.0, 1.0, 1.0], atol=1e-10)
def test_preserves_dtype_as_float64(self):
raw = np.array([[1, 2, 3], [4, 5, 6]], dtype=np.int32)
result = preprocess_pipeline(raw)
assert result.dtype == np.float64
def test_single_row_input(self):
raw = np.array([[1.0, 2.0, 3.0]])
result = preprocess_pipeline(raw)
assert result.shape == (1, 3)
assert_allclose(result.sum(axis=1), [1.0], atol=1e-10)
def test_large_random_input_stability(self):
rng = np.random.RandomState(123)
raw = rng.randn(1000, 50)
result = preprocess_pipeline(raw)
assert np.all(np.isfinite(result))
assert_allclose(result.sum(axis=1), np.ones(1000), atol=1e-8)
def test_pipeline_is_deterministic(self):
raw = np.random.RandomState(99).randn(20, 5)
result1 = preprocess_pipeline(raw)
result2 = preprocess_pipeline(raw)
assert_array_equal(result1, result2)
Integration tests like these validate that the contract between functions holds. The test_large_random_input_stability test is particularly important: it catches numerical instability that only manifests at scale, such as overflow, underflow, or catastrophic cancellation.
Testing Performance and Memory
Numerical code often has performance requirements. A test that verifies correctness but ignores speed can let a regression slip through when someone accidentally introduces a Python-level loop over array elements. Use pytest-benchmark to guard against this.
pip install pytest-benchmark
# tests/test_stats_benchmark.py
import numpy as np
from my_numpy_app.stats import zscore, moving_average
def test_zscore_benchmark(benchmark):
rng = np.random.RandomState(0)
data = rng.randn(10000, 100)
result = benchmark(zscore, data, axis=0)
assert result.shape == data.shape
def test_moving_average_benchmark(benchmark):
rng = np.random.RandomState(0)
data = rng.randn(1_000_000)
result = benchmark(moving_average, data, window=100)
assert result.shape == (999_901,)
Run with pytest tests/test_stats_benchmark.py --benchmark-only. The benchmark plugin stores baseline timings and can compare runs with --benchmark-compare, making it easy to spot regressions during code review.
Best Practices for Testing NumPy Applications
1. Use the Right Tolerance
Choose rtol and atol values that match the numerical precision of your computation. For double-precision arithmetic, atol=1e-10 or atol=1e-8 is usually safe. For algorithms that accumulate error (like iterative solvers), you may need atol=1e-6 or looser. Never use exact equality (==) on floating-point arrays unless you have a specific reason.
2. Seed Your Random Generators
Always use explicit RandomState or default_rng with a fixed seed in tests. This ensures reproducibility and makes failures debuggable. Avoid np.random.seed() at the module level because it mutates global state and can cause test-order dependencies.
# Good: explicit, local, reproducible
rng = np.random.default_rng(42)
data = rng.standard_normal((100, 10))
# Bad: global state, order-dependent
np.random.seed(42)
data = np.random.randn(100, 10)
3. Test Edge Cases Explicitly
Numerical code fails on edge cases that business logic rarely encounters. Always test with:
- Empty arrays (
np.array([])) - Single-element arrays
- Arrays containing zeros, NaNs, and infinities
- Arrays with constant columns or rows
- Very large and very small values (potential overflow/underflow)
- Integer inputs where float output is expected
4. Separate Pure Logic from I/O
Design your NumPy functions to accept arrays and return arrays, with no file I/O or network calls inside. This makes them trivially testable. Put I/O operations in thin wrapper functions that you test separately, perhaps with temporary files.
import tempfile
import os
import numpy as np
def test_save_and_load_roundtrip():
arr = np.random.RandomState(0).randn(10, 5)
with tempfile.TemporaryDirectory() as tmpdir:
path = os.path.join(tmpdir, "data.npy")
np.save(path, arr)
loaded = np.load(path)
np.testing.assert_array_equal(arr, loaded)
5. Use Parametrize for Matrix Tests
When you want to test the same logic across multiple input configurations, use @pytest.mark.parametrize instead of duplicating test functions.
@pytest.mark.parametrize("shape", [(10,), (5, 5), (2, 3, 4)])
@pytest.mark.parametrize("dtype", [np.float32, np.float64])
def test_zscore_works_across_shapes_and_dtypes(shape, dtype):
rng = np.random.default_rng(0)
data = rng.standard_normal(shape).astype(dtype)
result = zscore(data, axis=-1)
assert result.shape == shape
assert result.dtype == data.dtype
6. Assert Invariants, Not Implementation Details
Test what the function guarantees, not how it achieves it. For zscore, the invariant is "output has zero mean and unit standard deviation along the specified axis." Do not assert that the implementation calls np.mean internally — that couples the test to the implementation and makes refactoring painful.
7. Document Expected Behavior with Tests
Well-named test functions serve as executable documentation. A test called test_zscore_returns_zeros_for_constant_input tells a reader exactly what happens with constant input, without requiring them to read the source. This is especially valuable for numerical functions whose behavior on edge cases may be non-obvious.
Continuous Integration Considerations
When running NumPy tests in CI, be aware of platform-specific floating-point differences. The same operation can produce slightly different results on x86 versus ARM, or across BLAS implementations. Use generous tolerances in CI and tighter tolerances locally if needed. Also pin your NumPy version in CI to catch version-related regressions early.
A minimal GitHub Actions workflow for running your tests might look like:
name: tests
on: [push, pull_request]
jobs:
test:
runs-on: ubuntu-latest
strategy:
matrix:
python-version: ["3.10", "3.11", "3.12"]
steps:
- uses: actions/checkout@v4
- uses: actions/setup-python@v5
with:
python-version: ${{ matrix.python-version }}
- run: pip install numpy pytest hypothesis pytest-benchmark
- run: pytest tests/ -v
Conclusion
Testing NumPy applications requires a mindset shift from traditional software testing. Floating-point tolerance, broadcasting behavior, shape preservation, and numerical stability all demand specialized assertions and strategies. By combining numpy.testing for precise array comparisons, Hypothesis for property-based exploration of edge cases, integration tests for pipeline correctness, and benchmark tests for performance regressions, you build a safety net that catches the subtle bugs unique to numerical computing. The upfront investment in a thorough test suite pays dividends every time you refactor, upgrade NumPy, or onboard a new contributor — your tests become the executable specification of how your numerical code should behave under every condition that matters.