Desarrollo de Software Científico

Scientific Software Development.

Parte 2: Aspectos Prácticos con Python

En este artículo seguiremos los principios de TDD para desarrollar Software Científico, como se explica en la primera entrega de esta serie, para desarrollar un filtro de detección de bordes conocido como filtro Sobel.

En el primer artículo, hablamos de lo importante —y complicado— que puede ser desarrollar métodos de prueba confiables para software que aborda problemas complejos, como los que se encuentran comúnmente en el software científico. También vimos cómo superar esos problemas siguiendo un ciclo de desarrollo inspirado en TDD, pero adaptado para la computación científica. A continuación, reproduzco una versión abreviada de estas instrucciones.

Ciclo de implementación

1. Reunir los requisitos
2. Esbozar el diseño
3. Implementar pruebas iniciales
4. Implementar la versión alfa
5. Construir una biblioteca de oráculos
6. Revisar pruebas
7. Implementar el perfilado

Ciclo de optimización

1. Optimizar
2. Reimplementar

Ciclo de nuevo método

1. Implementar nuevo método
2. Validar contra oráculos curados previos

Empezando: El Filtro Sobel

En este artículo, seguiremos las instrucciones anteriores para desarrollar una función que aplique el filtro Sobel. El filtro Sobel es una herramienta comúnmente utilizada en visión por computadora para detectar o mejorar los bordes en imágenes. ¡Sigue leyendo para ver algunos ejemplos!

Comenzando con el primer paso, reunimos algunos requisitos. Seguiremos la formulación estándar del filtro Sobel descrita en este artículo. En pocas palabras, el operador Sobel consiste en convolucionar la imagen A con los siguientes dos kernels de 3 × 3, elevar al cuadrado cada píxel en las imágenes resultantes, sumarlos y tomar la raíz cuadrada punto por punto. Si Ax y Ay son las imágenes resultantes de las convoluciones, el resultado del filtro Sobel S es √(Ax² + Ay²).

Requisitos

Queremos que esta función tome cualquier matriz 2D y genere otra matriz 2D. Es posible que deseemos que funcione en cualquier par de ejes de un ndarray. Incluso es posible que deseemos que funcione en más (o menos) de dos ejes. Es posible que tengamos especificaciones sobre cómo manejar los bordes de la matriz.

• La Computación Cuántica podría recibir un impulso a partir del desc...
• Transistores de madera arraigan
• La Administración Biden considera imponer más restricciones en la v...

Por ahora, recordemos mantenerlo simple y comenzar con una implementación 2D. Pero antes de hacer eso, esbozaremos el diseño.

Esbozar el Diseño

Comenzamos con un diseño simple, aprovechando las anotaciones de Python. Recomiendo encarecidamente anotar tanto como sea posible y utilizar mypy como 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:    """Calcular el filtro Sobel de una imagen    Parámetros    ----------    arr : NDArray        Imagen de entrada    axes : Tuple[int, int], opcional        Ejes sobre los cuales calcular el filtro, por defecto (-2, -1)    Devoluciones    -------    NDArray        Salida    """    # Solo acepta matrices 2D    if arr.ndim != 2:        raise NotImplementedError    # Asegura que el eje[0] sea el primer eje y el eje[1] sea el segundo    # eje. El oscuro `normalize_axis_index` convierte índices negativos en    # índices entre 0 y arr.ndim - 1.    if any(        normalize_axis_index(ax, arr.ndim) != i        for i, ax in zip(range(2), axes)    ):        raise NotImplementedError    pass

Implementar Pruebas

Esta función no hace mucho. Pero está documentada, anotada y sus limitaciones actuales están incorporadas en ella. Ahora que tenemos un diseño, pasamos inmediatamente a las pruebas.

Primero, notamos que las imágenes vacías (todos ceros) no tienen bordes. Por lo tanto, también deben devolver ceros. De hecho, cualquier imagen constante también debería devolver ceros. Vamos a incorporar eso en nuestras primeras pruebas. También veremos cómo podemos usar pruebas de monkey para probar muchas matrices.

