Parte 2 - Arrays de NumPy

¿Qué es NumPy?

  • Una biblioteca para Python: ndarray + ufunc
  • Los arrays multidimensionales (ndarray) nos permiten almacenar datos de manera estructurada
  • Las funciones universales (ufunc) nos permiten operar con esos datos de manera eficiente

Python está organizado en módulos, que son archivos con extensión .py que contienen funciones, variables y otros objetos, y paquetes, que son conjuntos de módulos. Cuando queremos utilizar objetos que están definidos en un módulo tenemos que importarlo, y una vez que lo hemos hecho podemos usar el operador . para ir descendiendo en la jerarquía de paquetes y acceder al objeto que necesitamos. Por ejemplo, de esta manera importamos NumPy:

In [1]:
import numpy

Y de esta manera accedemos a la función norm, que calcula la norma (o módulo) de un array:

In [2]:
numpy.linalg.norm
Out[2]:
<function numpy.linalg.linalg.norm>

La función norm está dentro del paquete linalg, que a su vez está dentro del paquete NumPy.

La convención para importar NumPy siempre es esta:

In [3]:
import numpy as np

Lo que hacemos es crear un alias al paquete NumPy de nombre np. Es simplemente una forma de abreviar el código. Esta forma de separar las funciones en paquetes (que se llaman espacios de nombres o namespaces) conduce a una mayor legibilidad del código y a la supresión de ambigüedades.

Para encontrar ayuda sobre cierto tema podemos usar la función lookfor:

In [4]:
np.lookfor("solve")
Search results for 'solve'
--------------------------
numpy.linalg.solve
    Solve a linear matrix equation, or system of linear scalar equations.
numpy.linalg.lstsq
    Return the least-squares solution to a linear matrix equation.
numpy.linalg.tensorsolve
    Solve the tensor equation ``a x = b`` for x.
numpy.linalg._umath_linalg.solve
    solve the system a x = b, on the last two dimensions, broadcast to the rest.
numpy.linalg._umath_linalg.solve1
    solve the system a x = b, for b being a vector, broadcast in the outer dimensions.
numpy.distutils.misc_util.njoin
    Join two or more pathname components +
numpy.distutils.misc_util.minrelpath
    Resolve `..` and '.' from path.
numpy.distutils.system_info.UmfpackNotFoundError
    UMFPACK sparse solver (http://www.cise.ufl.edu/research/sparse/umfpack/)
numpy.linalg.pinv
    Compute the (Moore-Penrose) pseudo-inverse of a matrix.
numpy.linalg.cholesky
    Cholesky decomposition.
numpy.linalg.tensorinv
    Compute the 'inverse' of an N-dimensional array.
numpy.linalg.LinAlgError
    Generic Python-exception-derived object raised by linalg functions.
numpy.polynomial.hermite.hermfit
    Least squares fit of Hermite series to data.
numpy.polynomial.laguerre.lagfit
    Least squares fit of Laguerre series to data.
numpy.polynomial.legendre.legfit
    Least squares fit of Legendre series to data.
numpy.polynomial.chebyshev.chebfit
    Least squares fit of Chebyshev series to data.
numpy.polynomial.hermite_e.hermefit
    Least squares fit of Hermite series to data.
numpy.polynomial.polynomial.polyfit
    Least-squares fit of a polynomial to data.

Constantes y funciones matemáticas

Además de arrays, NumPy contiene también constantes y funciones matemáticas de uso cotidiano.

In [5]:
np.e
Out[5]:
2.718281828459045
In [6]:
np.pi
Out[6]:
3.141592653589793
In [7]:
np.log(2)
Out[7]:
0.69314718055994529

¿Qué es exactamente un array?

Un array de NumPy es una colección de N elementos, igual que una secuencia de Python (por ejemplo, una lista). Tiene las mismas propiedades que una secuencia y alguna más. Para crear un array, la forma más directa es pasarle una secuencia a la función np.array.

In [8]:
np.array([1, 2, 3])
Out[8]:
array([1, 2, 3])

Los arrays de NumPy son homogéneos, es decir, todos sus elementos son del mismo tipo. Si le pasamos a np.array una secuencia con objetos de tipos diferentes, promocionará todos al tipo con más información. Para acceder al tipo del array, podemos usar el atributo dtype.

In [9]:
a = np.array([1, 2, 3.0])
print(a.dtype)
float64
In [10]:
np.array([1, 2, "3"])
Out[10]:
array(['1', '2', '3'], 
      dtype='<U1')
**Nota**: Si NumPy no entiende el tipo de datos o construimos un array con argumentos incorrectos devolverá un array con `dtype` `object`. Estos arrays rara vez son útiles y su aparición suele ser signo de que algo falla en nuestro programa.

NumPy intentará automáticamente construir un array con el tipo adecuado teniendo en cuenta los datos de entrada, aunque nosotros podemos forzarlo.

In [11]:
np.array([1, 2, 3], dtype=float)
Out[11]:
array([ 1.,  2.,  3.])
In [12]:
np.array([1, 2, 3], dtype=complex)
Out[12]:
array([ 1.+0.j,  2.+0.j,  3.+0.j])

También podemos convertir un array de un tipo a otro utilizando el método .astype.

In [13]:
a
Out[13]:
array([ 1.,  2.,  3.])
In [14]:
a.astype(int)
Out[14]:
array([1, 2, 3])

Motivo: eficiencia

  • Los bucles son costosos
  • Eliminar bucles: vectorizar operaciones
  • Los bucles se ejecutan en Python, las operaciones vectorizadas en C
  • Las operaciones entre arrays de NumPy se realizan elemento a elemento

Ejemplo:

$$ a_{ij} = b_{ij} + c_{ij} $$
In [15]:
N, M = 100, 100
a = np.empty(10000).reshape(N, M)
b = np.random.rand(10000).reshape(N, M)
c = np.random.rand(10000).reshape(N, M)
In [16]:
%%timeit
for i in range(N):
    for j in range(M):
        a[i, j] = b[i, j] + c[i, j]
100 loops, best of 3: 8.75 ms per loop
In [17]:
%%timeit
a = b + c
100000 loops, best of 3: 8.34 µs per loop

¡1000 veces más rápido! Se hace fundamental vectorizar las operaciones y aprovechar al máximo la velocidad de NumPy.

Indexación de arrays

Una de las herramientas más importantes a la hora de trabajar con arrays es el indexado. Consiste en seleccionar elementos aislados o secciones de un array. Nosotros vamos a ver la indexación básica, pero existen técnicas de indexación avanzada que convierten los arrays en herramientas potentísimas.

In [18]:
a = np.array([
    [1, 2, 3],
    [4, 5, 6]
])
a
Out[18]:
array([[1, 2, 3],
       [4, 5, 6]])

Los índices se indican entre corchetes justo después del array. Recuerda que en Python la indexación empieza en 0. Si recuperamos el primer elemento de un array de dos dimensiones, obtenemos la primera fila.

In [19]:
a[0]
Out[19]:
array([1, 2, 3])

En vez de usar a[0][0] para recuperar el primer elemento de la primera fila, podemos abreviar aún más la sintaxis:

In [20]:
a[0, 0]
Out[20]:
1

No solo podemos recuperar un elemento aislado, sino también porciones del array, utilizando la sintaxis [<inicio>:<final>:<salto>].

In [21]:
a[0, 1:3]
Out[21]:
array([2, 3])
In [22]:
a[0, ::2]
Out[22]:
array([1, 3])

Creación de arrays

Muchos métodos y muy variados

  • A partir de datos existentes: array, copy
  • Unos y ceros: empty, eye, ones, zeros, *_like
  • Rangos: arange, linspace, logspace, meshgrid
  • Aleatorios: rand, randn

Unos y ceros

  • empty(shape) crea un array con «basura», equivalente a no inicializarlo, ligeramente más rápido que zeros o ones
  • eye(N, M=None, k=0) crea un array con unos en una diagonal y ceros en el resto
  • identity(n) devuelve la matriz identidad
  • Las funciones *_like(a) construyen arrays con el mismo tamaño que uno dado
In [23]:
np.identity(5).astype(int)
Out[23]:
array([[1, 0, 0, 0, 0],
       [0, 1, 0, 0, 0],
       [0, 0, 1, 0, 0],
       [0, 0, 0, 1, 0],
       [0, 0, 0, 0, 1]])
In [24]:
_.shape
Out[24]:
(5, 5)

Si la función recibe como argumento shape, debemos pasarle el número de filas y columnas como una tupla (es decir, encerrado entre paréntesis).

In [25]:
np.zeros((3, 4))
Out[25]:
array([[ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.]])
**Nota**: Un error muy típico es tratar de construir un array llamando a la función con dos argumentos, como se ejemplifica en la celda siguiente. Esto produce un error, porque NumPy espera un solo argumento: una tupla con el número de filas y el número de columnas. Es conveniente asegurarse de cuál es el convenio en cada caso porque no siempre hay consistencia interna.
In [26]:
np.zeros(3, 4)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-26-06beb765944a> in <module>()
----> 1 np.zeros(3, 4)

TypeError: data type not understood
In [27]:
np.ones((3, 4))
Out[27]:
array([[ 1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.]])
In [28]:
i3 = np.identity(3)
i3
Out[28]:
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])
In [29]:
i3.shape
Out[29]:
(3, 3)
In [30]:
np.ones(i3.shape)
Out[30]:
array([[ 1.,  1.,  1.],
       [ 1.,  1.,  1.],
       [ 1.,  1.,  1.]])

Si en lugar de pasar directamente la forma del array ya sabemos que queremos crear un array con la misma forma que otro, podemos usar las funciones *_like, que reciben un array en vez de una tupla.

In [31]:
np.ones_like(i3)
Out[31]:
array([[ 1.,  1.,  1.],
       [ 1.,  1.,  1.],
       [ 1.,  1.,  1.]])

Rangos

  • linspace(start, stop, num=50) devuelve números equiespaciados dentro de un intervalo
  • logspace(start, stop, num=50, base=10.0) devuelve números equiespaciados según una escala logarítmica
  • meshgrid(x1, x2, ...) devuelve matrices de n-coordenadas
In [32]:
np.linspace(0, 1, num=10)
Out[32]:
array([ 0.        ,  0.11111111,  0.22222222,  0.33333333,  0.44444444,
        0.55555556,  0.66666667,  0.77777778,  0.88888889,  1.        ])
In [33]:
np.logspace(0, 3)
Out[33]:
array([    1.        ,     1.1513954 ,     1.32571137,     1.52641797,
           1.75751062,     2.02358965,     2.32995181,     2.6826958 ,
           3.0888436 ,     3.55648031,     4.09491506,     4.71486636,
           5.42867544,     6.25055193,     7.19685673,     8.28642773,
           9.54095476,    10.98541142,    12.64855217,    14.56348478,
          16.76832937,    19.30697729,    22.22996483,    25.59547923,
          29.47051703,    33.93221772,    39.06939937,    44.98432669,
          51.79474679,    59.63623317,    68.6648845 ,    79.06043211,
          91.0298178 ,   104.81131342,   120.67926406,   138.94954944,
         159.98587196,   184.20699693,   212.09508879,   244.20530945,
         281.1768698 ,   323.74575428,   372.75937203,   429.19342601,
         494.17133613,   568.9866029 ,   655.12855686,   754.31200634,
         868.51137375,  1000.        ])

La función np.meshgrid se utiliza mucho a la hora de representar funciones en dos dimensiones, y crea dos arrays: uno varía por filas y otro por columnas. Combinándolos, podemos evaluar la función en un cuadrado.

In [37]:
x = np.linspace(0, 1, num=5)
y = np.linspace(0, 1, num=5)

