Développement de logiciels scientifiques

Software Development for Science

Partie 2: Aspects Pratiques avec Python

Dans cet article, nous suivrons les principes du TDD pour développer un logiciel scientifique, tel que présenté dans la première partie de cette série, afin de développer un filtre de détection de contours connu sous le nom de filtre de Sobel.

Dans le premier article, nous avons parlé de l’importance – et de la complexité – de développer des méthodes de test fiables pour les logiciels qui résolvent souvent des problèmes complexes dans le domaine scientifique. Nous avons également vu comment surmonter ces problèmes en suivant un cycle de développement inspiré du TDD, mais adapté au calcul scientifique. Je reproduis ci-dessous une version abrégée de ces instructions.

Cycle d’implémentation

1. Rassembler les exigences
2. Esquisser la conception
3. Mettre en place les tests initiaux
4. Implémenter votre version alpha
5. Construire une bibliothèque oracle
6. Revoir les tests
7. Implémenter le profilage

Cycle d’optimisation

1. Optimiser
2. Réimplémenter

Cycle de nouvelle méthode

1. Implémenter une nouvelle méthode
2. Valider par rapport aux oracles précédemment validés

Pour commencer: Le filtre de Sobel

Dans cet article, nous suivrons les instructions ci-dessus pour développer une fonction qui applique le filtre de Sobel. Le filtre de Sobel est un outil couramment utilisé en vision par ordinateur pour détecter ou améliorer les contours dans les images. Continuez à lire pour voir quelques exemples !

En commençant par la première étape, nous rassemblons quelques exigences. Nous suivrons la formulation standard du filtre de Sobel décrite dans cet article. En résumé, l’opérateur de Sobel consiste à convoluer l’image A avec les deux matrices 3 × 3 suivantes, à élever au carré chaque pixel des images résultantes, à les additionner et à prendre la racine carrée point par point. Si Ax et Ay sont les images résultant des convolutions, le résultat du filtre de Sobel S est √(Ax² + Ay²).

• L’informatique quantique pourrait bénéficier de la découverte du Q-Silicon
• Transistors en bois prennent racine
• L’Administration Biden envisage de renforcer les restrictions sur la vente de puces d’intelligence artificielle à la Chine.

Exigences

Nous voulons que cette fonction prenne un tableau 2D quelconque et génère un autre tableau 2D. Nous voulons peut-être qu’elle fonctionne sur n’importe quelles deux axes d’un ndarray. Nous voulons peut-être même qu’elle fonctionne sur plus (ou moins) de deux axes. Nous pouvons avoir des spécifications sur la manière de gérer les bords du tableau.

Pour l’instant, rappelons-nous de rester simple et de commencer par une implémentation en 2D. Mais avant de faire cela, esquissons la conception.

Esquisser la conception

Nous commençons par une conception simple, en utilisant les annotations de Python. Je recommande vivement d’annoter autant que possible, et d’utiliser mypy comme linter.

from typing import Tuplefrom numpy.core.multiarray import normalize_axis_indexfrom numpy.typing import NDArraydef sobel(arr: NDArray, axes: Tuple[int, int] = (-2, -1)) -> NDArray:    """Calculer le filtre de Sobel d'une image    Paramètres    ----------    arr : NDArray        Image d'entrée    axes : Tuple[int, int], optionnel        Axes sur lesquels calculer le filtre, par défaut (-2, -1)    Renvoie    -------    NDArray        Sortie    """    # N'accepte que des tableaux 2D    if arr.ndim != 2:        raise NotImplementedError    # Assurez-vous que l'axe[0] est le premier axe et l'axe[1] est le deuxième    # axe. L'obscure `normalize_axis_index` convertit les indices négatifs en    # indices entre 0 et arr.ndim - 1.    if any(        normalize_axis_index(ax, arr.ndim) != i        for i, ax in zip(range(2), axes)    ):        raise NotImplementedError    pass

Implémentez les Tests

