NumPy

Ejecutar este documento en forma dinámica: Binder

Introducción

NumPy es un paquete de Scipy. Comentaremos brevemente cómo está compuesto Scipy para que se comprenda la potencialidad y el contexto de esta herramienta. Luego profundizaremos en NumPy.

Scipy es un acrónimo de "Scientific Computing Tools for Python" y se basa en los siguientes paquetes:

  • Python.

  • NumPy, diseñado para computación científica es su paquete principal para cálculo numérico. Define arreglos numéricos, matrices y operaciones básicas con éstos.

  • La librería Scipy, es una colección de algoritmos numéricos que incluyen: procesamiento de señales, optimización, estadística, procesamiento de imágenes, etc.

  • Matplotlib, paquete para visualización con calidad de publicación. Existen otros paquetes para visualización (una introducción en GMGcode).

También incluye otros paquetes especializados. Algunos de ellos son:

  • pandas, provee herramientas para el manejo sencillo de datos y estructuras de datos.

  • SymPy, para cálculo algebráico y matemática simbólica.

  • scikit-image es una colección de algoritmos para procesamiento de imágenes.

  • h5py manejo de datos en formato HDF5.

Ahora si, NumPy

Instalación

Primero se debe instalar, para esto ver las instrucciones aquí. Es recomendable hacerlo directamente de repositorios debian, ubuntu, etc. de tal manera que se instalen correctamente las dependencias. Como ya lo hemos mencionado, existen distribuciones de Python que instalan por defecto todas las herramientas necesarias (ejemplo: Anaconda).

Aspectos básicos

El objeto principal de NumPy es un arreglo multidimensional homogéneo. Es una tabla de elementos (números del mismo tipo) indexados como una tupla de enteros positivos. Por ejemplo:

In [ ]:
import numpy as np #importo numpy y lo denomino np
a = np.array([(1.5,2,3), (4,5,6)])
print(a)

Lo primero que hicimos fue importar el módulo NumPy como una variable de nombre np, esta es una variable del tipo módulo (<type 'module'>). A esta variable la usaremos para llamar todas las funciones y objetos que tiene NumPy. En el ejemplo anterior accedimos a la clase array (se llama ndarray) como np.array.

En NumPy se denomina axes a las dimensiones. En la variable a hemos guardado un arreglo que tiene dos ejes (axes), el primero tiene longitud 2 y el segundo 3.

La clase ndarray tiene varios atributos. Algunos de ellos son:

In [ ]:
#ndarray.ndim
ndimensiones = a.ndim
print('Tiene dos axes:',ndimensiones)

#ndarray.shape
forma = a.shape
print('Indica la dimensión de cada eje en una tupla:',forma)
print(type(forma))
forma[1]#Se puede acceder al valor de cada elemento de la tupla.

#ndarray.size
tamanio = a.size
print('Cantidad de elementos que puede almacenar:',tamanio)

#ndarray.dtype
print('Tipo de dato dentro del arreglo "a":',a.dtype)
b = np.array([0,1,2])
print('Tipo de dato dentro del arreglo "b":',b.dtype)
a[1][2]

Crear un arreglo

Existen varias formas de crear un arreglo NumPy:

  1. Como hicimos antes, Python list o tuplas.
  2. Usando funciones específicas de NumPy (ejemplos: arange, linspace, etc.).
  3. Leer datos desde un archivo.

1. Escribiendo los valores, desde listas o tuplas

Debe remarcarse que en los arreglos el tipo de elemento se determina cuando se crea el arreglo. En el ejemplo, a es float64 porque puse algún número con coma, y b es int64 porque los puse sin coma. También se puede especificar de la siguiente manera:

In [ ]:
c = np.array([(2.+1.j,1),(0+2.j,1-1j)],dtype=complex)
print('c =\n',c)
print('Parte real:\n',c.real)
print('Parte imaginaria:\n',c.imag)
c.imag.dtype

Algunos tipos de datos comunes son: int, float, complex, bool, object, etc. Además, podemos definir el tamaño (cantindad de bits) con: int64, int16, float128, complex128.

Otra cosa que se puede hacer es convertir una lista a un NumPy array:

In [ ]:
d = [2.1, 2.4, 2.3, 1.9]
dnp = np.asarray(d)
print('"d" es tipo:',type(d))
print('"dnp" es tipo:',type(dnp))
print(dnp.mean()) #método de la clase ndarray que calcula el promedio

