El paquete (módulo) numpy es usado en casi todos los cálculos numéricos usando Python. Es un paquete que provee a Python de estructuras de datos vectoriales, matriciales, y de rango mayor, de alto rendimiento. Está implementado en C y Fortran, de modo que cuando los cálculos son vectorizados (formulados con vectores y matrices), el rendimiento es muy bueno.
Para usar numpy necesitamos importar el módulo usando, por ejemplo:
In [276]:
import numpy as np
En el paquete numpy la terminología usada para vectores, matrices y conjuntos de datos de dimensión mayor es la de un arreglo.
numpyExisten varias formas para inicializar nuevos arreglos de numpy, por ejemplo desde
arange, linspace, etc.Por ejemplo, para crear nuevos arreglos de matrices y vectores desde listas Python podemos usar la función numpy.array.
In [277]:
# un vector: el argumento de la función array es una lista de Python
v = np.array([1,2,3,4])
v
Out[277]:
In [278]:
# una matriz: el argumento de la función array es una lista anidada de Python
M = np.array([[1, 2], [3, 4]])
M
Out[278]:
Los objetos v y M son ambos del tipo ndarray que provee el módulo numpy.
In [279]:
type(v), type(M)
Out[279]:
La diferencia entre los arreglos v y M es sólo su forma. Podemos obtener información de la forma de un arreglo usando la propiedad ndarray.shape
In [280]:
v.shape
Out[280]:
In [281]:
M.shape
Out[281]:
El número de elementos de un arreglo puede obtenerse usando la propiedad ndarray.size:
In [282]:
M.size
Out[282]:
Equivalentemente, podemos usar las funciones numpy.shape y numpy.size
In [283]:
np.shape(M)
Out[283]:
In [284]:
np.size(M)
Out[284]:
Hasta el momento el arreglo numpy.ndarray luce como una lista Python (anidada). Entonces, ¿por qué simplemente no usar listas para hacer cálculos en lugar de crear un tipo nuevo de arreglo?
Existen varias razones:
numpy usando lenguajes compilados (se usan C y Fortran).Usando la propiedad dtype (tipo de dato) de un ndarray, podemos ver qué tipo de dato contiene un arreglo:
In [285]:
M.dtype
Out[285]:
Se obtiene un error si intentamos asignar un valor de un tipo equivocado a un elemento de un arreglo numpy:
In [286]:
M[0,0] = "hola"
Si lo deseamos, podemos definir explícitamente el tipo de datos de un arreglo cuando lo creamos, usando el argumento dtype:
In [287]:
M = np.array([[1, 2], [3, 4]], dtype=complex)
M
Out[287]:
Algunos tipos comunes que pueden ser usados con dtype son: int, float, complex, bool, object, etc.
Podemos también definir explícitamente el número de bit de los tipos de datos, por ejemplo: int64, int16, float64, complex64.
In [288]:
# crea un arreglo con valores en un rango
x = np.arange(0, 10, 1) # argumentos: desde, hasta, paso
x
Out[288]:
In [289]:
x = np.arange(-1, 1, 0.1)
x
Out[289]:
In [290]:
# usando linspace, ambos puntos finales SON incluidos. Formato: (desde, hasta, número de elementos)
np.linspace(0, 10, 25)
Out[290]:
In [291]:
# logspace también incluye el punto final. Por defecto base=10
np.logspace(0, 10, 11, base=np.e) # produce e elevado a cada valor en linspace(0, 10, 11), e.d.[e**0, e**1,...,e**10]
Out[291]:
In [292]:
x, y = np.mgrid[0:5, 0:5] # similar a meshgrid en MATLAB
In [293]:
x
Out[293]:
In [294]:
y
Out[294]:
In [295]:
from numpy import random
In [296]:
# números aleatorios uniformes en [0,1]
random.rand(2,5)
Out[296]:
In [297]:
# números aleatorios con distribución estándar normal (distribución gaussiana de media 0 y varianza 1).
random.randn(2,5)
Out[297]:
In [298]:
# una matriz diagonal
np.diag([1,2,3])
Out[298]:
In [300]:
# diagonal desplazada desde la diagonal principal
np.diag([1,2,3], k=1)
Out[300]:
In [301]:
np.zeros((3,3))
Out[301]:
In [302]:
np.ones((3,3))
Out[302]:
Un formato muy común para archivos de datos es el de valores separados por comas, o formatos relacionados, como por ejemplo TSV (tab-separated values, valores separados por tabs). Para leer datos desde tales archivos a un arreglo Numpy podemos usar la función numpy.genfromtxt. Por ejemplo,
In [303]:
data = np.genfromtxt('stockholm_td_adj.dat') # asigna los datos del archivo al arreglo data
In [304]:
data.shape
Out[304]:
In [305]:
# ¿qué hace esta línea?. La respuesta mas adelante
%matplotlib inline
import matplotlib.pyplot as plt
In [306]:
fig, ax = plt.subplots(figsize=(14,4))
ax.plot(data[:,0]+data[:,1]/12.0+data[:,2]/365, data[:,5])
ax.axis('tight')
ax.set_title('temperaturas en Estocolmo')
ax.set_xlabel(u'año')
ax.set_ylabel(u'temperatura (°C)');
Usando numpy.savetxt podemos almacenar un arreglo Numpy a un archivo en formato CSV:
In [307]:
M = random.rand(3,3)
M
Out[307]:
In [308]:
np.savetxt("matriz-aleatoria.csv", M)
In [309]:
np.savetxt("matriz-aleatoria.csv", M, fmt='%.5f') # fmt especifica el formato
In [310]:
M.itemsize # los bits de cada elemento
Out[310]:
In [311]:
M.nbytes # número de bytes
Out[311]:
In [312]:
M.ndim # número de dimensiones
Out[312]:
In [313]:
# v es un vector, tiene por lo tanto sólo una dimensión, y requiere un índice
v[0]
Out[313]:
In [314]:
# M es una matriz, es decir un arreglo bidimensional, requiere dos índices
M[1,1]
Out[314]:
Si omitimos un índice de una arreglo multidimensional Numpy entrega la fila completa (o, en general, al arreglo de dimensión N-1 correspondiente)
In [315]:
M
Out[315]:
In [316]:
M[1]
Out[316]:
Puede obtenerse lo mismo usando : en el lugar de un índice:
In [317]:
M[1,:] # fila 1
Out[317]:
In [318]:
M[:,1] # columna 1
Out[318]:
Podemos asignar nuevos valores a los elementos de un arreglo usando el indexado:
In [319]:
M[0,0] = 1
In [320]:
M
Out[320]:
In [321]:
# también funciona para filas y columnas completas
M[1,:] = 0
M[:,2] = -1
In [322]:
M
Out[322]:
In [323]:
A = np.array([1,2,3,4,5])
A
Out[323]:
In [324]:
A[1:3]
Out[324]:
Los cortes de índices son mutables: si se les asigna un nuevo valor el arreglo original es modificado:
In [325]:
A[1:3] = [-2,-3]
A
Out[325]:
Podemos omitir cualquiera de los tres parámetros en M[desde:hasta:paso]:
In [326]:
A[::] # desde, hasta y paso asumen los valores por defecto
Out[326]:
In [327]:
A[::2] # el paso es 2, desde y hasta se asumen desde el comienzo hasta el fin del arreglo
Out[327]:
In [328]:
A[:3] # primeros tres elementos
Out[328]:
In [329]:
A[3:] # elementos desde el índice 3
Out[329]:
Los índices negativos se cuentan desde el fin del arreglo (los índices positivos desde el comienzo):
In [330]:
A = np.array([1,2,3,4,5])
In [331]:
A[-1] # el último elemento del arreglo
Out[331]:
In [332]:
A[-3:] # los últimos 3 elementos
Out[332]:
El corte de índices funciona exactamente del mismo modo para arreglos multidimensionales:
In [333]:
A = np.array([[n+m*10 for n in range(5)] for m in range(5)])
A
Out[333]:
In [334]:
# un bloque parte del arreglo original
A[1:4, 1:4]
Out[334]:
In [335]:
# elemento por medio
A[::2, ::2]
Out[335]:
In [336]:
indices_fila = [1, 2, 3]
A[indices_fila]
Out[336]:
In [337]:
indices_col = [1, 2, -1] # recuerde que el índice -1 corresponde al último elemento
A[indices_fila, indices_col]
Out[337]:
Podemos también usar máscaras de índices: Si la máscara de índice es un arreglo Numpy con tipo de dato booleano (bool), entonces un elemento es seleccionado (True) o no (False) dependiendo del valor de la máscara de índice en la posición de cada elemento:
In [338]:
B = np.array([n for n in range(5)])
B
Out[338]:
In [339]:
masc_fila = np.array([True, False, True, False, False])
B[masc_fila]
Out[339]:
In [340]:
# lo mismo
masc_fila = np.array([1,0,1,0,0], dtype=bool)
B[masc_fila]
Out[340]:
Esta característica es muy útil para seleccionar en forma condicional elementos de un arreglo, usando por ejemplo los operadores de comparación:
In [341]:
x = np.arange(0, 10, 0.5)
x
Out[341]:
In [342]:
masc = (5 < x) * (x < 7.5)
masc
Out[342]:
In [343]:
x[masc]
Out[343]:
In [344]:
indices = np.where(masc)
indices
Out[344]:
In [345]:
x[indices] # este indexado es equivalente al indexado fancy x[masc]
Out[345]:
In [346]:
np.diag(A)
Out[346]:
In [347]:
np.diag(A, -1)
Out[347]:
In [348]:
v2 = np.arange(-3,3)
v2
Out[348]:
In [349]:
indices_fila = [1, 3, 5]
v2[indices_fila] # indexado fancy
Out[349]:
In [350]:
v2.take(indices_fila)
Out[350]:
Pero la función take también funciona sobre listas y otros objetos:
In [351]:
np.take([-3, -2, -1, 0, 1, 2], indices_fila)
Out[351]:
In [352]:
cuales = [1, 0, 1, 0]
posibilidades = [[-2,-2,-2,-2], [5,5,5,5]]
np.choose(cuales, posibilidades)
Out[352]:
In [353]:
v1 = np.arange(0, 5)
In [354]:
2*v1
Out[354]:
In [355]:
v1 + 2
Out[355]:
In [356]:
A * 2, A + 2
Out[356]:
In [357]:
A * A # multiplicación elemento a elemento
Out[357]:
In [358]:
v1 * v1
Out[358]:
Si multiplicamos arreglos con formas compatibles, obtenemos una multiplicación elemento a elemento de cada fila:
In [359]:
A.shape, v1.shape
Out[359]:
In [360]:
A * v1
Out[360]:
In [361]:
np.dot(A, A)
Out[361]:
Desde Python 3.5 se pueden multiplicar matrices usando el operador @
In [362]:
A @ A
Out[362]:
In [363]:
A @ v1
Out[363]:
In [364]:
v1 @ v1
Out[364]:
Alternativamente, podemos transformar el arreglo al tipo matrix. Esto cambia el comportamiento de los operadores aritméticos estándar +, -, * al de álgebra de matrices.
In [365]:
M = np.matrix(A)
v = np.matrix(v1).T # aplica la traspuesta, convirtiéndolo en vector columna
In [366]:
v
Out[366]:
In [367]:
M*M
Out[367]:
In [368]:
M*v
Out[368]:
In [369]:
# productor interior
v.T * v
Out[369]:
In [370]:
# con objetos matriciales, el álgebra matricial estándar es usada
v + M*v
Out[370]:
Si intentamos sumar, restar, o multiplicar objetos con formas incompatibles, obtendremos un error:
In [371]:
v = np.matrix([1,2,3,4,5,6]).T
In [372]:
np.shape(M), np.shape(v)
Out[372]:
In [373]:
M * v
Vea también las funciones relacionadas: inner, outer, cross, kron, tensordot. Por ejemplo, introduzca help(kron).
Antes hemos usado .T para transponer un vector v. Podemos también usar la función transpose para conseguir el mismo resultado.
Otras funciones matemáticas que transforman objetos matriciales son:
In [374]:
C = np.matrix([[1j, 2j], [3j, 4j]])
C
Out[374]:
In [375]:
np.conjugate(C)
Out[375]:
Hermítico conjugado: transpuesta + conjugado
In [376]:
C.H
Out[376]:
Podemos extraer las partes reales e imaginarias de un arreglo con elementos complejos usando real y imag:
In [377]:
np.real(C) # lo mismo que: C.real
Out[377]:
In [378]:
np.imag(C) # lo mismo que: C.imag
Out[378]:
Podemos también extraer el módulo y el argumento complejo
In [379]:
np.angle(C+1) # Atención usuarios de MATLAB, se usa angle en lugar de arg
Out[379]:
In [380]:
abs(C)
Out[380]:
In [381]:
np.linalg.inv(C) # equivalente a C.I
Out[381]:
In [382]:
C.I * C
Out[382]:
In [383]:
np.linalg.det(C)
Out[383]:
In [384]:
np.linalg.det(C.I)
Out[384]:
In [385]:
# recuerde, los datos de la temperatura están almacenados en la variable data
np.shape(data)
Out[385]:
In [386]:
# la temperatura está almacenada en la columna 3
np.mean(data[:,3])
#data[:,3].mean()
Out[386]:
La temperatura diaria promedio en Estocolmo en los últimos 200 años ha sido aproximadamente 6.2 C.
In [387]:
np.std(data[:,3]), np.var(data[:,3])
Out[387]:
In [388]:
# valor mínimo del promedio diario de temperatura. Lo mismo que min(data[:,3])
data[:,3].min()
Out[388]:
In [389]:
# valor máximo del promedio diario de temperatura. Lo mismo que max(data[:,3])
data[:,3].max()
Out[389]:
In [390]:
d = 1 + np.arange(0, 10)
d
Out[390]:
In [391]:
# suma todos los elementos
np.sum(d)
Out[391]:
In [392]:
# multiplica todos los elementos
np.prod(d)
Out[392]:
In [393]:
# suma acumulativa
np.cumsum(d)
Out[393]:
In [394]:
# producto acumulativo
np.cumprod(d)
Out[394]:
In [395]:
# lo mismo que: diag(A).sum()
np.trace(A)
Out[395]:
In [396]:
!head -n 3 stockholm_td_adj.dat
El formato de los datos es: año, mes, día, temperatura promedio diaria, mínima, máxima, lugar.
Si estamos interesados sólo en la temperatura promedio de un mes particular, Febrero por ejemplo, podemos crear una máscara de índice y seleccionar sólo los datos de ese mes usando:
In [397]:
np.unique(data[:,1]) # la columna mes asume valores entre 1 y 12
Out[397]:
In [398]:
masc_feb = (data[:,1] == 2) # los paréntesis () son opcionales
masc_feb
Out[398]:
In [399]:
# los datos de temperatura están en la columna 3
np.mean(data[masc_feb,3])
Out[399]:
Estas funciones ponen a nuestra disposición herramientas muy poderosas para procesar datos. Por ejemplo, para extraer las temperaturas promedio mensuales para cada mes del año sólo necesitamos unas pocas líneas de código:
In [400]:
meses = np.arange(1,13)
media_mensual = [np.mean(data[data[:,1] == mes, 3]) for mes in meses]
fig, ax = plt.subplots()
ax.bar(meses, media_mensual)
ax.set_xlabel("Mes")
ax.set_ylabel("Temperatura promedio mensual");
In [401]:
m = random.rand(3,3)
m
Out[401]:
In [402]:
# máximo global
m.max()
Out[402]:
In [403]:
# máximo en cada columna
m.max(axis=0)
Out[403]:
In [404]:
# máximo en cada fila
m.max(axis=1)
Out[404]:
Muchas otras funciones y métodos de las clases array y matrix aceptan el argumento (optional) axis.
In [405]:
A
Out[405]:
In [406]:
n, m = A.shape
In [407]:
B = A.reshape((1,n*m))
B
Out[407]:
In [408]:
B[0,0:5] = 5 # modifica el arreglo
B
Out[408]:
In [409]:
A # la variable original es también cambiada. B sólo constituye una forma distinta de ver los mismos datos
Out[409]:
Podemos también usar la función flatten para transformar un arreglo multidimensional a un vector. A diferencia de reshape esta función crea una copia de los datos.
In [410]:
B = A.flatten()
B
Out[410]:
In [411]:
B[0:5] = 10
B
Out[411]:
In [412]:
A # ahora A no ha cambiado, ya que los datos de B han sido duplicados desde A
Out[412]:
In [413]:
v = np.array([1,2,3])
v
Out[413]:
In [414]:
np.shape(v)
Out[414]:
In [415]:
# crea una matriz con el vector v como su columna
v[:, np.newaxis]
Out[415]:
In [416]:
# matriz columna
v[:,np.newaxis].shape
Out[416]:
In [417]:
# matriz fila
v[np.newaxis,:]
Out[417]:
In [418]:
v[np.newaxis,:].shape
Out[418]:
In [419]:
a = np.array([[1, 2], [3, 4]])
a
Out[419]:
In [420]:
# repite cada elemente 3 veces
np.repeat(a, 3)
Out[420]:
In [421]:
# repite la matriz 3 veces
np.tile(a, 3)
Out[421]:
In [422]:
b = np.array([[5, 6]])
In [423]:
np.concatenate((a, b), axis=0)
Out[423]:
In [424]:
np.concatenate((a, b.T), axis=1)
Out[424]:
In [425]:
np.vstack((a,b))
Out[425]:
In [426]:
np.hstack((a,b.T))
Out[426]:
In [427]:
A = np.array([[1, 2], [3, 4]])
A
Out[427]:
In [428]:
# ahora B apunta al mismo arreglo que A
B = A
In [429]:
# cambiar B afecta a A
B[0,0] = 10
B
Out[429]:
In [430]:
A
Out[430]:
Si queremos evitar este comportamiento, para así obtener un nuevo objecto B copiado desde A, pero totalmente independiente de A, necesitamos realizar una "copia produnfa" ("deep copy") usando la función copy:
In [431]:
B = np.copy(A)
In [432]:
# ahora A no cambia si modificamos B
B[0,0] = -5
B
Out[432]:
In [433]:
A
Out[433]:
Generalmente, deseamos evitar iterar sobre los elementos de un arreglo donde sea posible (a cualquier precio!). La razón es que en un lenguaje interpretado como Python (o MATLAB), las iteraciones son realmente lentas comparadas con las operaciones vectorizadas.
Sin embargo, algunas veces es ineludible. En tales cases el bucle Python for es la forma más conveniente para iterar sobre un arreglo:
In [434]:
v = np.array([1,2,3,4])
for elemento in v:
print(elemento)
In [435]:
M = np.array([[1,2], [3,4]])
for fila in M:
print("fila", fila)
for elemento in fila:
print(elemento)
Cuando necesitamos iterar sobre cada elemento de un arreglo y modificar sus elementos, es conveniendo usar la función enumerate para obtener tanto el elemento como su índice en el bucle for:
In [436]:
for ind_fila, fila in enumerate(M):
print("ind_fila", ind_fila, "fila", fila)
for ind_colu, elemento in enumerate(fila):
print("ind_colu", ind_colu, "elemento", elemento)
# actualiza la matriz: eleva al cuadrado cada elemento
M[ind_fila, ind_colu] = elemento**2
In [437]:
# cada elemento en M está ahora al cuadrado
M
Out[437]:
Como se ha mencionado en varias ocasiones, para obtener un buen rendimiento deberíamos tratar de evitar realizar bucles sobre los elementos de nuestros vectores y matrices, y en su lugar usar algoritmos vectorizados. El primer paso para convertir un algoritmo escalar a uno vectorizado es asegurarnos que las funciones que escribamos funcionen con argumentos vectoriales.
In [438]:
def Theta(x):
"""
implementación escalar de la función escalón de Heaviside.
"""
if x >= 0:
return 1
else:
return 0
In [439]:
Theta(np.array([-3,-2,-1,0,1,2,3]))
Ok, eso no funcionó porque no definimos la función Theta de modo que pueda manejar argumentos vectoriales...
Para obtener una version vectorizada de Theta podemos usar la función vectorize de Numpy. En muchos casos, puede vectorizar automáticamente una función:
In [440]:
Theta_vec = np.vectorize(Theta)
In [441]:
Theta_vec(np.array([-3,-2,-1,0,1,2,3]))
Out[441]:
Podemos también implementar la función de modo que desde el comienzo acepte un argumento vectorial (esto requiere más esfuerzo, pero mejorar el rendimiento):
In [442]:
def Theta(x):
"""
Implementación preparada para vectores de la función escalón de Heaviside.
"""
return 1 * (x >= 0)
In [443]:
Theta(np.array([-3,-2,-1,0,1,2,3]))
Out[443]:
In [444]:
# y también funciona para escalares!
Theta(-1.2), Theta(2.6)
Out[444]:
In [445]:
M
Out[445]:
In [446]:
if (M > 5).any():
print("al menos un elemento de M es mayor que 5")
else:
print("ningún elemento de M es mayor que 5")
In [447]:
if (M > 5).all():
print("todos los elementos de M son mayores que 5")
else:
print("no todos los elementos de M son mayores que 5")
Como los arreglos de Numpy son de tipo estático, el tipo de un arrglo no puede ser cambiado luego de que es creado. Sin embarg, podemos convertir explícitamente un arreglo de un tipo a otro usando las funciones astype (ver también las funciones similares asarray). Esto crea un nuevo arreglo de un nuevo tipo:
In [448]:
M.dtype
Out[448]:
In [449]:
M2 = M.astype(float)
M2
Out[449]:
In [450]:
M2.dtype
Out[450]:
In [451]:
M3 = M.astype(bool)
M3
Out[451]:
In [452]:
# Esta celda da el estilo al notebook
from IPython.core.display import HTML
css_file = './css/aeropython.css'
HTML(open(css_file, "r").read())
Out[452]:
In [ ]: