Numpy contiguo

En una charla que dí sobre numpy comenté que los datos en un numpy array son homogéneos, es decir, del mismo tipo, y contiguos en memoria.

Una persona entre los asistentes me comentó que no creía que fueran contiguos en memoria al 100%. Que creía que la situación podía ser de la siguiente forma:

Yo le comenté que pensaba que eso daría un MemoryError ya que no había espacio suficiente para poder meter todo el array de forma contigua en RAM. Pero uno nunca puede estar seguro de todo y por eso le dije que ‘creía’.

Todo esto que viene a continuación es un intento de explicación de porqué los numpy arrays parece que son 100% contiguos en memoria a día de hoy.

Evidencias

Primero de todo voy a crear un numpy array de una dimensión:

import numpy as np

arr = np.arange(10)

print(arr)

import numpy as np

arr = np.arange(10)

print(arr)

Lo anterior mostrará:

[0 1 2 3 4 5 6 7 8 9]

1	[0 1 2 3 4 5 6 7 8 9]

El objeto arr dispone de mucha información, además de los datos en sí. Si inspecciono algunos de los metadatos de un numpy array puedo obtener la siguiente información interesante:

print(f'posición del inicio del array en memoria: {arr.data}')

1	print(f'posición del inicio del array en memoria: {arr.data}')

Lo anterior debería mostrar algo como:

posición del inicio del array en memoria: <memory at 0x0000023E3E195D08>

1	posición del inicio del array en memoria: <memory at 0x0000023E3E195D08>

Lo anterior nos informa de una única dirección en memoria y siempre es así (hasta hoy, en el futuro ya veremos, luego más sobre esto). Si el array no fuera contiguo necesitaríamos más información para poder acceder a todos los elementos.

Puedo intentar ver si los strides me dan algo de información sobre cómo tendría que hacer un salto intermedio si el array no fuera contiguo. Veamos la información que nos ofrece numpy sobre los strides:

print(f'Strides del array: {arr.strides}')

1	print(f'Strides del array: {arr.strides}')

Lo anterior mostrará:

Strides del array: (4,)

1	Strides del array: (4,)

Si miramos la documentación sobre los strides nos dice:

Tuple of bytes to step in each dimension when traversing an array.

Que traducido libremente vendría a decir:

Tupla indicando el número de bytes que debemos ‘caminar’ en cada dimensión para recorrer el array en cada una de sus dimensiones.

Nuevamente, no vemos ningún indicio de que los strides contengan información en el hipotético caso de un array que no sea 100% contiguo en memoria.

Pero en el PEP-3118 (buffer-protocol) hablan sobre suboffsets que podrían ser usados potencialmente por numpy o PIL/Pillow para acceder a ‘cosas’ no contiguas en memoria dando esa información de los suboffsets.

El protocolo buffer define cómo diferentes objetos pueden acceder a información en memoria y compartir esta información y es una de las bases de que numpy pueda ser muy eficiente al saltarse muchos pasos y poder operar a un nivel más bajo (C) sin necesidad de crear copias de datos.

¿Qué es esto de los `suboffsets`?

print(arr.data.suboffsets)

1	print(arr.data.suboffsets)

Lo anterior nos dará:

()

()

En el caso de arr me devuelve algo nulo. Si buscamos suboffsets en el código de numpy parece que siempre lo inicializa a lo mismo, NULL, y en los tests siempre hacen algo similar a assert_equal(arr.suboffsets, ()). Viendo eso parece que en numpy no se ha prestado mucha atención a la posibilidad de añadir soporte para suboffsets cuando se usa el protocolo buffer para crear el objeto ndarray.

Aquí podéis ver un PR todavía abierto (en el momento de escribir este artículo) donde se habla de añadir soporte a esto. También hay un issue abierto (en el momento de escribir este artículo) hablando sobre esto mismo. También, en la documentación sobre el protocolo buffer de CPython hablan sobre arrays complejos citando explícitamente a numpy y a PIL/Pillow pero parece que el tema de los suboffsets se restringe a PIL/Pillow.

Parece que, potencialmente, podríamos tener numpy arrays no contiguos en memoria pero parece que a día de hoy (en el momento de escribir este artículo) no es así. El PEP-3118 es de 2006 y parece que, después de 14 años desde el PEP, no ha habido nadie que haya encontrado una necesidad muy clara para implementar esto.

Como dijo Jake VanderPlas en su blog en este excelente artículo:

it is fundamentally the Buffer Protocol and related NumPy functionality that make Python useful as a scientific computing platform

En castellano:

Básicamente, son el protocolo buffer y la funcionalidad asociada que proporciona Numpy lo que hace que Python sea útil como plataforma de cálculo científico.

Demuéstrame que están contiguos en memoria

El MemoryView que me devuelve arr.data dispone del método tobytes que me permite convertir los datos que tiene el numpy array en un churro de bytes:

churro_bytes = arr.data.tobytes()

print(churro_bytes)

churro_bytes = arr.data.tobytes()

print(churro_bytes)

Lo anterior dara:

b'\x00\x00\x00\x00\x01\x00\x00\x00\x02\x00\x00\x00\x03\x00\x00\x00\x04\x00\x00\x00\x05\x00\x00\x00\x06\x00\x00\x00\x07\x00\x00\x00\x08\x00\x00\x00\t\x00\x00\x00'

1	b'\x00\x00\x00\x00\x01\x00\x00\x00\x02\x00\x00\x00\x03\x00\x00\x00\x04\x00\x00\x00\x05\x00\x00\x00\x06\x00\x00\x00\x07\x00\x00\x00\x08\x00\x00\x00\t\x00\x00\x00'

Necesitaré también saber como están estructurados los bytes para poder reconstruir los datos. Para ello voy a usar el objeto dtype accesible desde arr:

print(arr.dtype.byteorder)
print(arr.dtype.kind)
print(len(arr))

print(arr.dtype.byteorder)

print(arr.dtype.kind)

print(len(arr))

Que mostrará en pantalla:

=
i
10

Lo primero está relacionado con esto: https://docs.python.org/3/library/struct.html#byte-order-size-and-alignment

Lo segundo está relacionado con esto: https://docs.python.org/3/library/struct.html#format-characters

Lo tercero me dice cuantas veces se repetirá la estructura para reconstruir cada uno de los datos.

Si lo junto todo de la siguiente forma:

formato_simple = f'{arr.dtype.byteorder}{arr.dtype.kind}'

print(formato_simple)

formato_simple = f'{arr.dtype.byteorder}{arr.dtype.kind}'

print(formato_simple)

obtenemos:

=i

Voy a usar la biblioteca struct para convertir esos bytes nuevamente en la información consumible para mortales. Para ello necesitaré saber interpretar qué me estan queriendo decir esos bytes:

import struct

for pack in struct.iter_unpack(formato_simple, churro_bytes):
    print(pack)

import struct

for pack in struct.iter_unpack(formato_simple, churro_bytes):

print(pack)

Lo anterior mostrará:

(0,)
(1,)
(2,)
(3,)
(4,)
(5,)
(6,)
(7,)
(8,)
(9,)

(0,)

(1,)

(2,)

(3,)

(4,)

(5,)

(6,)

(7,)

(8,)

(9,)

Usando la función iter_unpack puedo ir reconstruyendo uno a uno los elementos. Si lo quiero hacer todo de una vez puedo usar directamente la función unpack pero en el formato necesito indicar la longitud total por lo que modifico el valor del formato para indicar que tenemos 10 elementos:

formato = f'{arr.dtype.byteorder}{len(arr)}{arr.dtype.kind}'

print(formato)

print(struct.unpack(formato, churro_bytes))

formato = f'{arr.dtype.byteorder}{len(arr)}{arr.dtype.kind}'

print(formato)

print(struct.unpack(formato, churro_bytes))

Lo anterior mostrará:

=10i
(0, 1, 2, 3, 4, 5, 6, 7, 8, 9)

1 2	=10i (0, 1, 2, 3, 4, 5, 6, 7, 8, 9)

Pero lo anterior no me asegura para nada que eso sea contiguo en memoria. El método tobytes, podría haber hecho algo de magia para devolverme datos no contiguos representados como un único churro. Es cierto, aunque es fácil de ver que no es así pero he encontrado una excusa para hablar algo de los dtypes 😛

Veamos ahora a ver si os convenzo. En este caso voy a necesitar la biblioteca ctypes. En concreto voy a usar la función string_at:

from ctypes import string_at

1	from ctypes import string_at

La documentación oficial dice:

This function returns the C string starting at memory address address as a bytes object. If size is specified, it is used as size, otherwise the string is assumed to be zero-terminated.

que traducido sería algo como:

Esta función devuelve la cadena C que comienza en la dirección de memoria address como un objeto bytes. Si se especifica size se usará como tamaño si no se asume que la cadena termina en cero.

Parece que si obtengo la posición del dato en memoria y su tamaño en bytes podría convertir esa información. Pues voy a recorrer desde el punto en que sé donde comienzan los datos y sabiendo el tamaño que debe ocupar cada elemento veremos si obtenemos lo mismo que en los datos originales.

La posición inicial de la información la podemos obtener de la interfaz del array:

_donde = arr.__array_interface__['data']

print(_donde)
`
# Me quedo solo con la primera parte
# que es la que me interesa ahora
donde = _donde[0]

_donde = arr.__array_interface__['data']

print(_donde)

# Me quedo solo con la primera parte

# que es la que me interesa ahora

donde = _donde[0]

El print mostrará:

(2466278768880, False)

1	(2466278768880, False)

Lo que nos interesa aquí es el primer número, que indica la dirección del primer dato en memoria.

Por otro lado, ¿cuántos bytes usa cada elemento?

size = arr.itemsize
print(size)

1 2	size = arr.itemsize print(size)

Serán:

4

Cada elemento hace uso de 4 bytes. Vamos a usar la función atring_at para ver si puedo reproducir el valor del primer elemento de arr (recordad que el primer valor es 0):

print(string_at(donde, size))

1	print(string_at(donde, size))

Lo anterior muestra:

b'\x00\x00\x00\x00'

1	b'\x00\x00\x00\x00'

Lo anterior me daría la representación en bytes. Puedo volver a usar la función struct.unpack y el formato_simple (ahora estoy desempaquetando solo un elemento, no los diez):

print(struct.unpack(formato_simple, string_at(donde, size)))

1	print(struct.unpack(formato_simple, string_at(donde, size)))

Eso dará:

(0,)

(0,)

¡¡¡Ta-daaaaaa!!! Parece que hemos reconstruido el valor solo conociendo su posición en memoria y su tamaño en bytes. ¿Los demás son contiguos?

print(
    'iteración, '
    'posición en memoria, '
    'valor del array, '
    'valor reconstruido'
)
for i in range(len(arr)):
    _donde = donde + i * size
    print(
        i, 
        _donde, 
        arr[i],
        struct.unpack(formato_simple, string_at(_donde, size))
    )

print(

'iteración, '

'posición en memoria, '

'valor del array, '

'valor reconstruido'

)

for i in range(len(arr)):

_donde = donde + i * size

print(

_donde,

arr[i],

struct.unpack(formato_simple, string_at(_donde, size))

)

Lo anterior dará:

iteración, posición en memoria, valor del array, valor reconstruido
0 2466278768880 0 (0,)
1 2466278768884 1 (1,)
2 2466278768888 2 (2,)
3 2466278768892 3 (3,)
4 2466278768896 4 (4,)
5 2466278768900 5 (5,)
6 2466278768904 6 (6,)
7 2466278768908 7 (7,)
8 2466278768912 8 (8,)
9 2466278768916 9 (9,)

iteración, posición en memoria, valor del array, valor reconstruido

0 2466278768880 0 (0,)

1 2466278768884 1 (1,)

2 2466278768888 2 (2,)

3 2466278768892 3 (3,)

4 2466278768896 4 (4,)

5 2466278768900 5 (5,)

6 2466278768904 6 (6,)

7 2466278768908 7 (7,)

8 2466278768912 8 (8,)

9 2466278768916 9 (9,)

Conclusiones

Parece que, potencialmente, los numpy arrays podrían ser discontínuos pero no es así. El buffer protocol es una de las cosas que ha hecho que (C)Python sea una herramienta importante en el mundo científico por su facilidad de interactuar con código muy eficiente escrito en lenguajes de programación de más bajo nivel y de evitar la necesidad de hacer copias de la información y poder compartir memoria.

No sé si esto contestará a la persona que hizo aquella objeción pero yo he aprendido unas cuantas cosas escribiendo este artículo 🙂

Evidencias

¿Qué es esto de los suboffsets?

Demuéstrame que están contiguos en memoria

Conclusiones

Deja una respuesta Cancelar la respuesta

¿Qué es esto de los `suboffsets`?