Otro ejemplo convirtiendo una lista (esta puede servir cuando debo cargar muchos datos):

In [ ]:
e = []#inicializo una lista vacía
for n in range(0,30,5):
    e.append(n)
    
print(e)
print(type(e))
enp = np.asarray(e)
print(type(enp))
print(enp.shape)

2. Utilizando funciones de NumPy

Cuando debo escribir muchos números o inicializar un arreglo grande es conveniente usar funciones que tiene NumPy.

  • Ejemplos de creación de vectores:
In [ ]:
#Función arange (pariente de range), cuidado el límite superior no está incluido
x = np.arange(0.0, 10.0, 0.5) # argumentos: start, stop, step

#Función linspace, ambos límites están incluidos
y = np.linspace(0, 10, 25)

#Función logspace
pot1, pot2 = 2, 6
z = np.logspace(pot1, pot2, 5)#Aquí debo poner las potencias (10^pot1 10^pot2), puedo agregar también la base "base=e"
print(z)
  • Ejemplos de creación de matrices:

Matrices de todos ceros o todos unos:

In [ ]:
Ceros = np.zeros((2,2))
print(Ceros)
Unos = np.ones((3,3))
print(Unos)
Unosvec = np.ones(5)
print(Unosvec)

También se pueden crear matrices/vectores copiando el tamaño de una matriz/vector determinado. Por ejemplo, si quisieramos un arreglo del tamaño del vector de a:

In [ ]:
Ceros_como_a = np.zeros_like(a)
print(Ceros_como_a)
Unos_como_a = np.ones_like(a)
print(Unos_como_a)

Matrices diagonales:

In [ ]:
#ndarray.diag
A = np.diag([1,2,3])
print(A)

También se pueden crear otros tipos de matrices diagonales:

In [ ]:
np.diag([1,2,3], k=-1) 

También se pueden crear matrices o vectores aleatorios, para eso hay que utilizar la función random. Por ejemplo:

In [ ]:
Aruido = np.random.rand(5,5)
print(Aruido)

Note que en el ejemplo anterior llamamos la clase random que tiene un método llamado rand que genera números aleatorios en el intervalo (0,1) con una distribución aproximadamente Uniforme. El algoritmo utiliza el generador Mersenne Twister. Con una distribución Normal Estandar (media cero y varianza uno):

In [ ]:
k = np.random.randn(40)
print(k)
print('Media muestral:',k.mean())
print('Varianza muestral:',(k.std())**2)

como se ve, calculamos media y varianza muestrales y dan parecidas a las poblacionales.

3. Leer datos desde un archivo

Existen diferentes funciones para hacerlo:

a) Cargar un archivo de texto con loadtxt

b) Formato Numpy

a) Comenzaremos cargando un tipo de archivo del tipo .csv (comma separated value). En este tipo de archivo los datos están en columnas separadas por coma. La función de NumPy loadtxt carga directamente las columnas y las guarda un arreglos. Para ejemplificar tomaremos datos de un ensayo de tracción de una probeta de fundición de sección rectangular con dimensiones 4,8 mm $\times$ 16,9 mm.

In [ ]:
!head datos.dat #funciona en jupyter, en consola sólo head file.
F,L = np.loadtxt('datos.dat', usecols=(0,1),comments='#',delimiter=',', unpack=True)
epsilon = 100.*(L-L[0])/L[0] #Normalizando la deformación en %
A0 = 4.8e-3*16.9e-3 #Sección 
sigma = F/(A0*1e6) #Obteniendo el esfuerzo en MPa

Note que la función tiene como argumentos que columnas selecciono, con qué caracter comienzan los comentarios, cuál es el separador de los datos (en este caso coma). Se puede utilizar una función un poco más sofisticada que permite operar cuando faltan datos se llama genfromtxt. Podemos dibujar los datos utilizando matplotlib de la siguiente forma:

In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(6,6))
ax.plot(epsilon,sigma,'-o')
ax.axis('tight')
ax.set_title('Ensayo de tracción de pieza de fundición')
ax.set_xlabel(r'$\varepsilon$ %')
ax.set_ylabel(r'$\sigma$ [MPa]');

b) También se puede cargar un archivo del formato específico de NumPy. Para esto, primero vamos a salvar un archivo con ese formato y luego lo cargaremos.

In [ ]:
np.save("vectorAleatorio.npy", k)
!file vectorAleatorio.npy

Ahora lo cargamos:

In [ ]:
vecAleatorio = np.load("vectorAleatorio.npy")

También lo podemos graficar.

In [ ]:
fig, ax = plt.subplots(figsize=(6,6))
ax.plot(vecAleatorio,'-^',label= 'Distribución Normal')
ax.axis('tight')
ax.set_title('Vector aleatorio Normal')
ax.set_xlabel('Índice')
ax.set_ylabel('Valor')
plt.legend(framealpha=0.5,loc=2);

Manipulación de arreglos

Acceso mediante índices

Si tomamos el vector anterior, podemos acceder a cada elemento mediante corchetes y el índice, por ejemplo:

In [ ]:
vecAleatorio[8]

Si en cambio es una matriz como la que definimos antes M:

In [ ]:
M = np.random.rand(5,5)
print(M)
print(M[0,0])
print(M[0][0])#Es equivalente a la notación anterior
print(M[0,1])

O podríamos mostrar una fila:

In [ ]:
print(M[4])
M[4].shape

Una forma un poco más correcta es:

In [ ]:
M[4,:]#Fila

Para una columna:

In [ ]:
M[:,0]#Columna

Se pueden asignar valores a determinados elementos o conjunto de elementos (filas, columnas o parte de ellas).

In [ ]:
M[1,1] = 100.0#Al índice 1,1
print('M:\n',M)
M[:,0] = 0.0#Agrego ceros a toda la primer columna
print('M:\n',M)
M[2:4,3] = -10.0#Cambiar los elementos del tercero al quinto de la cuarta columna
print('M:\n',M)

Rebanado

Con vectores podemos hacer lo siguiente:

In [ ]:
a = np.arange(10)
print('Original a:',a)
a[2:9:3] # [start:end:step]
In [ ]:
a[:5]#El último no está incluido
In [ ]:
a[3:]
In [ ]:
a[-1]#Último elemento
In [ ]:
a[::-1]#Invierto
In [ ]:
a[::-2]#Invierto y remuestreo
In [ ]:
a[-3:]#Últimos tres elementos

Para matrices, más o menos similar:

In [ ]:
A = np.array([[n+m*10 for n in range(5)] for m in range(5)])
A
In [ ]:
A[::2, ::2]

Indexado Fancy

Se denomina indexado fancy cuando un arreglo o una lista se utiliza en lugar de un índice. Algunos ejemplos de este tipo de indexado:

In [ ]:
indices_filas = [2, 3]#Tercer y cuarta fila
A[indices_filas]
In [ ]:
indices_columnas= [1, -1]#Segunda columna y última
A[:,indices_columnas]

También se puede enmascarar algunos elementos con variables booleanas:

In [ ]:
b = np.array([n for n in range(5)])
b
In [ ]:
mascara = np.array([True, False, True, False, False])
b[mascara]

Mejor que pasarle a mano estas variables se pueden utilizar inecuaciones:

In [ ]:
mascara = b < 2
b[mascara]
In [ ]:
c = np.array([n for n in range(20)])
print(c)
mascara = (c < 10) & (c % 2 == 0) #Menores que diez divisibles por dos
c[mascara]

Operaciones por elemento

Aritméticas

Suma con un escalar:

In [ ]:
A+100

Suma de elemento a elemento:

In [ ]:
A+A

Producto elemento a elemento (cuidado no es producto de matrices):

In [ ]:
A*A

Para multiplicar matrices se usar matmul):

In [ ]:
B = np.eye(2)
print(B)
C = np.ones((2, 2)) * 2
print(C)
np.matmul(B,C)

Para comparar elementos se pueden utilizar <, >, ==, etc.:

In [ ]:
a = np.array([1, 2, 3, 4])
b = np.array([4, 2, 2, 4])
a==b

También se puede comparar el arreglo completo:

In [ ]:
np.array_equal(a,b)
In [ ]:
np.array_equal(a,a)

También existen otras operaciones lógicas como logical_or y logical_and para esto los arreglos deben ser booleanos:

In [ ]:
a = np.array([1, 1, 0, 0], dtype=bool)
b = np.array([1, 0, 1, 0], dtype=bool)
np.logical_or(a, b)

