Arreglos de NumPy¶
Numpy es el paquete fundamental para cómputo científico en Python. Se trata de una librería para Python que proporciona un arreglo multidimensional, varios objetos derivados (tales como arreglos enmascarados y matrices), y varias rutinas para operaciones rápidas sobre arreglos, incluyendo manipulación matemática, lógica y de la “forma” de los arreglos; ordenamiento, selección, I/O, transformadas discretas de Fourier, álgebra lineal básica, operaciones estadísticas básicas, simulación aleatoria y mucho más.
Para ilustrar la importancia de los arreglos de NumPy, considere el siguiente problema: multiplicar cada elemento de una secuencia unidimensional con el correspondiente elemento de otra secuencia de la misma longitud. Si los datos están almacenados en dos listas de Python, a y b, podríamos iterar sobre cada elemento, de la siguiente manera:
c = []
for i in range(len(a)):
c.append(a[i]*b[i])
Esto produce el resultado correcto, pero si a y b contienen cada uno un millón de números, pagaremos el precio por las ineficiencias de los ciclos en Python. Podríamos conseguir la misma tarea mucho más rápido en C al escribir (por simplicidad no escribimos las declaraciones de variables, inicializaciones, la reservación del espacio de memoria, etc.):
for (i = 0; i < numero_elementos; i++) {
c[i] = a[i]*b[i];
}
Utilizando arreglos de Numpy a y b, podemos realizar la misma operación con:
c = a * b
casi a la misma velocidad que con C, pero con la simplicidad en la programación que se espera de código escrito en Python.
Un ejemplo más realista en C¶
El siguiente ejemplo muestra el producto punto de dos vectores en C (se está utilizando el framework “Accelerate”, específico para Mac OS X, que incluye una librería para BLAS, es decir “Basic Linear Algebra”):
#include <stdio.h>
#include <stdlib.h>
#include <Accelerate/Accelerate.h>
void runSDOTSample (void) {
float *X, *Y;
int N, i;
float resultFloat;
// Initialize the inputs and allocate space
N = 10;
X = ( float* ) malloc ( sizeof ( float ) * N );
Y = ( float* ) malloc ( sizeof ( float ) * N );
for ( i = 0; i < N; i++ ) {
X[i] = 1.0;
Y[i] = 2.0;
}
resultFloat = cblas_sdot ( N, X, 1, Y, 1 );
printf("the result of cblas_sdot is %4.4f\n", resultFloat);
//deallocate space
free ( X );
free ( Y );
}
int main(int argc, const char * argv[]) {
runSDOTSample();
return 0;
}
Práctica 4¶
Imagine que desea sumar dos vectores llamados a y b. El vector a contiene los cuadrados de los enteros 0 hasta n (excluyendo n); por ejemplo, si n = 3, entonces a = (0, 1, 4). El vector b contiene los cubos de los enteros 0 hasta n (excluyendo n); por ejemplo, si n = 3, entonces b = (0, 1, 8). ¿Cómo haría eso con Python únicamente?
Sumando vectores utilizando Python “puro”¶
La siguiente función resuelve el problema de la suma de vectores utilizando Python sin NumPy:
def pythonsum(n):
a = range(n)
b = range(n)
c = []
for i in range(len(a)):
a[i] = i ** 2
b[i] = i ** 3
c.append(a[i] + b[i])
return c
Sumando vectores utilizando NumPy¶
La siguiente función resuelve el mismo problema con NumPy:
def numpysum(n):
a = numpy.arange(n) ** 2
b = numpy.arange(n) ** 3
c = a + b
return c
Comparando los tiempos de ejecución¶
Con ayuda del siguiente código, compare los tiempos de ejecución para n = 1000, 2000, 3000.
- Utilice Emacs para editar su código.
- Ejecute su programa desde la línea de comandos.
- Utilice Emacs para ejecutar su programa (pida instrucciones a su profesor).
import sys
from datetime import datetime
import numpy
def numpysum(n):
a = numpy.arange(n) ** 2
b = numpy.arange(n) ** 3
c = a + b
return c
def pythonsum(n):
a = range(n)
b = range(n)
c = []
for i in range(len(a)):
a[i] = i ** 2
b[i] = i ** 3
c.append(a[i] + b[i])
return c
size = int(sys.argv[1])
start = datetime.now()
c = pythonsum(size)
delta = datetime.now() - start
print "The last 2 elements of the sum", c[-2:]
print "PythonSum elapsed time in microseconds", delta.microseconds
start = datetime.now()
c = numpysum(size)
delta = datetime.now() - start
print "The last 2 elements of the sum", c[-2:]
print "NumPySum elapsed time in microseconds", delta.microseconds
Arreglos multidimensionales¶
Los arreglos de NumPy son homogéneos, es decir, los elementos del arreglo tienen que ser del mismo tipo. La ventaja es que ahora es fácil determinar el tamaño que se requiere para almacenar el arreglo. Los arreglos de NumPy se indexan como en Python, a partir de 0.
Utilizando el siguiente código, averigüe cuál es el tipo de los elementos del arreglo a (las respuestas son diferentes dependiendo de si utiliza Python de 32 bits o de 64 bits):
>>> a = arange(5)
>>> a.dtype
Además del tipo de los elementos del arreglo, también es importante conocer su “forma”:
>>> a
>>> array([0, 1, 2, 3, 4])
>>> a.shape
>>> (5,)
El atributo shape (forma) del arreglo es una tupla, en este caso, una tupla de 1 elemento, que contiene el tamaño del arreglo en cada dimensión.
Creando un arreglo multidimensional¶
Muestre la forma y el tipo de los elementos del arreglos siguientes:
>>> m1 = array([arange(2), arange(2)])
>>> m2 = array([ [1, 2], [3, 4] ])
Ahora cree una matriz de tamaño \(3 \times 3\) con los siguientes elementos
Compruebe que también puede crearlo de la siguiente manera:
>>> m3 = arange(9).reshape(3, 3)
Los elementos de m3 pueden accederse mediante la siguiente sintáxis:
>>> m3[0, 0]
¿Cuál es la salida de m3[2, 1]?
Un poco de Álgebra Lineal con NumPy¶
El paquete numpy.linalg contiene funciones de álgebra lineal. Con este módulo, usted puede invertir matrices, calcular eigen-valores, resolver sistemas de ecuaciones lineales, calcular determinantes, entre otras cosas.
La inversa de una matriz¶
La inversa de una matriz cuadrada \(A\) de \(n \times n\) se denota por \(A^{-1}\), que cuando se multiplica por la matriz original, es igual a la matriz identidad. Esto puede escribirse como se indica a continuación:
La función inv del paquete numpy.linalg puede invertir una matriz por nosotros. Vamos a resolver un ejemplo:
Crea la matriz \(A\) con la función mat:
>>> A = numpy.mat("0 1 2; 1 0 3; 4 -3 8")
Invierte la matriz \(A\) con la función inv:
>>> A_inv = numpy.linalg.inv(A)
Si la matriz es singular, o no es cuadrada, resulta un error LinAlgError.
Verifica utilizando multiplicación matricial:
>>> print "Verifica\n", A * A_inv
Resolviendo un sistema de ecuaciones lineales¶
La función solve de numpy.linalg resuelve sistemas de ecuaciones de la forma \(Ax = b\), donde \(A\) es una matriz, \(b\) puede ser un arreglo unidimensional o bidimensional, y \(x\) es la variable desconocida.
Crea la matriz \(A\) y el vector \(b\):
>>> A = numpy.mat("1 -2 1; 0 2 -8; -4 5 9") >>> print "A\n", A >>> b = numpy.array([0, -8, 9]) >>> print "b\n", b
Invoque la función solve para resolver el sistema de ecuaciones lineales:
>>> x = numpy.linalg.solve(A, b) >>> print "Solucion: ", x
Verifique que la solución es correcta (dot realiza una multiplicación matricial):
>>> print "Verifica\n", numpy.dot(A, x)
Gráficas simples con Matplotlib¶
La visualización de una curva \(f(x)\) se realiza al graficar lineas rectas entre \(n\) puntos en la curva. Entre más puntos se usen, más suave es la forma de la gráfica. Com ejemplo, vamos a graficar la curva \(y = t^2 e^{-t^2}\) para valores de \(t\) entre 0 y 3. Primero generamos coordenadas igualmente espaciadas para \(t\), para este ejemplo, 51 valores (50 intervalos). Entonces calculamos los correspondientes valores para \(y\) en estos puntos, antes de invocar el comando plot(t, y) para obtener la gráfica. A continuación mostramos el programa:
from matplotlib.pylab import *
def f(t):
return t**2*exp(-t**2)
t = linspace(0, 3, 51) # 51 puntos entre 0 y 3
y = zeros(len(t)) # reservar espacio de memoria para y con elementos flotantes.
# Se puede hacer utilizando ciclos, pero es mas lento
# for i in xrange(len(t)):
# y[i] = f(t[i])
# Se obtiene una ejecucion mas rapida (y la sintaxis tambien es mas simple)
# actuando sobre todo el arreglo t
y = f(t)
plot(t, y)
xlabel('t')
ylabel('y')
legend(['t^2*exp(t^2)'])
show()
(Source code, png, hires.png, pdf)