# test_zero_constant.pyimport numpy as npimport pytest# Prueba de múltiples dtypes a la [email protected](    "dtype",    ["float16", "float32", "float64", "float128"],)def test_zero(dtype):    # Establecer semilla aleatoria    seed = int(np.random.rand() * (2**32 - 1))    np.random.seed(seed)    # Crear una matriz 2D de forma aleatoria y llenar con ceros    nx, ny = np.random.randint(3, 100, size=(2,))    arr = np.zeros((nx, ny), dtype=dtype)    # Aplicar la función sobel    arr_sob = sobel(arr)    # `assert_array_equal` debería fallar la mayoría de las veces.    # Solo funcionará cuando `arr_sob` sea idénticamente cero,    # lo cual generalmente no es el caso.    # ¡NO USAR!    # np.testing.assert_array_equal(    #     arr_sob, 0.0, err_msg=f"{seed=} {nx=}, {ny=}"    # )    # `assert_almost_equal` puede fallar cuando se usa con decimales altos.    # También se basa en la comprobación de float64, que puede fallar para    # tipos float128.    # ¡NO USAR!    # np.testing.assert_almost_equal(    #     arr_sob,    #     np.zeros_like(arr),    #     err_msg=f"{seed=} {nx=}, {ny=}",    #     decimal=4,    # )    # `assert_allclose` con tolerancia personalizada es mi método preferido    # El 10 es arbitrario y depende del problema. Si un método    # que sabes que es correcto no pasa, aumenta a 100, etc.    # Si la tolerancia necesaria para que las pruebas pasen es demasiado alta, asegúrate de que el método sea realmente correcto.    tol = 10 * np.finfo(arr.dtype).eps    err_msg = f"{seed=} {nx=}, {ny=} {tol=}"  # Registrar semillas y otra información    np.testing.assert_allclose(        arr_sob,        np.zeros_like(arr),        err_msg=err_msg,        atol=tol,  # rtol es inútil para 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 es inútil para desired=zeros    )

Este fragmento de código se puede ejecutar desde la línea de comandos con

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

Versión Alfa

Por supuesto, nuestras pruebas fallarán actualmente, pero son suficientes por ahora. Vamos a implementar una primera versión.

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    # Definir nuestros núcleos de filtro. Observa que heredan el tipo de entrada, por lo    # que un input de float32 nunca tiene que ser convertido a float64 para su cálculo.    # ¿Pero puedes ver dónde usar otro dtype para Gx y Gy podría tener sentido para algunos dtypes de entrada?    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,    )    # Crear la matriz de salida y llenar con ceros    s = np.zeros_like(arr)    for ix in range(1, arr.shape[0] - 1):        for iy in range(1, arr.shape[1] - 1):            # Multiplicación elemento a elemento seguido de suma, también conocido como convolución            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

Con esta nueva función, todas nuestras pruebas deberían pasar y deberíamos obtener una salida como esta:

$ pytest -qq -s -x -vv --durations=0 test_zero_constant.py........======================================== slowest durations =========================================0.09s llamada     t_049988eae7f94139a7067f142bf2852f.py::test_constant[float16]0.08s llamada     t_049988eae7f94139a7067f142bf2852f.py::test_zero[float64]0.06s llamada     t_049988eae7f94139a7067f142bf2852f.py::test_constant[float128]0.04s llamada     t_049988eae7f94139a7067f142bf2852f.py::test_zero[float128]0.04s llamada     t_049988eae7f94139a7067f142bf2852f.py::test_constant[float64]0.02s llamada     t_049988eae7f94139a7067f142bf2852f.py::test_constant[float32]0.01s llamada     t_049988eae7f94139a7067f142bf2852f.py::test_zero[float16](17 duraciones < 0.005s ocultas. Use -vv para mostrar estas duraciones.)8 pasadas en 0.35s

Hasta ahora tenemos:

1. Requisitos recopilados de nuestro problema.
2. Realizado un diseño inicial.
3. Implementado algunas pruebas.
4. Implementado la versión alfa, que pasa estas pruebas.

Estas pruebas proporcionan verificación para futuras versiones, así como una biblioteca de oráculos muy rudimentaria. Pero si bien las pruebas unitarias son excelentes para capturar desviaciones mínimas de los resultados esperados, no son excelentes para diferenciar resultados incorrectos de resultados cuantitativamente diferentes, pero aún correctos. Supongamos que mañana queremos implementar otro operador de tipo Sobel, que tiene un kernel más largo. No esperamos que los resultados coincidan exactamente con nuestra versión actual, pero sí esperamos que las salidas de ambas funciones sean cualitativamente muy similares.

Además, es una excelente práctica probar muchas entradas diferentes en nuestras funciones para tener una idea de cómo se comportan para diferentes entradas. Esto garantiza que validemos científicamente los resultados.

Scikit-image tiene una excelente biblioteca de imágenes que podemos usar para crear oráculos.

# !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)    # Los datos se obtienen mediante llamadas a funciones    if isinstance(attr, Callable):        try:            # Descargar los datos            data = attr()            # Asegurarse de que es una matriz 2D            if isinstance(data, np.ndarray) and data.ndim == 2:                # Convertir de varios tipos de enteros a float32 para evaluar mejor la precisión                bwimages[attrname] = data.astype(np.float32)        except:            continue# Aplicar sobel a las imágenesbwimages_sobel = {k: sobel(v) for k, v in bwimages.items()}

Una vez que tenemos los datos, podemos representarlos gráficamente.

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")

