Developing Scientific Software. Part 2: Practical Aspects with Python | by Carlos Costa, Ph.D. | Jul, 2023

Half 2: Sensible Points with Python

On this article we’ll observe the tenets of TDD for growing Scientific Software program as specified by the first installment of this series to develop an edge detection filter often known as the Sobel filter.

Within the first article, we talked about how vital — and difficult — it may be to develop dependable testing strategies for software program which the advanced issues usually present in scientific software program. We additionally noticed the right way to overcome these points by following a growth cycle impressed by TDD, however tailored for scientific computing. I reproduce a shortened model of those directions under.

Implementation cycle

1. Collect necessities
2. Sketch the design
3. Implement preliminary exams
4. Implement your alpha model
5. Construct an oracle library
6. Revisit exams
7. Implement profiling

Optimization cycle

1. Optimize
2. Reimplement

New technique cycle

1. Implement new technique
2. Validate towards earlier curated oracles

On this article, we’ll observe the above directions to develop a perform which applies the Sobel filter. The Sobel filter is a generally used pc imaginative and prescient software to detect or improve edges in pictures. Preserve studying to see some examples!

Beginning with step one, we collect some necessities. We are going to observe the usual formulation of the Sobel filter described in this article. Merely put, the Sobel operator consists of convolving picture A with the next two 3 × 3 kernels, squaring every pixel within the ensuing pictures, summing them and taking the point-by-point sq. root. If Ax and Ay are the photographs ensuing from the convolutions, the results of the Sobel filter S is √(Ax² + Ay²).

We wish this perform to take any 2D array and generate one other 2D array. We’d need it to function on any two axes of an ndarray. We’d even need it to work on extra (or much less) than two axes. We’d have specs for the right way to deal with the perimeters of the array.

For now let’s keep in mind to maintain it easy, and begin with a 2D implementation. However earlier than we do this, let’s sketch the design.

We begin with a easy design, leveraging Python’s annotations. I extremely advocate annotating as a lot as doable, and utilizing mypy as a linter.

from typing import Tuple
from numpy.core.multiarray import normalize_axis_index
from numpy.typing import NDArray
def sobel(arr: NDArray, axes: Tuple[int, int] = (-2, -1)) -> NDArray:
"""Compute the Sobel filter of a picture
Parameters
----------
arr : NDArray
Enter picture
axes : Tuple[int, int], optionally available
Axes over which to compute the filter, by default (-2, -1)
Returns
-------
NDArray
Output
"""
# Solely accepts 2D arrays
if arr.ndim != 2:
elevate NotImplementedError
# Make sure that the axis[0] is the primary axis, and axis[1] is the second
# axis. The obscure `normalize_axis_index` converts unfavorable indices to
# indices between 0 and arr.ndim - 1.
if any(
normalize_axis_index(ax, arr.ndim) != i
for i, ax in zip(vary(2), axes)
):
elevate NotImplementedError
cross

This perform doesn’t do a lot. However it’s documented, annotated, and its present limitations are baked into it. Now that we’ve a design, we instantly pivot to exams.

First, we discover that vacant pictures (all zeroes) don’t have any edges. So they need to output zeroes as properly. In reality, any fixed picture must also return zeros. Let’s bake that into out first exams. We may also see how we will use monkey testing to check many arrays.

# test_zero_constant.py
import numpy as np
import pytest
# Take a look at a number of dtypes directly
@pytest.mark.parametrize(
"dtype",
["float16", "float32", "float64", "float128"],
)
def test_zero(dtype):
# Set random seed
seed = int(np.random.rand() * (2**32 - 1))
np.random.seed(seed)
# Create a 2D array of random form and fill with zeros
nx, ny = np.random.randint(3, 100, measurement=(2,))
arr = np.zeros((nx, ny), dtype=dtype)
# Apply sobel perform
arr_sob = sobel(arr)
# `assert_array_equal` ought to fail many of the occasions.
# It can solely work when `arr_sob` is identically zero,
# which is often not the case.
# DO NOT USE!
# np.testing.assert_array_equal(
#     arr_sob, 0.0, err_msg=f"{seed=} {nx=}, {ny=}"
# )
# `assert_almost_equal` can fail when used with excessive decimals.
# It additionally depends on float64 checking, which could fail for
# float128 sorts.
# DO NOT USE!
# np.testing.assert_almost_equal(
#     arr_sob,
#     np.zeros_like(arr),
#     err_msg=f"{seed=} {nx=}, {ny=}",
#     decimal=4,
# )
# `assert_allclose` with customized tolerance is my most popular technique
# The ten is bigoted and is dependent upon the issue. If a technique
# which you recognize to be appropriate doesn't cross, enhance to 100, and so forth.
# If the tolerance wanted to make the exams cross is simply too excessive, make
# positive the strategy is definitely appropriate.
tol = 10 * np.finfo(arr.dtype).eps
err_msg = f"{seed=} {nx=}, {ny=} {tol=}"  # Log seeds and different information
np.testing.assert_allclose(
arr_sob,
np.zeros_like(arr),
err_msg=err_msg,
atol=tol,  # rtol is ineffective for desired=zeros
)
@pytest.mark.parametrize(
"dtype", ["float16", "float32", "float64", "float128"]
)
def test_constant(dtype):
seed = int(np.random.rand() * (2**32 - 1))
np.random.seed(seed)
nx, ny = np.random.randint(3, 100, measurement=(2,))
fixed = np.random.randn(1).merchandise()
arr = np.full((nx, ny), fill_value=fixed, dtype=dtype)
arr_sob = sobel(arr)
tol = 10 * np.finfo(arr.dtype).eps
err_msg = f"{seed=} {nx=}, {ny=} {tol=}"
np.testing.assert_allclose(
arr_sob,
np.zeros_like(arr),
err_msg=err_msg,
atol=tol,  # rtol is ineffective for desired=zeros
)

This code snippet could be run from the command-line with

$ pytest -qq -s -x -vv --durations=0 test_zero_constant.py

In fact our exams will at present fail, however they’re sufficient for now. Let’s implement a primary model.

from typing import Tuple
import numpy as np
from numpy.core.multiarray import normalize_axis_index
from numpy.typing import NDArray
def sobel(arr: NDArray, axes: Tuple[int, int] = (-2, -1)) -> NDArray:
if arr.ndim != 2:
elevate NotImplementedError
if any(
normalize_axis_index(ax, arr.ndim) != i
for i, ax in zip(vary(2), axes)
):
elevate NotImplementedError
# Outline our filter kernels. Discover they inherit the enter kind, so
# {that a} float32 enter by no means must be solid to float64 for computation.
# However are you able to see the place utilizing one other dtype for Gx and Gy would possibly make
# sense for some enter dtypes?
Gx = np.array(
[[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]],
dtype=arr.dtype,
)
Gy = np.array(
[[-1, -2, -1], [0, 0, 0], [1, 2, 1]],
dtype=arr.dtype,
)
# Create the output array and fill with zeroes
s = np.zeros_like(arr)
for ix in vary(1, arr.form[0] - 1):
for iy in vary(1, arr.form[1] - 1):
# Pointwise multiplication adopted by sum, aka convolution
s1 = (Gx * arr[ix - 1 : ix + 2, iy - 1 : iy + 2]).sum()
s2 = (Gy * arr[ix - 1 : ix + 2, iy - 1 : iy + 2]).sum()
s[ix, iy] = np.hypot(s1, s2)  # np.sqrt(s1**2 + s2**2)
return s

With this new perform, all our exams ought to cross, and we must always get an output like this:

$ pytest -qq -s -x -vv --durations=0 test_zero_constant.py
........
======================================== slowest durations =========================================
0.09s name     t_049988eae7f94139a7067f142bf2852f.py::test_constant[float16]
0.08s name     t_049988eae7f94139a7067f142bf2852f.py::test_zero[float64]
0.06s name     t_049988eae7f94139a7067f142bf2852f.py::test_constant[float128]
0.04s name     t_049988eae7f94139a7067f142bf2852f.py::test_zero[float128]
0.04s name     t_049988eae7f94139a7067f142bf2852f.py::test_constant[float64]
0.02s name     t_049988eae7f94139a7067f142bf2852f.py::test_constant[float32]
0.01s name     t_049988eae7f94139a7067f142bf2852f.py::test_zero[float16]
(17 durations < 0.005s hidden.  Use -vv to point out these durations.)
8 handed in 0.35s

We’ve got thus far:

1. Gathered necessities of our drawback.
2. Sketched an preliminary design.
3. Carried out some exams.
4. Carried out the alpha model, which passes these exams.

These exams present verification for future variations, in addition to a really rudimentary oracle library. However whereas unit exams are nice at capturing minute deviations from anticipated outcomes, they don’t seem to be nice at differentiating fallacious outcomes from quantitatively totally different — however nonetheless appropriate — outcomes. Suppose tomorrow we wish to implement one other Sobel-type operator, which has an extended kernel. We don’t anticipate the outcomes to match precisely to our present model, however we do anticipate the outputs of each capabilities to be qualitatively very comparable.

As well as, it’s glorious observe to strive many various inputs to our capabilities to get a way of how they behave for various inputs. This ensures that we scientifically validate the outcomes.

Scikit-image has a wonderful library of pictures, which we will use to create oracles.

# !pip installscikit-image pooch
from typing import Dict, Callable
import numpy as np
import skimage.knowledge
bwimages: Dict[str, np.ndarray] = {}
for attrname in skimage.knowledge.__all__:
attr = getattr(skimage.knowledge, attrname)
# Information are obtained by way of perform calls
if isinstance(attr, Callable):
strive:
# Obtain the information
knowledge = attr()
# Guarantee it's a 2D array
if isinstance(knowledge, np.ndarray) and knowledge.ndim == 2:
# Convert from varied int sorts to float32 to higher
# assess precision
bwimages[attrname] = knowledge.astype(np.float32)
besides:
proceed
# Apply sobel to pictures
bwimages_sobel = {okay: sobel(v) for okay, v in bwimages.gadgets()}

As soon as we’ve the information, we will plot it.

import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
def create_colorbar(im, ax):
divider = make_axes_locatable(ax)
cax = divider.append_axes("proper", measurement="5%", pad=0.1)
cb = ax.get_figure().colorbar(im, cax=cax, orientation="vertical")
return cax, cb
for identify, knowledge in bwimages.gadgets():
fig, axs = plt.subplots(
1, 2, figsize=(10, 4), sharex=True, sharey=True
)
im = axs[0].imshow(knowledge, side="equal", cmap="grey")
create_colorbar(im, axs[0])
axs[0].set(title=identify)
im = axs[1].imshow(bwimages_sobel[name], side="equal", cmap="grey")
create_colorbar(im, axs[1])
axs[1].set(title=f"{identify} sobel")

These look actually good! I might advocate storing these, each as arrays and as figures which I can shortly test towards for a brand new model. I extremely advocate HD5F for array storage. You can even arrange a Sphinx Gallery to instantly generate the figures them when updating documentation, that approach you will have a reproducible determine construct that you should use to test towards earlier variations.

After the outcomes have been validated, you may retailer them on disk and use them as a part of your unit testing. One thing like this:

oracle_library = [(k, v, bwimages_sobel[k]) for okay, v in bwimages.gadgets()]
# save_to_disk(oracle_library, ...)

# test_oracle.py
import numpy as np
import pytest
from numpy.typing import NDArray
# oracle_library = read_from_disk(...)
@pytest.mark.parametrize("identify,enter,output", oracle_library)
def test_oracles(identify: str, enter: NDArray, output: NDArray):
output_new = sobel(enter)
tol = 10 * np.finfo(enter.dtype).eps
mean_avg_error = np.abs(output_new - output).imply()
np.testing.assert_allclose(
output_new,
output,
err_msg=f"{identify=} {tol=} {mean_avg_error=}",
atol=tol,
rtol=tol,
)

Computing the Sobel filter for these datasets took some time! So the subsequent step is to profile the code. I like to recommend making a benchmark_xyz.py file for every check, and re-run them repeatedly to probe for efficiency regressions. This could even be a part of your unit testing, however we gained’t go thus far on this instance. One other concept is to make use of the oracles above for benchmarking.

There are numerous methods of timing code execution. To get the system-wide, “actual” elapsed time, use perf_counter_ns from the built-in time module to measure time in nanoseconds. In a Jupyter pocket book, the built-in %%timeit cell magic occasions a sure cell execution. The decorator under is impressed by this cell magic and can be utilized to time any perform.

import time
from functools import wraps
from typing import Callable, Non-compulsory
def sizeof_fmt(num, suffix="s"):
for unit in ["n", "μ", "m"]:
if abs(num) < 1000:
return f"{num:3.1f} {unit}{suffix}"
num /= 1000
return f"{num:.1f}{suffix}"
def timeit(
func_or_number: Non-compulsory[Callable] = None,
quantity: int = 10,
) -> Callable:
"""Apply to a perform to time its executions.
Parameters
----------
func_or_number : Non-compulsory[Callable], optionally available
Operate to be embellished or `quantity` argument (see under), by
default None
quantity : int, optionally available
Variety of occasions the perform will run to acquire statistics, by
default 10
Returns
-------
Callable
When fed a perform, returns the embellished perform. In any other case
returns a decorator.
Examples
--------
.. code-block:: python
@timeit
def my_fun():
cross
@timeit(100)
def my_fun():
cross
@timeit(quantity=3)
def my_fun():
cross
"""
if isinstance(func_or_number, Callable):
func = func_or_number
quantity = quantity
elif isinstance(func_or_number, int):
func = None
quantity = func_or_number
else:
func = None
quantity = quantity
def decorator(f):
@wraps(f)
def wrapper(*args, **kwargs):
runs_ns = np.empty((quantity,))
# Run first and measure retailer the outcome
start_time = time.perf_counter_ns()
outcome = f(*args, **kwargs)
runs_ns[0] = time.perf_counter_ns() - start_time
for i in vary(1, quantity):
start_time = time.perf_counter_ns()
f(*args, **kwargs)  # With out storage, sooner
runs_ns[i] = time.perf_counter_ns() - start_time
time_msg = f"{sizeof_fmt(runs_ns.imply())} ± "
time_msg += f"{sizeof_fmt(runs_ns.std())}"
print(
f"{time_msg} per run (imply ± std. dev. of {quantity} runs)"
)
return outcome
return wrapper
if func will not be None:
return decorator(func)
return decorator

Placing our perform to the check:

arr_test = np.random.randn(500, 500)
sobel_timed = timeit(sobel)
sobel_timed(arr_test);
# 3.9s ± 848.9 ms per run (imply ± std. dev. of 10 runs)

Not too quick 🙁

When investigating slowness or efficiency regressions, it is usually doable to trace the runtime of every line. The line_profiler library is a superb useful resource for this. It may be used by way of Jupyter cell magic, or utilizing the API. Right here is an API instance:

from line_profiler import LineProfiler
lp = LineProfiler()
lp_wrapper = lp(sobel)
lp_wrapper(arr_test)
lp.print_stats(output_unit=1)  # 1 for seconds, 1e-3 for milliseconds, and so forth.

Right here is an instance output:

Timer unit: 1 s
Complete time: 4.27197 s
File: /tmp/ipykernel_521529/1313985340.py
Operate: sobel at line 8
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
8                                           def sobel(arr: NDArray, axes: Tuple[int, int] = (-2, -1)) -> NDArray:
9                                               # Solely accepts 2D arrays
10         1          0.0      0.0      0.0      if arr.ndim != 2:
11                                                   elevate NotImplementedError
12                                           
13                                               # Make sure that the axis[0] is the primary axis, and axis[1] is the second
14                                               # axis. The obscure `normalize_axis_index` converts unfavorable indices to
15                                               # indices between 0 and arr.ndim - 1.
16         1          0.0      0.0      0.0      if any(
17                                                   normalize_axis_index(ax, arr.ndim) != i
18         1          0.0      0.0      0.0          for i, ax in zip(vary(2), axes)
19                                               ):
20                                                   elevate NotImplementedError
21                                           
22         1          0.0      0.0      0.0      Gx = np.array(
23         1          0.0      0.0      0.0          [[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]],
24         1          0.0      0.0      0.0          dtype=arr.dtype,
25                                               )
26         1          0.0      0.0      0.0      Gy = np.array(
27         1          0.0      0.0      0.0          [[-1, -2, -1], [0, 0, 0], [1, 2, 1]],
28         1          0.0      0.0      0.0          dtype=arr.dtype,
29                                               )
30         1          0.0      0.0      0.0      s = np.zeros_like(arr)
31       498          0.0      0.0      0.0      for ix in vary(1, arr.form[0] - 1):
32    248004          0.1      0.0      2.2          for iy in vary(1, arr.form[1] - 1):
33    248004          1.8      0.0     41.5              s1 = (Gx * arr[ix - 1 : ix + 2, iy - 1 : iy + 2]).sum()
34    248004          1.7      0.0     39.5              s2 = (Gy * arr[ix - 1 : ix + 2, iy - 1 : iy + 2]).sum()
35    248004          0.7      0.0     16.8              s[ix, iy] = np.hypot(s1, s2)
36         1          0.0      0.0      0.0      return s

Lot’s of time is spend contained in the innermost loop. NumPy prefers vectorized math, as it may well then depend on compiled code. Since we’re utilizing specific for loops, it is smart that this perform could be very gradual.

As well as, you will need to be good about reminiscence allocations inside loops. NumPy is considerably good about allocating small quantities of reminiscence inside loops, so the traces defining s1 or s2 may not be allocating new reminiscence. However additionally they might be, or there might be another reminiscence allocation that’s occurring beneath the hood that we’re not conscious of. I due to this fact advocate additionally profiling reminiscence. I like utilizing Memray for that, however there are others like Fil and Sciagraph. I might additionally keep away from memory_profiler which (very sadly!) is not maintained. Additionally Memray is far more highly effective. Right here is an instance of Memray in motion:

$ # Create sobel.bin which holds the profiling data
$ memray run -fo sobel.bin --trace-python-allocators sobel_run.py
Writing profile outcomes into sobel.bin
Memray WARNING: Correcting image for aligned_alloc from 0x7fc5c984d8f0 to 0x7fc5ca4a5ce0
[memray] Efficiently generated profile outcomes.
Now you can generate studies from the saved allocation information.
Some instance instructions to generate studies:
python3 -m memray flamegraph sobel.bin

$ # Generate flame graph
$ memray flamegraph -fo sobel_flamegraph.html --temporary-allocations sobel.bin
⠙ Calculating excessive watermark... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸  99% 0:00:0100:01
⠏ Processing allocation information... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸  98% 0:00:0100:01
Wrote sobel_flamegraph.html

$ # Present reminiscence tree
$ memray tree --temporary-allocations sobel.bin
⠧ Calculating excessive watermark... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸ 100% 0:00:0100:01
⠧ Processing allocation information... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸ 100% 0:00:0100:01
Allocation metadata
-------------------
Command line arguments: 
'memray run -fo sobel.bin --trace-python-allocators sobel_run.py'
Peak reminiscence measurement: 11.719MB
Variety of allocations: 15332714
Largest 10 allocations:
-----------------------
📂 123.755MB (100.00 %) <ROOT>  
└── [[3 frames hidden in 2 file(s)]]
└── 📂 123.755MB (100.00 %) _run_code  /usr/lib/python3.10/runpy.py:86
├── 📂 122.988MB (99.38 %) <module>  sobel_run.py:40
│   ├── 📄 51.087MB (41.28 %) sobel  sobel_run.py:35
│   ├── [[1 frames hidden in 1 file(s)]]
│   │   └── 📄 18.922MB (15.29 %) _sum  
│   │       lib/python3.10/site-packages/numpy/core/_methods.py:49
│   └── [[1 frames hidden in 1 file(s)]]
│       └── 📄 18.921MB (15.29 %) _sum  
│           lib/python3.10/site-packages/numpy/core/_methods.py:49
...

Now that we’ve a working alpha and a few profiling capabilities, we’ll leverage the SciPy library to acquire a a lot sooner model.

from typing import Tuple
import numpy as np
from numpy.core.multiarray import normalize_axis_index
from numpy.typing import NDArray
from scipy.sign import convolve2d
def sobel_conv2d(
arr: NDArray, axes: Tuple[int, int] = (-2, -1)
) -> NDArray:
if arr.ndim != 2:
elevate NotImplementedError
if any(
normalize_axis_index(ax, arr.ndim) != i
for i, ax in zip(vary(2), axes)
):
elevate NotImplementedError
# Create the kernels as a single, advanced array. Permits us to make use of
# np.abs as an alternative of np.hypot to calculate the magnitude.
G = np.array(
[[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]],
dtype=arr.dtype,
)
G = G + 1j * np.array(
[[-1, -2, -1], [0, 0, 0], [1, 2, 1]],
dtype=arr.dtype,
)
s = convolve2d(arr, G, mode="similar")
np.absolute(s, out=s)  # In-place abs
return s.actual

sobel_timed = timeit(sobel_conv2d)
sobel_timed(arr_test)
# 14.3 ms ± 1.71 ms per loop (imply ± std. dev. of 10 runs)

A lot better! However is it proper?

The photographs look very comparable, however should you discover the colour scale, they don’t seem to be the identical. Working the exams flags a small imply common error. Fortunately, we at the moment are well-equipped at detecting quantitative and qualitative variations.

After investigating this bug, we attribute it to the totally different boundary situations. Trying into convolve2d documentation tells us that the enter array is padded with zeroes earlier than making use of the kernel. Within the alpha, we padded the output!

Which one is true? Arguably the SciPy implementation makes extra sense. On this case we must always undertake the brand new model because the oracle, and if required, repair the alpha model to match it. That is widespread in scientific software program growth: new data of the right way to do issues higher adjustments the oracles and the exams.

On this case, the repair is simple, pad the array with zeros previous to processing it.

def sobel_v2(arr: NDArray, axes: Tuple[int, int] = (-2, -1)) -> NDArray:
# ...
arr = np.pad(arr, (1,))  # After padding, it's formed (nx + 2, ny + 2)
s = np.zeros_like(arr)
for ix in vary(1, arr.form[0] - 1):
for iy in vary(1, arr.form[1] - 1):
s1 = (Gx * arr[ix - 1 : ix + 2, iy - 1 : iy + 2]).sum()
s2 = (Gy * arr[ix - 1 : ix + 2, iy - 1 : iy + 2]).sum()
s[ix - 1, iy - 1] = np.hypot(s1, s2)  # Alter indices
return s

As soon as we corrected out perform, we will replace the oracles and exams which depend on them.

We noticed the right way to put in observe a number of the software program growth concepts explored in this article. We additionally noticed some instruments which you should use to develop high-quality, high-performance code.

I recommend you strive a few of these concepts by yourself tasks. Particularly, observe profiling code and enhancing it. The Sobel filter perform we coded could be very inefficient, I recommend attempting to enhance it.

For instance, strive CPU parallelization with a JIT compiler similar to Numba, porting the inside loop into Cython, or implementing a CUDA GPU perform with Numba or CuPy. Be at liberty to take a look at my series on coding CUDA kernels with Numba.