Funciones

Existen varias funciones implementadas en numpy. Trigonométricas, estadísticas (ya vimos mean(), std()), etc.

In [ ]:
t = np.linspace(0., 10., 100)
f = 0.5 #Frecuencia 0.5 Hz
y = np.sin(2*np.pi*f*t)
fig, ax = plt.subplots(figsize=(6,6))
ax.plot(t,y)
ax.axis('tight')
ax.set_title('Función trigonométricas')
ax.set_xlabel('tiempo [s]')
ax.set_ylabel('y');

Algunas operaciones conocidas como reducciones:

In [ ]:
x = np.array([1, 2, 3, 4])
x.sum()#también funciona np.sum(x)

Se pueden sumar los diferentes ejes:

In [ ]:
a = np.array([[1.0, 2.0], [3.0, 4.0]])
print('Matriz original:\n',a)
print('Por columnas:\n',a.sum(axis=0))
print('Por filas:\n',a.sum(axis=1))

Otro tipo de funciones y reducciones (también se pueden hacer por eje):

In [ ]:
a.min(axis=0)
In [ ]:
x = np.array([10, 20, 8, 9])
x.argmin()

Operaciones lógicas y comparaciones:

In [ ]:
np.any(x==8)
In [ ]:
np.all(x==8)

Se puede manipular la forma de un arreglo. Aquí van unas funciones útiles.

Aplanamiento:

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

Reshape:

In [ ]:
b = a.ravel()
b = b.reshape((2, 3))
b

Algunas cosas más raras:

In [ ]:
a = np.arange(4*3*2).reshape(4, 3, 2)
a

Se puede ordenar:

In [ ]:
a = np.array([[4, 3, 5], [1, 2, 1]])
b = np.sort(a, axis=1)
b

Álgebra lineal

También se pueden realizar algunos cálculos de álgebra lineal. Para esto se deben acceder a algunas funciones de np.linalg. Aquí un ejemplo:

$Ax = b$

In [ ]:
A = np.array([[1.0, 2.0], [3.0, 4.0]])
print('Matriz original:\n',A)
print('Matriz transpuesta:\n',A.transpose())#También funciona A.T
print('Inversa de la Matriz:\n',np.linalg.inv(A))
print('Verificamos:\n',np.matmul(A,np.linalg.inv(A)))
b = np.array([[5.], [7.]])
print('Solución del sistema lineal (despejamos x):\n',np.linalg.solve(A, b))

Histogramas

Se pueden construir histogramas muy facilmente con ayuda de matplotlib. Un ejemplo de un histograma de una variable aleatoria con distribución Normal.

In [ ]:
mu, sigma = 2, 0.5
v = np.random.normal(mu,sigma,10000)
plt.hist(v, bins=50, density=1);
In [ ]:
(n, bins) = np.histogram(v, bins=50, density=True)  # Con la versión de NumPy
plt.plot(.5*(bins[1:]+bins[:-1]), n);

Interpolación (fitting)

Con NumPy también podemos realizar una interpolación con polinomios y otras funciones. Continuaremos con el ejemplo del ensayo de tracción y estimaremos el módulo de elasticidad. Utilizaremos la función polyfit.

In [ ]:
F,L = np.loadtxt('datos.dat', usecols=(0,1),comments='#',delimiter=',', unpack=True)
epsilon = 100.*(L-L[0])/L[0] #Normalizando la deformación en %
A0 = 4.8e-3*16.9e-3 #Sección 
sigma = F/(A0*1e6) #Obteniendo el esfuerzo en MPa
fig, ax = plt.subplots(figsize=(6,6))
ax.plot(epsilon,sigma,'o')
ax.axis('tight')
ax.set_title('Ensayo de tracción de pieza de fundición')
ax.set_xlabel(r'$\varepsilon$ %')
ax.set_ylabel(r'$\sigma$ [MPa]')
p = np.polyfit(epsilon[:5], sigma[:5], 1)
sigma_est = np.polyval(p, np.arange(-0.5,.2,0.01))
ax.plot(np.arange(-0.5,.2,0.01),sigma_est,'r');
ax.set_xlim(0,1)
print('Módulo de elasticidad E =',p[1]/(100*1e-3),' GPa')
In [ ]: