USM Numérica

Librería numpy

Objetivos

Conocer la librería numpy y su aplicación para cálculo numérico.
Aprender las diferencias de utilización entre numpy.matrix y numpy.array.

0.1 Instrucciones

Las instrucciones de instalación y uso de un ipython notebook se encuentran en el siguiente link.

Después de descargar y abrir el presente notebook, recuerden:

Desarrollar los problemas de manera secuencial.
Guardar constantemente con Ctr-S para evitar sorpresas.
Reemplazar en las celdas de código donde diga FIX_ME por el código correspondiente.
Ejecutar cada celda de código utilizando Ctr-Enter

0.2 Licenciamiento y Configuración

Ejecutar la siguiente celda mediante Ctr-S.



In [2]:

    
"""
IPython Notebook v4.0 para python 3.0
Librerías adicionales: numpy, matplotlib
Contenido bajo licencia CC-BY 4.0. Código bajo licencia MIT. 
(c) Sebastian Flores, Christopher Cooper, Alberto Rubio, Pablo Bunout.
"""
# Configuración para recargar módulos y librerías dinámicamente
%reload_ext autoreload
%autoreload 2

# Configuración para graficos en línea
%matplotlib inline

# Configuración de estilo
from IPython.core.display import HTML
HTML(open("./style/style.css", "r").read())









    Out[2]:

Contenido

Overview de Numpy y Scipy
Librería Numpy
- Arreglos vs Matrices
- Axis
- Funciones basicas.
- Input y Output
Tips

1. Overview de numpy y scipy

¿Cual es la diferencia entre numpy y scipy?

In an ideal world, NumPy would contain nothing but the array data type and the most basic operations: indexing, sorting, reshaping, basic elementwise functions, et cetera. All numerical code would reside in SciPy. However, one of NumPy’s important goals is compatibility, so NumPy tries to retain all features supported by either of its predecessors. Thus NumPy contains some linear algebra functions, even though these more properly belong in SciPy. In any case, SciPy contains more fully-featured versions of the linear algebra modules, as well as many other numerical algorithms. If you are doing scientific computing with python, you should probably install both NumPy and SciPy. Most new features belong in SciPy rather than NumPy. Link stackoverflow

Por ser python un lenguaje open-source, existen miles de paquetes disponibles creados por individuos o comunidades. Éstos pueden estar disponibles en un repositorio como github o bitbucket, o bien estar disponibles en el repositorio oficial de python: pypi. En un inicio, cuando no existía una librerías de cálculo científico oficial, varios candidatos proponían soluciones:

numpy: tenía una excelente representación de vectores, matrices y arreglos, implementados en C y llamados fácilmente desde python
scipy: proponía linkear a librerías ya elaboradas de calculo científico de alto rendimiento en C o fortran, permitiendo ejecutar rápidamente desde python.

Ambos projectos fueron creciendo en complejidad y alcance, y en vez de competir, decidieron dividir tareas y unificar fuerzas para proponer una plataforma de cálculo científico que reemplazara completamente otros programas.

numpy: Corresponde a lo relacionado con la estructura de los datos (arrays densos y sparse, matrices, constructores especiales, lectura de datos regulares, etc.), pero no las operaciones en sí. Por razones históricas y de compatibilidad, tiene algunos algoritmos, pero en realidad resulta más consistente utilizar los algoritmos de scipy.
scipy: Corresponde a la implementación numérica de diversos algoritmos de corte científicos: algebra lineal, estadística, ecuaciones diferenciales ordinarias, interpolacion, integracion, optimización, análisis de señales, entre otros.

OBSERVACIÓN IMPORTANTE:

Las matrices y arrays de numpy deben contener variables con el mismo tipo de datos: sólo enteros, sólo flotantes, sólo complejos, sólo booleanos o sólo strings. La uniformicidad de los datos es lo que permite acelerar los cálculos con implementaciones en C a bajo nivel.

2. Librería Numpy

Siempre importaremos la librería numpy de la siguiente forma:

import numpy as np

Todas las funciones y módulos de numpy quedan a nuestro alcance a 3 carácteres de distancia:

np.array([1,4,9,16])
np.linspace(0.,1.,100)

Evite a todo costo utilizar lo siguiente:

from numpy import *



In [ ]:

    
import numpy as np
print np.version.version # Si alguna vez tienen problemas, verifiquen su version de numpy

Importante

Ipython notebook es interactivo y permite la utilización de tabulación para ofrecer sugerencias o enseñar ayuda (no solo para numpy, sino que para cualquier código en python).

Pruebe los siguientes ejemplos:



In [ ]:

    
# Presionar tabulacción con el cursor despues de np.arr
np.arr



In [ ]:

    
# Presionar Ctr-Enter para obtener la documentacion de la funcion np.array usando "?"
np.array?



In [ ]:

    
# Presionar Ctr-Enter
%who



In [ ]:

    
x = 10
%who

2. Librería Numpy

2.1 Array vs Matrix

Por defecto, la gran mayoria de las funciones de numpy y de scipy asumen que se les pasará un objeto de tipo array.

Veremos las diferencias entre los objetos array y matrix, pero recuerden utilizar array mientras sea posible.

Matrix

Una matrix de numpy se comporta exactamente como esperaríamos de una matriz:

Pros:

Multiplicación utiliza el signo * como es esperable.
Resulta natural si lo único que haremos es algebra lineal.

Contras:

Todas las matrices deben estar completamente alineadas para poder operar correctamente.
Operaciones elementwise son mas dificiles de definir/acceder.
Están exclusivamente definidas en 2D: un vector fila o un vector columna siguen siendo 2D.



In [ ]:

    
# Operaciones con np.matrix
A = np.matrix([[1,2],[3,4]])
B = np.matrix([[1, 1],[0,1]], dtype=float)
x = np.matrix([[1],[2]])
print "A =\n", A
print "B =\n", B
print "x =\n", x

print "A+B =\n", A+B
print "A*B =\n", A*B
print "A*x =\n", A*x
print "A*A = A^2 =\n", A**2
print "x.T*A =\n", x.T * A

2.1 Array vs Matrix

Array

Un array de numpy es simplemente un "contenedor" multidimensional.

Pros:

Es multidimensional: 1D, 2D, 3D, ...
Resulta consistente: todas las operaciones son element-wise a menos que se utilice una función específica.

Contras:

Multiplicación maticial utiliza la función dot()



In [ ]:

    
# Operaciones con np.matrix
A = np.array([[1,2],[3,4]])
B = np.array([[1, 1],[0,1]], dtype=float)
x = np.array([1,2]) # No hay necesidad de definir como fila o columna!
print "A =\n", A
print "B =\n", B
print "x =\n", x

print "A+B =\n", A+B
print "AoB = (multiplicacion elementwise) \n", A*B
print "A*B = (multiplicacion matricial, v1) \n", np.dot(A,B)
print "A*B = (multiplicacion matricial, v2) \n", A.dot(B)
print "A*A = A^2 = (potencia matricial)\n", np.linalg.matrix_power(A,2)
print "AoA = (potencia elementwise)\n", A**2
print "A*x =\n", np.dot(A,x)
print "x.T*A =\n", np.dot(x,A) # No es necesario transponer.

Desafío 1: matrix vs array

Sean $$ A = \begin{pmatrix} 1 & 0 & 1 \\ 0 & 1 & 1\end{pmatrix} $$ y $$ B = \begin{pmatrix} 1 & 0 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 1\end{pmatrix} $$

Cree las matrices utilizando np.matrix y multipliquelas en el sentido matricial. Imprima el resultado.
Cree las matrices utilizando np.array y multipliquelas en el sentido matricial. Imprima el resultado.



In [ ]:

    
# 1: Utilizando matrix
A = np.matrix([]) # FIX ME
B = np.matrix([]) # FIX ME
print "np.matrix, AxB=\n", #FIX ME

# 2: Utilizando arrays
A = np.array([]) # FIX ME
B = np.array([]) # FIX ME
print "np.matrix, AxB=\n", #FIX ME

2.2 Indexación y Slicing

Los arrays se indexan de la forma "tradicional".

Para un array unidimensional: sólo tiene una indexacion. ¡No es ni fila ni columna!
Para un array bidimensional: primera componente son las filas, segunda componente son las columnas. Notación respeta por tanto la convención tradicional de matrices.
Para un array tridimensional: primera componente son las filas, segunda componente son las columnas, tercera componente la siguiente dimension.

Respecto a los índices de los elementos, éstos comienzan en cero, como en C. Además, es posible utilizar índices negativos, que como convención asignan -1 al último elemento, al -2 el penúltimo elemento, y así sucesivamente.

Por ejemplo, si a = [2,3,5,7,11,13,17,19], entonces a[0] es el valor 2 y a[1] es el valor 3, mientras que a[-1] es el valor 19 y a[-2] es el valor 17.

Ademas, en python existe la "slicing notation":

a[start:end] : items desde índice start hasta end-1
a[start:] : items desde índice start hasta el final del array
a[:end] : items desde el inicio hasta el índice end-1
a[:] : todos los items del array (una copia nueva)
a[start:end:step] : items desde start hasta pasar end (sin incluir) con paso step



In [ ]:

    
x = np.arange(9)   # "Vector" con valores del 0 al 8 
print "x = ", x
print "x[:] = ", x[:]
print "x[5:] = ", x[5:]
print "x[:8] = ", x[:8]
print "x[:-1] = ", x[:-1]
print "x[1:-1] = ", x[1:-1]
print "x[1:-1:2] = ", x[1:-1:2]

A = x.reshape(3,3) # Arreglo con valores del 0 al 8, en 3 filas y 3 columnas.
print "\n"
print "A = \n", A
print "primera fila de A\n", A[0,:]
print "ultima columna de A\n", A[:,-1]
print "submatriz de A\n", A[:2,:2]

Observación

Cabe destacar que al tomar slices (subsecciones) de un arreglo obtenemos siempre un arreglo con menores dimensiones que el original.
Esta notación es extremadamente conveniente, puesto que nos permite manipular el array sin necesitar conocer el tamaño del array y escribir de manera compacta las fórmulas numéricas.

Por ejemplo, implementar una derivada numérica es tan simple como sigue.



In [ ]:

    
def f(x):
    return 1 + x**2

x = np.array([0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]) # O utilizar np.linspace!
y = f(x) # Tan facil como llamar f sobre x
dydx = ( y[1:] - y[:-1] ) / ( x[1:] - x[:-1] )
x_aux = 0.5*(x[1:] + x[:-1])
# To plot
fig = plt.figure(figsize=(12,8))
plt.plot(x, y, '-s', label="f")
plt.plot(x_aux, dydx, '-s', label="df/dx")
plt.legend(loc="upper left")
plt.show()

Desafío 2: Derivación numérica

Implemente el cálculo de la segunda derivada, que puede obtenerse por diferencias finitas centradas mediante $$ \frac{d f(x_i)}{dx} = \frac{1}{\Delta x^2} \Big( f(x_{i+1}) -2 f(x_{i}) + f(x_{i-1}) \Big)$$



In [ ]:

    
def g(x):
    return 1 + x**2 + np.sin(x)

x = np.linspace(0,1,10)
y = g(x) 
d2ydx2 = 0 * x # FIX ME
x_aux = 0*d2ydx2 # FIX ME
# To plot
fig = plt.figure(figsize=(12,8))
plt.plot(x, y, label="f")
plt.plot(x_aux, d2ydx2, label="d2f/dx2")
plt.legend(loc="upper left")
plt.show()

2. Librería Numpy

2.2 Funciones Básicas

Algunas funciones básicas que es conveniente conocer son las siguientes:

shape: Entrega las dimensiones del arreglo. Siempre es una tupla.
len: Entrega el número de elementos de la primera dimensión del arreglo. Siempre es un entero.
ones: Crea un arreglo con las dimensiones provistas e inicializado con valores 1. Por defecto array 1D.
zeros: Crea un arreglo con las dimensiones provistas e inicializado con valores 1. Por defecto array 1D.
eye: Crea un arreglo con las dimensiones provistas e inicializado con 1 en la diagonal. Por defecto array 2D.



In [ ]:

    
# arrays 1d
A = np.ones(3)
print "A = \n", A
print "A.shape =", A.shape
print "len(A) =", len(A)
B = np.zeros(3)
print "B = \n", B
print "B.shape =", B.shape
print "len(B) =", len(B)
C = np.eye(1,3)
print "C = \n", C
print "C.shape =", C.shape
print "len(C) =", len(C)
# Si queremos forzar la misma forma que A y B
C = np.eye(1,3).flatten() # o np.eye(1,3)[0,:] 
print "C = \n", C
print "C.shape =", C.shape
print "len(C) =", len(C)



In [ ]:

    
# square arrays
A = np.ones((3,3))
print "A = \n", A
print "A.shape =", A.shape
print "len(A) =", len(A)
B = np.zeros((3,3))
print "B = \n", B
print "B.shape =", B.shape
print "len(B) =", len(B)
C = np.eye(3) # Or  np.eye(3,3)
print "C = \n", C
print "C.shape =", C.shape
print "len(C) =", len(C)



In [ ]:

    
# fat 2d array
A = np.ones((2,5))
print "A = \n", A
print "A.shape =", A.shape
print "len(A) =", len(A)
B = np.zeros((2,5))
print "B = \n", B
print "B.shape =", B.shape
print "len(B) =", len(B)
C = np.eye(2,5)
print "C = \n", C
print "C.shape =", C.shape
print "len(C) =", len(C)

2. Librería Numpy

2.2 Funciones Básicas

Algunas funciones básicas que es conveniente conocer son las siguientes:

reshape: Convierte arreglo a nueva forma. Numero de elementos debe ser el mismo.
linspace: Regresa un arreglo con valores linealmente espaciados.
diag(x): Si x es 1D, regresa array 2D con valores en diagonal. Si x es 2D, regresa valores en la diagonal.
sum: Suma los valores del arreglo. Puede hacerse en general o a lo largo de un axis.
mean: Calcula el promedio de los valores del arreglo. Puede hacerse en general o a lo largo de un axis.
std: Calcula la desviación estándar de los valores del arreglo. Puede hacerse en general o a lo largo de un axis.



In [ ]:

    
x = np.linspace(0., 1., 6)
A = x.reshape(3,2)
print "x = \n", x
print "A = \n", A

print "np.diag(x) = \n", np.diag(x)
print "np.diag(B) = \n", np.diag(A)

print ""
print "A.sum() = ", A.sum()
print "A.sum(axis=0) = ", A.sum(axis=0)
print "A.sum(axis=1) = ", A.sum(axis=1)

print ""
print "A.mean() = ", A.mean()
print "A.mean(axis=0) = ", A.mean(axis=0)
print "A.mean(axis=1) = ", A.mean(axis=1)

print ""
print "A.std() = ", A.std()
print "A.std(axis=0) = ", A.std(axis=0)
print "A.std(axis=1) = ", A.std(axis=1)

Desafío 3

Complete el siguiente código:

Se le provee un array A cuadrado
Calcule un array B como la multiplicación element-wise de A por sí misma.
Calcule un array C como la multiplicación matricial de A y B.
Imprima la matriz C resultante.
Calcule la suma, promedio y desviación estándar de los valores en la diagonal de C.
Imprima los valores anteriormente calculados.



In [ ]:

    
A = np.outer(np.arange(3),np.arange(3))
print A
# FIX ME
# FIX ME
# FIX ME
# FIX ME
# FIX ME

Desafío 4

Implemente la regla de integración trapezoidal



In [ ]:

    
def mi_funcion(x):
    f = 1 + x + x**3 + x**5 + np.sin(x)
    return f

N = 5
x = np.linspace(-1,1,N)
y = mi_funcion(x)
# FIX ME
I = 0 # FIX ME
# FIX ME
print "Area bajo la curva: %.3f" %I

# Ilustración gráfica
x_aux = np.linspace(x.min(),x.max(),N**2)
fig = plt.figure(figsize=(12,8))
fig.gca().fill_between(x, 0, y, alpha=0.25)
plt.plot(x_aux, mi_funcion(x_aux), 'k')
plt.plot(x, y, 'r.-')
plt.show()

2. Librería Numpy

2.5 Inputs y Outputs

Numpy permite leer datos en formato array con la función loadtxt. Existen variados argumentos opcionales, pero los mas importantes son:

skiprows: permite saltarse lineas en la lectura.
dtype: declarar que tipo de dato tendra el array resultante



In [ ]:

    
# Ejemplo de lectura de datos
data = np.loadtxt("data/cherry.txt")
print data.shape
print data



In [ ]:

    
# Ejemplo de lectura de datos, saltandose 11 lineas y truncando a enteros
data_int = np.loadtxt("data/cherry.txt", skiprows=11).astype(int)
print data_int.shape
print data_int

2. Librería Numpy

2.5 Inputs y Outputs

Numpy permite guardar datos de manera sencilla con la función savetxt: siempre debemos dar el nombre del archivo y el array a guardar.

Existen variados argumentos opcionales, pero los mas importantes son:

header: Línea a escribir como encabezado de los datos
fmt: Formato con el cual se guardan los datos (%d para enteros, %.5f para flotantes con 5 decimales, %.3E para notación científica con 3 decimales, etc.).



In [ ]:

    
# Guardando el archivo con un header en español
encabezado = "Diametro Altura Volumen (Valores truncados a numeros enteros)"
np.savetxt("data/cherry_int.txt", data_int, fmt="%d", header=encabezado)

Revisemos si el archivo quedó bien escrito. Cambiaremos de python a bash para utilizar los comandos del terminal:



In [ ]:

    
%%bash 
cat data/cherry_int.txt

Desafío 5

Lea el archivo data/cherry.txt
Escale la matriz para tener todas las unidades en metros o metros cubicos.
Guarde la matriz en un nuevo archivo data/cherry_mks.txt, con un encabezado apropiado y 2 decimales de precisión para el flotante (pero no en notación científica).



In [ ]:

    
# Leer datos
#FIX_ME#

# Convertir a mks
#FIX_ME#

# Guardar en nuevo archivo
#FIX_ME#

2. Librería Numpy

2.6 Selecciones de datos

Existen 2 formas de seleccionar datos en un array A:

Utilizar máscaras de datos, que corresponden a un array con las mismas dimensiones del array A, pero de tipo booleano. Todos aquellos elementos True del array de la mascara serán seleccionados.
Utilizar un array con valores enteros. Los valores del array indican los valores que desean conservarse.

2.6 Máscaras

Observe que el array regresado siempre es unidimensional puesto que no es posible garantizar que se mantenga la dimensión original del array.



In [ ]:

    
x = np.linspace(0,42,10)
print "x = ", x
print "x.shape = ", x.shape

print "\n"
mask_x_1 = x>10
print "mask_x_1 = ", mask_x_1
print "x[mask_x_1] = ", x[mask_x_1]
print "x[mask_x_1].shape = ", x[mask_x_1].shape

print "\n"
mask_x_2 = x > x.mean()
print "mask_x_2 = ", mask_x_2
print "x[mask_x_2] = ", x[mask_x_2]
print "x[mask_x_2].shape = ", x[mask_x_2].shape



In [ ]:

    
A = np.linspace(10,20,12).reshape(3,4)
print "\n"
print "A = ", A
print "A.shape = ", A.shape

print "\n"
mask_A_1 = A>13
print "mask_A_1 = ", mask_A_1
print "A[mask_A_1] = ", A[mask_A_1]
print "A[mask_A_1].shape = ", A[mask_A_1].shape

print "\n"
mask_A_2 = A > 0.5*(A.min()+A.max())
print "mask_A_2 = ", mask_A_2
print "A[mask_A_2] = ", A[mask_A_2]
print "A[mask_A_2].shape = ", A[mask_A_2].shape



In [ ]:

    
T = np.linspace(-100,100,24).reshape(2,3,4)
print "\n"
print "T = ", T
print "T.shape = ", T.shape

print "\n"
mask_T_1 = T>=0
print "mask_T_1 = ", mask_T_1
print "T[mask_T_1] = ", T[mask_T_1]
print "T[mask_T_1].shape = ", T[mask_T_1].shape

print "\n"
mask_T_2 =  1 - T + 2*T**2 < 0.1*T**3
print "mask_T_2 = ", mask_T_2
print "T[mask_T_2] = ", T[mask_T_2]
print "T[mask_T_2].shape = ", T[mask_T_2].shape

2.6 Índices

Observe que es posible repetir índices, por lo que el array obtenido puede tener más elementos que el array original.

En un arreglo 2d, es necesario pasar 2 arrays, el primero para las filas y el segundo para las columnas.



In [ ]:

    
x = np.linspace(10,20,11)
print "x = ", x
print "x.shape = ", x.shape

print "\n"
ind_x_1 = np.array([1,2,3,5,7])
print "ind_x_1 = ", ind_x_1
print "x[ind_x_1] = ", x[ind_x_1]
print "x[ind_x_1].shape = ", x[ind_x_1].shape

print "\n"
ind_x_2 = np.array([0,0,1,2,3,4,5,6,7,-3,-2,-1,-1])
print "ind_x_2 = ", ind_x_2
print "x[ind_x_2] = ", x[ind_x_2]
print "x[ind_x_2].shape = ", x[ind_x_2].shape



In [ ]:

    
A = np.linspace(-90,90,10).reshape(2,5)
print "A = ", A
print "A.shape = ", A.shape

print "\n"
ind_row_A_1 = np.array([0,0,0,1,1])
ind_col_A_1 = np.array([0,2,4,1,3])
print "ind_row_A_1 = ", ind_row_A_1
print "ind_col_A_1 = ", ind_col_A_1
print "A[ind_row_A_1,ind_col_A_1] = ", A[ind_row_A_1,ind_col_A_1]
print "A[ind_row_A_1,ind_col_A_1].shape = ", A[ind_row_A_1,ind_col_A_1].shape

print "\n"
ind_row_A_2 = 1
ind_col_A_2 = np.array([0,1,3])
print "ind_row_A_2 = ", ind_row_A_2
print "ind_col_A_2 = ", ind_col_A_2
print "A[ind_row_A_2,ind_col_A_2] = ", A[ind_row_A_2,ind_col_A_2]
print "A[ind_row_A_2,ind_col_A_2].shape = ", A[ind_row_A_2,ind_col_A_2].shape

Desafío 6

La potencia de un aerogenerador, para $k$ una constante relacionada con la geometría y la eficiencia, $\rho$ la densidad del are, $r$ el radio del aerogenerador en metros y $v$ la velocidad el viento en metros por segundo, está dada por: $$ P = \begin{cases} k \ \rho \ r^2 \ v^3, 3 \leq v \leq 25\\ 0,\ eoc\end{cases}$$ Típicamente se considera una valor de $k=0.8$ y una densidad del aire de $\rho = 1.2$ [$kg/m^3$].

Calcule el número de aerogeneradores activos, la potencia promedio y la potencia total generada por los 11 generadores del parque Eólico Canela 1.

Los valores de radio del aerogenerador (en metros) y la velocidad del viento (en kilometros por hora) se indican a continuación en arreglos en el código numérico.



In [ ]:

    
import numpy as np
k = 0.8
rho = 1.2 # 
r_m = np.array([ 25.,  25.,  25.,  25.,  25.,  25.,  20.,  20.,  20.,  20.,  20.])
v_kmh = np.array([10.4, 12.6,  9.7,  7.2, 12.3, 10.8, 12.9, 13.0,  8.6, 12.6, 11.2]) # En kilometros por hora
P = 0 
n_activos = 0
P_mean = 0.0
P_total =  0.0
print "Existen %d aerogeneradores activos del total de %d" %(n_activos, r.shape[0])
print "La potencia promedio de los aeorgeneradores es {0:.2f} ".format(P_mean) 
print "La potencia promedio de los aeorgeneradores es " + str(P_total)

Tips

La práctica y la necesidad hace al maestro.
Preguntar: en línea y en persona, pero tratar de solucionar los problemas antes.

Enlaces útiles:

http://www.labri.fr/perso/nrougier/teaching/numpy.100/ : Lista con 100 recetas prácticas.s
http://sebastianraschka.com/Articles/2014_matlab_vs_numpy.html : Desde matlab a numpy.
http://pages.physics.cornell.edu/~myers/teaching/ComputationalMethods/python/arrays.html : Otros consejos.