¡Estos se ven muy bien! Recomendaría almacenarlos, tanto como matrices como figuras, para poder verificarlos rápidamente en caso de una nueva versión. Recomiendo encarecidamente usar HD5F para el almacenamiento de matrices. También puedes configurar una Sphinx Gallery para generar directamente las figuras al actualizar la documentación, de esta manera tendrás una construcción de figuras reproducible que puedes utilizar para comparar con versiones anteriores.

Después de validar los resultados, puedes almacenarlos en disco y utilizarlos como parte de tus pruebas unitarias. Algo así:

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,    )

Perfiles

¡Calcular el filtro Sobel para estos conjuntos de datos llevó tiempo! Por lo tanto, el siguiente paso es perfilar el código. Recomiendo crear un archivo benchmark_xyz.py para cada prueba y volver a ejecutarlos regularmente para detectar posibles regresiones de rendimiento. Esto incluso puede formar parte de tus pruebas unitarias, pero no iremos tan lejos en este ejemplo. Otra idea es utilizar los oráculos anteriores para realizar pruebas de rendimiento.

Existen muchas formas de medir el tiempo de ejecución del código. Para obtener el tiempo transcurrido real a nivel del sistema, utiliza perf_counter_ns del módulo integrado time para medir el tiempo en nanosegundos. En un cuaderno Jupyter, la magia de celda %%timeit integrada cronometra la ejecución de una determinada celda. El decorador a continuación está inspirado en esta magia de celda y se puede utilizar para medir el tiempo de cualquier función.

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:    """Aplica a una función para cronometrar sus ejecuciones.    Parámetros    ----------    func_or_number : Optional[Callable], opcional        Función a decorar o argumento `number` (ver más abajo), por        defecto None    number : int, opcional        Número de veces que se ejecutará la función para obtener estadísticas, por        defecto 10    Retorna    -------    Callable        Cuando se le proporciona una función, devuelve la función decorada.        De lo contrario, devuelve un decorador.    Ejemplos    --------    .. 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,))            # Ejecutar primero y medir almacenar el resultado            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)  # Sin almacenamiento, más rápido                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} por ejecución (media ± desviación estándar de {number} ejecuciones)"            )            return result        return wrapper    if func is not None:        return decorator(func)    return decorator

Poniendo a prueba nuestra función:

arr_test = np.random.randn(500, 500)sobel_timed = timeit(sobel)sobel_timed(arr_test);# 3.9s ± 848.9 ms por ejecución (media ± desv. estándar de 10 ejecuciones)

No es muy rápido 🙁

Cuando se investiga la lentitud o las regresiones de rendimiento, también es posible hacer un seguimiento del tiempo de ejecución de cada línea. La biblioteca line_profiler es un excelente recurso para esto. Puede ser utilizada mediante la magia de celdas Jupyter, o utilizando la API. Aquí hay un ejemplo de API:

from line_profiler import LineProfilerlp = LineProfiler()lp_wrapper = lp(sobel)lp_wrapper(arr_test)lp.print_stats(output_unit=1)  # 1 para segundos, 1e-3 para milisegundos, etc.

Aquí hay una salida de ejemplo:

Unidad de tiempo: 1 sTiempo total: 4.27197 sArchivo: /tmp/ipykernel_521529/1313985340.pyFunción: sobel en la línea 8Línea #      Veces         Tiempo  Por vez   % Tiempo  Contenido de la línea==============================================================     8                                           def sobel(arr: NDArray, axes: Tuple[int, int] = (-2, -1)) -> NDArray:     9                                               # Solo acepta matrices 2D    10         1          0.0      0.0      0.0      if arr.ndim != 2:    11                                                   raise NotImplementedError    12                                               13                                               # Asegura que el eje[0] sea el primer eje y que el eje[1] sea el segundo    14                                               # eje. El oscuro `normalize_axis_index` convierte los índices negativos en    15                                               # índices entre 0 y 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

Se pasa mucho tiempo dentro del bucle más interno. NumPy prefiere las operaciones matemáticas vectorizadas, ya que puede confiar en código compilado. Dado que estamos utilizando bucles explícitos, tiene sentido que esta función sea muy lenta.

Además, es importante ser inteligente en cuanto a las asignaciones de memoria dentro de los bucles. NumPy es algo inteligente al asignar pequeñas cantidades de memoria dentro de los bucles, por lo que las líneas que definen s1 o s2 podrían no estar asignando nueva memoria. Pero también podrían estarlo haciendo, o podría haber alguna otra asignación de memoria que esté ocurriendo internamente de la que no estamos conscientes. Por lo tanto, también recomiendo perfilar la memoria. Me gusta usar Memray para eso, pero hay otros como Fil y Sciagraph. También evitaría memory_profiler que (¡muy desafortunadamente!) ya no se mantiene. Además, Memray es mucho más potente. Aquí tienes un ejemplo de Memray en acción:

$ # Crear sobel.bin que contiene la información de perfilado$ memray run -fo sobel.bin --trace-python-allocators sobel_run.pyEscribiendo los resultados del perfil en sobel.binMemray ADVERTENCIA: Corrigiendo el símbolo para aligned_alloc de 0x7fc5c984d8f0 a 0x7fc5ca4a5ce0[memray] Se generaron exitosamente los resultados del perfil.Ahora puedes generar informes a partir de los registros de asignación almacenados.Algunos ejemplos de comandos para generar informes:python3 -m memray flamegraph sobel.bin

$ # Generar gráfico de llamas$ memray flamegraph -fo sobel_flamegraph.html --temporary-allocations sobel.bin⠙ Calculando el nivel máximo de uso de memoria... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸  99% 0:00:0100:01⠏ Procesando registros de asignación... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸  98% 0:00:0100:01Se escribió sobel_flamegraph.html

$ # Mostrar árbol de memoria$ memray tree --temporary-allocations sobel.bin⠧ Calculando el nivel máximo de uso de memoria... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸ 100% 0:00:0100:01⠧ Procesando registros de asignación... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸ 100% 0:00:0100:01Metadatos de asignación-------------------Argumentos de línea de comandos: 'memray run -fo sobel.bin --trace-python-allocators sobel_run.py'Tamaño máximo de memoria: 11.719MBNúmero de asignaciones: 15332714Las 10 asignaciones más grandes:-----------------------📂 123.755MB (100.00 %) <ROOT>  └── [[3 frames ocultos en 2 archivo(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 ocultos en 1 archivo(s)]]        │   │   └── 📄 18.922MB (15.29 %) _sum          │   │       lib/python3.10/site-packages/numpy/core/_methods.py:49        │   └── [[1 frames ocultos en 1 archivo(s)]]        │       └── 📄 18.921MB (15.29 %) _sum          │           lib/python3.10/site-packages/numpy/core/_methods.py:49...

Versión Beta y un Error

Ahora que tenemos una versión alfa funcional y algunas funciones de perfilado, aprovecharemos la biblioteca SciPy para obtener una versión mucho más rápida.

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        # Crear los núcleos como un único arreglo complejo. Nos permite usar    # np.abs en lugar de np.hypot para calcular la magnitud.    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)  # Absoluto en el lugar    return s.real

sobel_timed = timeit(sobel_conv2d)sobel_timed(arr_test)# 14.3 ms ± 1.71 ms por iteración (media ± desviación estándar de 10 ejecuciones)

¡Mucho mejor! Pero ¿es correcto?

Las imágenes se ven muy similares, pero si te fijas en la escala de colores, no son iguales. Al ejecutar las pruebas se detecta un pequeño error promedio. Afortunadamente, ahora estamos bien equipados para detectar diferencias cuantitativas y cualitativas.

Después de investigar este error, lo atribuimos a las diferentes condiciones de límites. Al mirar la documentación de convolve2d nos dice que la matriz de entrada se rellena con ceros antes de aplicar el kernel. En la versión alfa, rellenamos la salida.

¿Cuál es la correcta? Podría argumentarse que la implementación de SciPy tiene más sentido. En este caso, deberíamos adoptar la nueva versión como el oráculo, y si es necesario, corregir la versión alfa para que coincida con ella. Esto es común en el desarrollo de software científico: nueva información sobre cómo hacer las cosas mejor cambia los oráculos y las pruebas.

En este caso, la solución es sencilla, rellenar la matriz con ceros antes de procesarla.

def sobel_v2(arr: NDArray, axes: Tuple[int, int] = (-2, -1)) -> NDArray:    # ...    arr = np.pad(arr, (1,))  # Después de rellenar, tiene forma (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)  # Ajustar índices    return s

Una vez corregida nuestra función, podemos actualizar los oráculos y las pruebas que dependen de ellos.

Pensamientos Finales

Vimos cómo poner en práctica algunas de las ideas de desarrollo de software exploradas en este artículo. También vimos algunas herramientas que puedes usar para desarrollar código de alta calidad y alto rendimiento.

Te sugiero que pruebes algunas de estas ideas en tus propios proyectos. En particular, practica el perfilado de código y mejóralo. La función del filtro Sobel que codificamos es muy ineficiente, te sugiero que intentes mejorarla.

Por ejemplo, prueba la paralelización de CPU con un compilador JIT como Numba, traslada el bucle interno a Cython o implementa una función de GPU CUDA con Numba o CuPy. Siéntete libre de visitar mi serie sobre cómo programar núcleos CUDA con Numba.

---

• Hands On Tutorials
• Python
• Scientific Computing
• software development
• Test Driven Development