Can good quality software improve my scientific career?
In many computational and applied sciences, software can be seen just as the means to reach a certain result. It derives that, from a merely technical aspect, several researchers — hopefully not computer scientists — are much more interested in making their codes work, rather than making them good. Of course, the reason is understandable: researchers have to do research, not coding. But is it really true that not caring about the quality of scientific scripts is the smartest or quickest way to do research?
We dedicate this article to providing a few considerations about the benefits of scientific code that respects (some of) the main concepts of software engineering, not only from a technical point of view but also from an academic perspective. Yes, right to convince even those who just don’t care about the code. Maybe.
“My code runs without any error on my laptop, it gives me the number I’m expecting for my experiment, why the hell should I lose time to improve it?”
Lots of scientific scripts are mostly a translation of mathematical operations. Excluding the big, well-structured projects (luckily, there are plenty of them!) the large variety of scientific programs do not care about the code quality, probably discouraged by the required effort and time. But are you sure the initial investment will not produce extra benefits in a long-term scenario? Maybe it’s better with some examples:
“I’ve used this script for half of my Ph.D. program for collecting the results of my thesis, but after updating the system I get different numbers”
“During the test of new features in my algorithm, I’m not able to make it work nor restore the previous version”
“I have to run again some tests because of the review of a paper, but I don’t even remember how to change those parameters”
Three situations that are not unrealistic. They of course are not referring to critical issues, but present three contexts where actually the quality of software could have prevented days of debugging. We introduce here some suggestions regarding the code style that should be used in scientific scripts, a very simple introduction to testing, and a few words about documentation and its great impact on scientific code.
Software engineering is a very wide aspect, so we just focus on some aspects that in my opinion are to take into account, especially with scientific software. For a more detailed overview, we refer to [1, 2].
Code Style
On the practical side, we list here a few rules to take into account during the development of scientific code. Naturally, there are plenty of guidelines to consider during the implementation, but instead of providing dozens of bullets here we just try to highlight the ones more important and easier to really follow on ordinary development.
- Use functions: rather than thinking of the source code as a (long) list of instructions, try to semantically split them into minimal (and possibly hierarchical) pieces of code.
- Keep it simple: as already mentioned, a scientific code should be easy to understand in order to enable future methodological enhancements.
- Source code should be self-informative: we suggest using meaningful names for variables, constants, functions, classes,
such that reading the source code allows in principle to understand also the algorithm. - Follow the standards: use the code standards to uniform all the code and, consequently, improve its readability.
As a simple example, we show here a few lines of code that are voluntarily written with very poor quality. Despite the number of lines, it’s extremely difficult to understand the meaning of that code (if you are not so familiar with the algorithm itself), resulting even in a harder task when we want to edit/use/test it.
import numpy as np
ss = []
mus = []
x = np.linspace(-1, 1, 36)
for i in range(20):
m = np.random.uniform()
s = m*(1./np.cosh(x+3))+(1-m)*(2./np.cosh(x)*np.tanh(x))
ss.append(s)
mus.append(m)
out = np.linalg.svd(np.asarray(ss).T)
myc = np.dot(out[0][:, :10].T, np.asarray(ss).T)
from scipy.interpolate import RBFInterpolator
i = RBFInterpolator(np.asarray(mus).reshape(-1, 1), myc.T)
print(np.linalg.norm(ss[0]-(out[0][:, :10]*i(np.asarray(mus[0]).reshape(-1, 1)).flatten()).sum(axis=1)))But let’s see how things change when we force ourselves to improve the quality. Instructions are re-ordered, we have isolated three main methods that are helping to recognize the semantic blocks in our algorithm, variables have now meaningful names, etc.
Now the code, excluding some reshape methods, is readable even by non-experts, improving a lot the reusability of that snippet!
import numpy as np
from scipy.interpolate import RBFInterpolator
def hf_func(x, mu):
f1 = 1./np.cosh(x+3)
f2 = 2./np.cosh(x)*np.tanh(x)
return mu*f1 + (1-mu)*f2
def offline(n, spatial_resolution=36, pod_rank=10):
""" Compute offline phase using `n` snapshots """
spatial_dof = np.linspace(-1, 1, spatial_resolution)
parameters = np.random.uniform(size=(n, 1))
snapshots = np.array([hf_func(spatial_dof, p) for p in parameters]).T
modes = np.linalg.svd(snapshots)[0][:, :pod_rank]
modal_coeffs = np.dot(modes.T, snapshots).T
interp = RBFInterpolator(parameters, modal_coeffs)
return parameters, snapshots, modes, interp
def train_error(parameters, snapshots, modes, interp, idx):
""" Compute the absolute l2 error obtained during online step at the `idx`
input parameter """
approx_modal_coeff = interp(parameters[idx].reshape(-1, 1)).flatten()
approx_sol = (modes * approx_modal_coeff).sum(axis=1)
return np.linalg.norm(snapshots[:, idx]-approx_sol)
print(train_error(*offline(20), idx=0))Testing
In producing data for scientific publications, the reliability of the software is of course another important feature. Things can become very frustrating discovering that your novel implementation returns awesome results for a specific problem, but it is not able to deal anymore with well-known benchmarks.
We need a way to test in a moment all the existing functionalities of the software. But how concretely do it? Let’s assume you are working on some code, which it’s already structured in several functions. We take as an example just a few lines of Python code that implement a simple norm computation.
def my_super_fancy_norm(vector):
return sqrt(sum([component**2 for component in vector]))Simple, isn’t it? Now, we have just to create another method (in a separate file), which basically check that, provided a specific input, our implementation return the expected output. We can simply do something like:
def test_my_norm():
input_ = numpy.array([1., 3., 5.4])
expected_result = 6.2577951388648065 # Manually computed
# You can also use other (well-tested) libraries to validate the result
#### expected_result = numpy.linalg.norm(input_)
numpy.testing.assert_equal(
my_super_fancy_norm(input_),
expected_result)In this way, after changing the my_super_fancy_norm code, we are quickly able to test its correctness, avoiding days and days of debugging.
There are also utilities like pytest that help in structuring and managing all the tests suite.
Documentation
In the scientific community, having (and maintaining!) informative documentation is a mandatory ingredient not only to provide the technical aspects of the software but also as the access point for new users to understand the algorithmic part. Practically, the documentation is a manual that illustrates all the details of the source code, like the expected input and output of any functions, as well as a guide for code employment for final users. But can I prepare that document without losing a lot of time?
The most efficient choice in order to create the documentation is usually by extracting it from the comments of the code, the so-called docstring. With the aim of illustrating how any function works, such comments are usually placed near functions prototypes: in this way, the documentation can be manipulated by simply editing the source code, being so a less impacting activity. A concrete example could be:
# from mathLab/PINA/pina/label_tensor.py
def extract(self, label_to_extract):
"""
Extract the subset of the original tensor by returning all the columns
corresponding to the passed `label_to_extract`.
:param label_to_extract: the label(s) to extract.
:type label_to_extract: str or iterable(str)
:raises TypeError: labels are not str
:raises ValueError: label to extract is not in the labels list
"""
...Several frameworks are available, depending on the used language, for this operation, giving the developers also the capability to aesthetic changes. Among the frameworks, we cite Sphinx and Doxygen.
References
[1] Wilson G, Aruliah DA, Brown CT, Chue Hong NP, Davis M, Guy RT, et al., Best practices for scientific computing, 2014, PLoS biology, doi:10.1371/journal.pbio.1001745
[2] N. Demo, M. Tezzele, G. Stabile, G. Rozza, Chapter 19: Scientific Software Development and Packages for Reduced Order Models in Computational Fluid Dynamics, 2022, Advanced Reduced Order Methods and Applications in Computational Fluid Dynamics, SIAM press Philadelphia, 2022, CS series Vol.27; doi:10.1137/1.9781611977257.ch19
[3] I. Sommerville, Software engineering: Seventh Edition, 2004, Pearson Education