xx, yy = np.meshgrid(x, y)
In [38]:
xx, yy
Out[38]:
(array([[ 0.  ,  0.25,  0.5 ,  0.75,  1.  ],
       [ 0.  ,  0.25,  0.5 ,  0.75,  1.  ],
       [ 0.  ,  0.25,  0.5 ,  0.75,  1.  ],
       [ 0.  ,  0.25,  0.5 ,  0.75,  1.  ],
       [ 0.  ,  0.25,  0.5 ,  0.75,  1.  ]]),
 array([[ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
       [ 0.25,  0.25,  0.25,  0.25,  0.25],
       [ 0.5 ,  0.5 ,  0.5 ,  0.5 ,  0.5 ],
       [ 0.75,  0.75,  0.75,  0.75,  0.75],
       [ 1.  ,  1.  ,  1.  ,  1.  ,  1.  ]]))
In [39]:
xx + 1j * yy
Out[39]:
array([[ 0.00+0.j  ,  0.25+0.j  ,  0.50+0.j  ,  0.75+0.j  ,  1.00+0.j  ],
       [ 0.00+0.25j,  0.25+0.25j,  0.50+0.25j,  0.75+0.25j,  1.00+0.25j],
       [ 0.00+0.5j ,  0.25+0.5j ,  0.50+0.5j ,  0.75+0.5j ,  1.00+0.5j ],
       [ 0.00+0.75j,  0.25+0.75j,  0.50+0.75j,  0.75+0.75j,  1.00+0.75j],
       [ 0.00+1.j  ,  0.25+1.j  ,  0.50+1.j  ,  0.75+1.j  ,  1.00+1.j  ]])

Operaciones con arrays

Las funciones universales (ufunc) operan sobre arrays de NumPy elemento a elemento y siguiendo las reglas de broadcasting.

  • Funciones matemáticas: sin, cos, sqrt, exp, ...
  • Operaciones lógicas: <, ~, ...
  • Funciones lógicas: all, any, isnan, allclose, ...
**Nota**: Las funciones matemáticas siempre devuelven el mismo tipo de datos de entrada
In [40]:
a = np.arange(2 * 3).reshape(2, 3)
a
Out[40]:
array([[0, 1, 2],
       [3, 4, 5]])
In [41]:
np.sqrt(a)
Out[41]:
array([[ 0.        ,  1.        ,  1.41421356],
       [ 1.73205081,  2.        ,  2.23606798]])
In [42]:
np.sqrt(np.arange(-3, 3))
-c:1: RuntimeWarning: invalid value encountered in sqrt
Out[42]:
array([        nan,         nan,         nan,  0.        ,  1.        ,
        1.41421356])
In [43]:
np.arange(-3, 3).astype(complex)
Out[43]:
array([-3.+0.j, -2.+0.j, -1.+0.j,  0.+0.j,  1.+0.j,  2.+0.j])
In [44]:
np.sqrt(_)
Out[44]:
array([ 0.00000000+1.73205081j,  0.00000000+1.41421356j,
        0.00000000+1.j        ,  0.00000000+0.j        ,
        1.00000000+0.j        ,  1.41421356+0.j        ])

Funciones de comparación

Las comparaciones devuelven un array de booleanos:

In [45]:
a = np.arange(6)
b = np.ones(6).astype(int)
a, b
Out[45]:
(array([0, 1, 2, 3, 4, 5]), array([1, 1, 1, 1, 1, 1]))
In [46]:
a < b
Out[46]:
array([ True, False, False, False, False, False], dtype=bool)
In [47]:
np.any(a < b)
Out[47]:
True
In [48]:
np.all(a < b)
Out[48]:
False
In [49]:
a = np.arange(6).astype(float)
b = np.ones(6)
a, b
Out[49]:
(array([ 0.,  1.,  2.,  3.,  4.,  5.]), array([ 1.,  1.,  1.,  1.,  1.,  1.]))

Las funciones isclose y allclose realizan comparaciones entre arrays especificando una tolerancia:

In [50]:
np.isclose(a, b, rtol=1e-6)
Out[50]:
array([False,  True, False, False, False, False], dtype=bool)
In [51]:
np.allclose(a, b, rtol=1e-6)
Out[51]:
False
**¡Importante!** Ni en Python ni en ningún otro lenguaje debemos hacer comparaciones exactas entre números de punto flotante. Las operaciones matemáticas con estos números producen casi siempre resultados poco intuitivos y hay que tener cuidado con ellas. Para una introducción a estas peculiaridades existe la web http://puntoflotante.org/.
In [52]:
0.1 + 0.2 + 0.3
Out[52]:
0.6000000000000001
In [53]:
0.3 + 0.2 + 0.1
Out[53]:
0.6
In [54]:
0.1 + 0.2 + 0.3 == 0.3 + 0.2 + 0.1
Out[54]:
False

Ejercicios

Ejercicio 1.

  1. Crear un array z1 3x4 lleno de ceros de tipo entero.
  2. Crear un array z2 3x4 lleno de ceros salvo la primera fila que serán todo unos.
  3. Crear un array z3 3x4 lleno de ceros salvo la última fila que será el rango entre 5 y 8.
In [55]:
a = np.zeros((3, 4))
a
Out[55]:
array([[ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.]])
In [56]:
a[0, :] = 1
a
Out[56]:
array([[ 1.,  1.,  1.,  1.],
       [ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.]])
In [57]:
b = np.zeros((3, 4))
b[-1] = np.arange(5, 9)
b
Out[57]:
array([[ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.],
       [ 5.,  6.,  7.,  8.]])

Ejercicio 2.

  1. Crea un vector de 10 elementos, siendo los impares unos y los pares doses.
  2. Crea un «tablero de ajedrez», con unos en las casillas negras y ceros en las blancas.
In [58]:
v = np.ones(10)
v
Out[58]:
array([ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.])
In [59]:
v[::2] = 2
v
Out[59]:
array([ 2.,  1.,  2.,  1.,  2.,  1.,  2.,  1.,  2.,  1.])
In [60]:
tablero = np.zeros((8, 8))
tablero
Out[60]:
array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])
In [61]:
tablero[1::2, ::2] = 1
tablero[::2, 1::2] = 1
tablero
Out[61]:
array([[ 0.,  1.,  0.,  1.,  0.,  1.,  0.,  1.],
       [ 1.,  0.,  1.,  0.,  1.,  0.,  1.,  0.],
       [ 0.,  1.,  0.,  1.,  0.,  1.,  0.,  1.],
       [ 1.,  0.,  1.,  0.,  1.,  0.,  1.,  0.],
       [ 0.,  1.,  0.,  1.,  0.,  1.,  0.,  1.],
       [ 1.,  0.,  1.,  0.,  1.,  0.,  1.,  0.],
       [ 0.,  1.,  0.,  1.,  0.,  1.,  0.,  1.],
       [ 1.,  0.,  1.,  0.,  1.,  0.,  1.,  0.]])

Ejercicio 3.

  1. Crea una matriz aleatoria 5x5 y halla los valores mínimo y máximo.
  2. Normaliza esa matriz entre 0 y 1.

Bibliografía