Cette fonction ne fait pas grand-chose. Mais elle est documentée, annotée et ses limitations actuelles sont intégrées. Maintenant que nous avons une conception, nous passons immédiatement aux tests.

Tout d’abord, nous remarquons que les images vides (toutes les valeurs sont zéro) n’ont pas de contours. Elles doivent donc également retourner des zéros. En fait, toute image constante doit également retourner des zéros. Intégrons cela dans nos premiers tests. Nous verrons également comment nous pouvons utiliser le monkey testing pour tester de nombreux tableaux.

# test_zero_constant.pyimport numpy as npimport pytest# Testez plusieurs types de données en même [email protected](    "dtype",    ["float16", "float32", "float64", "float128"],)def test_zero(dtype):    # Définissez la graine aléatoire    seed = int(np.random.rand() * (2**32 - 1))    np.random.seed(seed)    # Créez un tableau 2D de forme aléatoire et remplissez-le de zéros    nx, ny = np.random.randint(3, 100, size=(2,))    arr = np.zeros((nx, ny), dtype=dtype)    # Appliquez la fonction sobel    arr_sob = sobel(arr)    # `assert_array_equal` échouera la plupart du temps.    # Cela ne fonctionnera que lorsque `arr_sob` est identiquement zéro,    # ce qui n'est généralement pas le cas.    # NE PAS UTILISER!    # np.testing.assert_array_equal(    #     arr_sob, 0.0, err_msg=f"{seed=} {nx=}, {ny=}"    # )    # `assert_almost_equal` peut échouer lorsqu'il est utilisé avec un grand nombre de décimales.    # Il repose également sur la vérification float64, qui peut échouer pour    # les types float128.    # NE PAS UTILISER!    # np.testing.assert_almost_equal(    #     arr_sob,    #     np.zeros_like(arr),    #     err_msg=f"{seed=} {nx=}, {ny=}",    #     decimal=4,    # )    # `assert_allclose` avec une tolérance personnalisée est ma méthode préférée    # Le 10 est arbitraire et dépend du problème. Si une méthode    # que vous savez être correcte ne passe pas, augmentez à 100, etc.    # Si la tolérance nécessaire pour réussir les tests est trop élevée,    # assurez-vous que la méthode est effectivement correcte.    tol = 10 * np.finfo(arr.dtype).eps    err_msg = f"{seed=} {nx=}, {ny=} {tol=}"  # Enregistrez les graines et autres informations    np.testing.assert_allclose(        arr_sob,        np.zeros_like(arr),        err_msg=err_msg,        atol=tol,  # rtol est inutile pour 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, size=(2,))    constant = np.random.randn(1).item()    arr = np.full((nx, ny), fill_value=constant, 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 est inutile pour desired=zeros    )

Ce fragment de code peut être exécuté depuis la ligne de commande avec

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

Version Alpha

Évidemment, nos tests échoueront actuellement, mais ils sont suffisants pour l’instant. Implémentons une première version.

from typing import Tupleimport numpy as npfrom numpy.core.multiarray import normalize_axis_indexfrom numpy.typing import NDArraydef sobel(arr: NDArray, axes: Tuple[int, int] = (-2, -1)) -> NDArray:    if arr.ndim != 2:        raise NotImplementedError    if any(        normalize_axis_index(ax, arr.ndim) != i        for i, ax in zip(range(2), axes)    ):        raise NotImplementedError    # Définissez nos noyaux de filtre. Remarquez qu'ils héritent du type d'entrée, donc    # un input float32 n'a jamais besoin d'être converti en float64 pour les calculs.    # Mais pouvez-vous voir où l'utilisation d'un autre dtype pour Gx et Gy pourrait être    # logique pour certains types de données d'entrée?    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,    )    # Créez le tableau de sortie et remplissez-le de zéros    s = np.zeros_like(arr)    for ix in range(1, arr.shape[0] - 1):        for iy in range(1, arr.shape[1] - 1):            # Multiplication point par point suivie d'une somme, aussi appelée 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

Avec cette nouvelle fonction, tous nos tests devraient passer et nous devrions obtenir une sortie comme celle-ci :

$ pytest -qq -s -x -vv --durations=0 test_zero_constant.py........======================================== durées les plus lentes =========================================0.09s appel     t_049988eae7f94139a7067f142bf2852f.py::test_constant[float16]0.08s appel     t_049988eae7f94139a7067f142bf2852f.py::test_zero[float64]0.06s appel     t_049988eae7f94139a7067f142bf2852f.py::test_constant[float128]0.04s appel     t_049988eae7f94139a7067f142bf2852f.py::test_zero[float128]0.04s appel     t_049988eae7f94139a7067f142bf2852f.py::test_constant[float64]0.02s appel     t_049988eae7f94139a7067f142bf2852f.py::test_constant[float32]0.01s appel     t_049988eae7f94139a7067f142bf2852f.py::test_zero[float16](17 durées < 0.005s cachées. Utilisez -vv pour afficher ces durées.)8 réussis en 0.35s

Nous avons jusqu’à présent :

1. Rassemblé les exigences de notre problème.
2. Esquissé une conception initiale.
3. Mis en place certains tests.
4. Implémenté la version alpha, qui réussit ces tests.

Ces tests offrent une vérification pour les versions futures, ainsi qu’une bibliothèque oracle très rudimentaire. Mais bien que les tests unitaires soient excellents pour capturer de petites variations par rapport aux résultats attendus, ils ne sont pas excellents pour différencier les résultats incorrects des résultats quantitativement différents mais néanmoins corrects. Supposons que demain nous voulions implémenter un autre opérateur de type Sobel, qui a un noyau plus long. Nous ne nous attendons pas à ce que les résultats correspondent exactement à notre version actuelle, mais nous nous attendons à ce que les sorties des deux fonctions soient qualitativement très similaires.

De plus, il est excellent de pratiquer en essayant de nombreuses entrées différentes pour nos fonctions afin de comprendre comment elles se comportent pour différentes entrées. Cela garantit que nous validons les résultats de manière scientifique.

Scikit-image dispose d’une excellente bibliothèque d’images que nous pouvons utiliser pour créer des oracles.

# !pip installscikit-image poochfrom typing import Dict, Callableimport numpy as npimport skimage.databwimages: Dict[str, np.ndarray] = {}for attrname in skimage.data.__all__:    attr = getattr(skimage.data, attrname)    # Les données sont obtenues via des appels de fonction    if isinstance(attr, Callable):        try:            # Télécharger les données            data = attr()            # Assurer qu'il s'agit d'un tableau 2D            if isinstance(data, np.ndarray) and data.ndim == 2:                # Convertir différents types entiers en float32 pour mieux                # évaluer la précision                bwimages[attrname] = data.astype(np.float32)        except:            continue# Appliquer le sobel aux imagesbwimages_sobel = {k: sobel(v) for k, v in bwimages.items()}

Une fois que nous avons les données, nous pouvons les représenter graphiquement.

import matplotlib.pyplot as pltfrom mpl_toolkits.axes_grid1 import make_axes_locatabledef create_colorbar(im, ax):    divider = make_axes_locatable(ax)    cax = divider.append_axes("right", size="5%", pad=0.1)    cb = ax.get_figure().colorbar(im, cax=cax, orientation="vertical")    return cax, cbfor name, data in bwimages.items():    fig, axs = plt.subplots(        1, 2, figsize=(10, 4), sharex=True, sharey=True    )    im = axs[0].imshow(data, aspect="equal", cmap="gray")    create_colorbar(im, axs[0])    axs[0].set(title=name)    im = axs[1].imshow(bwimages_sobel[name], aspect="equal", cmap="gray")    create_colorbar(im, axs[1])    axs[1].set(title=f"{name} sobel")

Ils ont l’air vraiment bien ! Je recommande de les stocker à la fois en tant que tableaux et en tant que figures que je peux rapidement vérifier pour une nouvelle version. Je recommande vivement HD5F pour le stockage des tableaux. Vous pouvez également configurer une galerie Sphinx pour générer directement les figures lors de la mise à jour de la documentation, de cette manière vous disposez d’une construction de figure reproductible que vous pouvez utiliser pour vérifier par rapport aux versions précédentes.

Une fois les résultats validés, vous pouvez les stocker sur le disque et les utiliser dans le cadre de vos tests unitaires. Quelque chose comme ceci :

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

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

Profilage

Le calcul du filtre Sobel pour ces jeux de données a pris un certain temps ! La prochaine étape consiste donc à profiler le code. Je recommande de créer un fichier benchmark_xyz.py pour chaque test, et de les exécuter régulièrement pour détecter les régressions de performances. Cela peut même faire partie de vos tests unitaires, mais nous n’irons pas si loin dans cet exemple. Une autre idée est d’utiliser les oracles ci-dessus pour les tests de performance.

Il existe de nombreuses façons de mesurer le temps d’exécution du code. Pour obtenir le temps écoulé réel au niveau du système, utilisez perf_counter_ns du module intégré time pour mesurer le temps en nanosecondes. Dans un notebook Jupyter, la magie de cellule %%timeit intégrée mesure l’exécution d’une certaine cellule. Le décorateur ci-dessous est inspiré de cette magie de cellule et peut être utilisé pour mesurer le temps d’exécution de n’importe quelle fonction.

import timefrom functools import wrapsfrom typing import Callable, Optionaldef 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: Optional[Callable] = None,    number: int = 10,) -> Callable:    """Appliquer à une fonction pour mesurer son exécution.    Paramètres    ----------    func_or_number : Optional[Callable], facultatif        Fonction à décorer ou argument `number` (voir ci-dessous), par        défaut None    number : int, facultatif        Nombre de fois que la fonction sera exécutée pour obtenir des statistiques, par        défaut 10    Renvoie    -------    Callable        Lorsqu'il est appliqué à une fonction, renvoie la fonction décorée. Sinon,        renvoie un décorateur.    Exemples    --------    .. code-block:: python        @timeit        def my_fun():            pass        @timeit(100)        def my_fun():            pass        @timeit(number=3)        def my_fun():            pass    """    if isinstance(func_or_number, Callable):        func = func_or_number        number = number    elif isinstance(func_or_number, int):        func = None        number = func_or_number    else:        func = None        number = number    def decorator(f):        @wraps(f)        def wrapper(*args, **kwargs):            runs_ns = np.empty((number,))            # Exécuter la première fois et mesurer le résultat            start_time = time.perf_counter_ns()            result = f(*args, **kwargs)            runs_ns[0] = time.perf_counter_ns() - start_time            for i in range(1, number):                start_time = time.perf_counter_ns()                f(*args, **kwargs)  # Sans stockage, plus rapide                runs_ns[i] = time.perf_counter_ns() - start_time            time_msg = f"{sizeof_fmt(runs_ns.mean())} ± "            time_msg += f"{sizeof_fmt(runs_ns.std())}"            print(                f"{time_msg} par exécution (moyenne ± écart-type de {number} exécutions)"            )            return result        return wrapper    if func is not None:        return decorator(func)    return decorator

Mettre notre fonction à l’épreuve:

arr_test = np.random.randn(500, 500)sobel_timed = timeit(sobel)sobel_timed(arr_test);# 3.9s ± 848.9 ms par exécution (moyenne ± écart type de 10 exécutions)

Pas très rapide 🙁

Lors de l’investigation de la lenteur ou des régressions de performance, il est également possible de suivre le temps d’exécution de chaque ligne. La bibliothèque line_profiler est une excellente ressource pour cela. Elle peut être utilisée via la cellule magique de Jupyter, ou en utilisant l’API. Voici un exemple d’utilisation de l’API:

from line_profiler import LineProfilerlp = LineProfiler()lp_wrapper = lp(sobel)lp_wrapper(arr_test)lp.print_stats(output_unit=1)  # 1 pour les secondes, 1e-3 pour les millisecondes, etc.

Voici un exemple de sortie:

Unité de temps: 1 sTemps total: 4.27197 sFichier: /tmp/ipykernel_521529/1313985340.pyFonction: sobel à la ligne 8Ligne #      Hits         Time  Par Hit   % Temps  Contenu de la ligne==============================================================     8                                           def sobel(arr: NDArray, axes: Tuple[int, int] = (-2, -1)) -> NDArray:     9                                               # Accepte uniquement des tableaux 2D    10         1          0.0      0.0      0.0      if arr.ndim != 2:    11                                                   raise NotImplementedError    12                                               13                                               # Assurez-vous que l'axe[0] est le premier axe et l'axe[1] est le deuxième axe    14                                               # L'obscure `normalize_axis_index` convertit les indices négatifs en    15                                               # indices entre 0 et 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(range(2), axes)    19                                               ):    20                                                   raise 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 range(1, arr.shape[0] - 1):    32    248004          0.1      0.0      2.2          for iy in range(1, arr.shape[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

Beaucoup de temps est passé à l’intérieur de la boucle la plus interne. NumPy préfère les mathématiques vectorisées, car il peut ensuite s’appuyer sur un code compilé. Étant donné que nous utilisons des boucles explicites, il est logique que cette fonction soit très lente.

De plus, il est important d’être intelligent quant aux allocations de mémoire à l’intérieur des boucles. NumPy est quelque peu intelligent lorsqu’il s’agit d’allouer de petites quantités de mémoire à l’intérieur des boucles, donc les lignes définissant s1 ou s2 pourraient ne pas allouer de nouvelle mémoire. Mais cela pourrait également être le cas, ou il pourrait y avoir une autre allocation de mémoire qui se produit en arrière-plan dont nous ne sommes pas conscients. Je recommande donc également de profiler la mémoire. J’aime utiliser Memray pour cela, mais il y en a d’autres comme Fil et Sciagraph. Je recommande également d’éviter memory_profiler qui (très malheureusement !) n’est plus maintenu. De plus, Memray est beaucoup plus puissant. Voici un exemple de Memray en action :

$ # Créez sobel.bin qui contient les informations de profilage$ memray run -fo sobel.bin --trace-python-allocators sobel_run.pyÉcriture des résultats de profilage dans sobel.binMemray AVERTISSEMENT : Correction du symbole pour aligned_alloc de 0x7fc5c984d8f0 à 0x7fc5ca4a5ce0[memray] Résultats de profilage générés avec succès.Vous pouvez maintenant générer des rapports à partir des enregistrements d'allocation stockés.Quelques commandes d'exemple pour générer des rapports :python3 -m memray flamegraph sobel.bin

$ # Générer un graphique en flammes$ memray flamegraph -fo sobel_flamegraph.html --temporary-allocations sobel.bin⠙ Calcul de la ligne de flottaison maximale... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸  99% 0:00:0100:01⠏ Traitement des enregistrements d'allocation... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸  98% 0:00:0100:01Fichier sobel_flamegraph.html écrit

$ # Afficher l'arbre de mémoire$ memray tree --temporary-allocations sobel.bin⠧ Calcul de la ligne de flottaison maximale... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸ 100% 0:00:0100:01⠧ Traitement des enregistrements d'allocation... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸ 100% 0:00:0100:01Métadonnées d'allocation-------------------Arguments de ligne de commande : 'memray run -fo sobel.bin --trace-python-allocators sobel_run.py'Taille maximale de la mémoire : 11.719 MoNombre d'allocations : 15332714Les 10 plus grandes allocations :-----------------------📂 123.755 Mo (100.00 %) <ROOT>  └── [[3 frames cachés dans 2 fichiers]]    └── 📂 123.755 Mo (100.00 %) _run_code  /usr/lib/python3.10/runpy.py:86        ├── 📂 122.988 Mo (99.38 %) <module>  sobel_run.py:40        │   ├── 📄 51.087 Mo (41.28 %) sobel  sobel_run.py:35        │   ├── [[1 frame caché dans 1 fichier]]        │   │   └── 📄 18.922 Mo (15.29 %) _sum          │   │       lib/python3.10/site-packages/numpy/core/_methods.py:49        │   └── [[1 frame caché dans 1 fichier]]        │       └── 📄 18.921 Mo (15.29 %) _sum          │           lib/python3.10/site-packages/numpy/core/_methods.py:49...

Version bêta et un bogue

Maintenant que nous avons une version alpha fonctionnelle et quelques fonctions de profilage, nous utiliserons la bibliothèque SciPy pour obtenir une version beaucoup plus rapide.

from typing import Tupleimport numpy as npfrom numpy.core.multiarray import normalize_axis_indexfrom numpy.typing import NDArrayfrom scipy.signal import convolve2ddef sobel_conv2d(    arr: NDArray, axes: Tuple[int, int] = (-2, -1)) -> NDArray:    if arr.ndim != 2:        raise NotImplementedError    if any(        normalize_axis_index(ax, arr.ndim) != i        for i, ax in zip(range(2), axes)    ):        raise NotImplementedError        # Créez les noyaux en tant qu'ensemble unique et complexe. Cela nous permet d'utiliser    # np.abs au lieu de np.hypot pour calculer l'amplitude.    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="same")    np.absolute(s, out=s)  # Absolu en place    return s.real

sobel_timed = timeit(sobel_conv2d)sobel_timed(arr_test)# 14,3 ms ± 1,71 ms par boucle (moyenne ± écart type de 10 exécutions)

Beaucoup mieux ! Mais est-ce correct ?

Les images semblent très similaires, mais si vous observez l’échelle des couleurs, elles ne sont pas identiques. Les tests indiquent une petite erreur moyenne. Heureusement, nous sommes maintenant bien équipés pour détecter les différences quantitatives et qualitatives.

Après avoir enquêté sur ce bogue, nous l’attribuons aux différentes conditions aux limites. En consultant la documentation de convolve2d, nous apprenons que le tableau d’entrée est rembourré de zéros avant d’appliquer le noyau. Dans la version alpha, nous avons rembourré la sortie !

Laquelle est correcte ? On pourrait soutenir que l’implémentation de SciPy a plus de sens. Dans ce cas, nous devrions adopter la nouvelle version comme oracle, et si nécessaire, corriger la version alpha pour qu’elle corresponde à celle-ci. C’est courant dans le développement de logiciels scientifiques : de nouvelles informations sur la meilleure façon de faire changent les oracles et les tests.

Dans ce cas, la correction est simple, rembourrons le tableau de zéros avant de le traiter.

def sobel_v2(arr: NDArray, axes: Tuple[int, int] = (-2, -1)) -> NDArray:    # ...    arr = np.pad(arr, (1,))  # Après le rembourrage, il a une forme (nx + 2, ny + 2)    s = np.zeros_like(arr)    for ix in range(1, arr.shape[0] - 1):        for iy in range(1, arr.shape[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)  # Ajuster les indices    return s

Une fois que nous avons corrigé notre fonction, nous pouvons mettre à jour les oracles et les tests qui en dépendent.

Réflexions finales

Nous avons vu comment mettre en pratique certaines des idées de développement de logiciels explorées dans cet article. Nous avons également vu certains outils que vous pouvez utiliser pour développer du code de haute qualité et performant.

Je vous suggère d’essayer certaines de ces idées sur vos propres projets. En particulier, pratiquez le profilage du code et son amélioration. La fonction de filtre Sobel que nous avons codée est très inefficace, je vous suggère d’essayer de l’améliorer.

Par exemple, essayez la parallélisation CPU avec un compilateur JIT tel que Numba, portez la boucle interne en Cython, ou implémentez une fonction GPU CUDA avec Numba ou CuPy. N’hésitez pas à consulter ma série sur la programmation de noyaux CUDA avec Numba.

We will continue to update IPGirl; if you have any questions or suggestions, please contact us!