en construcción, revisado el 17/05/17
Alberto Ruiz, DIS-UM
Python es un lenguaje de propósito general que permite acceder cómodamente a una amplísima colección de bibliotecas útiles en todos los campos de la informática. Su uso en ciencia y tecnología es cada vez mayor.
Es un lenguaje interpretado, con tipos de datos dinámicos y manejo automático de memoria, que puede utilizarse tanto para escribir programas de forma tradicional como para experimentar en un entorno interactivo. Incluye las construcciones más importantes de programación funcional y también admite programación orientada a objetos.
Aunque su sintaxis es sencilla e intuitiva tiene algunas diferencias con respecto a otros lenguajes usados en el ámbito científico que inicialmente pueden inducir a confusión:
Los bloques de instrucciones de las construcciones condicionales, bucles y funciones se delimitan mediante la "indentación del código": no se utiliza "end" ni {
}
.
Los índices para acceder a los arrays o listas comienzan en 0 y acaban en tamaño-1. Las secuencias (range) usadas en bucles o list comprehensions no incluyen el límite superior.
Algunas funciones tienen la sintaxis tradicional f(x)
, g(x,a)
, mientras que otras se expresan como x.f()
, x.g(a)
, etc., indicando que el "objeto" x
se modifica de alguna forma.
Los arrays y las listas son "mutables": su asignación a otra variable no crea una copia del objeto original sino una "referencia" a través de la cual se puede modificar la estructura original.
Las funciones pueden leer directamente el valor de variables globales, pero para modificarlas hay que declararlas como global
. La asignación de variables dentro de una función crea variables locales.
Cadenas de caracteres:
s = 'Hola'
s
'Hola'
print(s)
Hola
type(s)
str
Se admiten diferentes tipos de delimitadores y cadenas multilínea.
type("Hola ")
str
type('''Hola''')
str
Variables lógicas:
c = 3 < 4
type(c)
bool
c and (2==1+1) or not (3 != 5)
True
Números reales aproximados con coma flotante de doble precisión:
x = 3.5
type(x)
float
Los enteros tienen tamaño ilimitado:
x = 20
type(x)
int
x**x
104857600000000000000000000
Variable compleja:
(1+1j)*(1-1j)
(2+0j)
import cmath
cmath.sqrt(-1)
1j
Se usan normalmente para almacenar unos pocos datos, posiblemente de diferentes tipos:
t = (2,'rojo')
t
(2, 'rojo')
t[0]
2
Son inmutables.
l = [1,2,3,4,5]
type(l)
list
l[0]
1
l[-1]
5
len(l)
5
l[1:3]
[2, 3]
[1,2,3] + [4,5]
[1, 2, 3, 4, 5]
Lo normal es guardar datos del mismo tipo aunque no es obligatorio.
s =[1,2,3,"hola",True]
type(s)
list
Las listas y arrays son referencias a objetos mutables:
x = 2
y = x
y = 100
x
2
x = [1,2]
y = x
y[0] = 100
x
[100, 2]
Para evitar este comportamiento podemos hacer una copia:
x = [1,2]
y = x.copy()
y[0] = 100
x
[1, 2]
import numpy as np
a = np.array([1,2,3])
b = a
c = a.copy()
b[1] = 100
a, c
(array([ 1, 100, 3]), array([1, 2, 3]))
Relacionado con esto, algunas operaciones modifican la estructura original y otras no:
m = np.array(range(10))
m.resize([2,5])
m
array([[0, 1, 2, 3, 4], [5, 6, 7, 8, 9]])
m = np.array(range(10))
np.resize(m,[2,5])
array([[0, 1, 2, 3, 4], [5, 6, 7, 8, 9]])
m
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
m = np.array(range(10))
m.reshape(2,5)
array([[0, 1, 2, 3, 4], [5, 6, 7, 8, 9]])
m
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
Es un array asociativo (el índice puede ser cualquier tipo (inmutable)). Es una estructura muy utilizada en Python.
d = {'lunes': 8, 'martes' : [1,2,3], 3: 5}
d['martes']
[1, 2, 3]
d.keys()
dict_keys([3, 'lunes', 'martes'])
d.values()
dict_values([5, 8, [1, 2, 3]])
Condiciones:
k = 7
if k%2 == 0:
print(k," es par")
else:
print(k," es impar")
print("me gustan los impares")
7 es impar me gustan los impares
En Python 3 las secuencias entregan los elementos "bajo demanda", no crean una lista completa.
range(3,8)
range(3, 8)
list(range(3,8))
[3, 4, 5, 6, 7]
Bucles:
for k in [1,2,3]:
print(k)
1 2 3
for k in range(5):
print(k)
0 1 2 3 4
k = 1
p = 1
while k < 5:
p = p*k
k = k+1
p
24
Hay muchos tipos de contenedores que pueden recorrerse con bucles: listas, arrays, diccionarios, etc.
def sp(n):
r = n**2+n+41
return r
sp(5)
71
Se pueden devolver varios resultados en una tupla:
import math
def ecsec(a,b,c):
d = math.sqrt(b**2- 4*a*c)
s1 = (-b+d)/2/a
s2 = (-b-d)/2/a
return (s1,s2)
ecsec(2,-6,4)
(2.0, 1.0)
Los paréntesis de la tupla son opcionales.
a,b = ecsec(1,-3,2)
b
1.0
Las variables globales son visibles dentro de las funciones y las asignaciones crean variables locales (a menos que el nombre se declare global
).
a = 5
b = 8
def f(x):
b = a+1
return b
print(f(3))
print(b)
6 8
Argumentos por omisión:
def incre(x,y=1):
return x + y
print(incre(5))
print(incre(5,3))
6 8
Argumentos por nombre:
incre(y=3, x=2)
5
Documentación:
# ? sum
help(sum)
Help on built-in function sum in module builtins: sum(iterable, start=0, /) Return the sum of a 'start' value (default: 0) plus an iterable of numbers When the iterable is empty, return the start value. This function is intended specifically for use with numeric values and may reject non-numeric types.
def fun(n):
"""Una función muy simple que calcula el triple de su argumento."""
return 3*n
help(fun)
Help on function fun in module __main__: fun(n) Una función muy simple que calcula el triple de su argumento.
Las funciones definidas en un archivo se pueden utilizar directamente haciendo un import
. Existe una convención para definir una función main
que se ejecuta cuando el archivo se arranca como programa y suele usarse para ejecutar tests.
En Python 3 las construcciones funcionales crean secuencias "bajo demanda".
map(sp,range(5))
<map at 0x7f768c0d2128>
for k in map(sp,range(5)):
print(k)
41 43 47 53 61
list(map(sp,range(5)))
[41, 43, 47, 53, 61]
list(filter(lambda x: x%2 == 1, range(10)))
[1, 3, 5, 7, 9]
List comprehensions:
[k**2 for k in range(10) if k >5 ]
[36, 49, 64, 81]
def divis(n):
return [k for k in range(2,n) if n%k==0]
divis(12)
[2, 3, 4, 6]
divis(1001)
[7, 11, 13, 77, 91, 143]
def perfect(n):
return sum(divis(n)) + 1 == n
perfect(4)
False
perfect(6)
True
def prime(n):
return divis(n)==[]
[k for k in range(2,21) if prime(k)]
[2, 3, 5, 7, 11, 13, 17, 19]
Función que construye funciones:
def mkfun(y):
return lambda x: x+y
f = mkfun(1)
g = mkfun(5)
print(f(10))
print(g(10))
11 15
fs = list(map(mkfun,range(1,6)))
print(fs[0](10))
print(fs[4](10))
11 15
Existe reduce
(fold) pero parece que no es muy "idiomático".
Hay muchos paquetes para producir gráficos espectaculares. Uno de los más conocidos es matplotlib
, que puede utilizarse con un interfaz muy parecido al de Matlab/Octave.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
x=np.linspace(0,2*np.pi,200)
plt.plot(np.sin(x))
[<matplotlib.lines.Line2D at 0x7f7678b9f2e8>]
plt.plot(np.cos(x),np.sin(x)); plt.axis('equal')
(-1.0998691589935401, 1.0999937694758828, -1.099965731583572, 1.099965731583572)
plt.plot(x,np.sin(x), x,np.cos(x));
plt.plot(x,np.sin(x),'r',x,np.sin(2*x),'k',[1,2.5],[-0.5,0],'bo');
plt.legend(['hola','bla','bla bla']);
plt.xlabel('x'); plt.ylabel('y'); plt.title('bonito plot'); plt.axis('tight');
plt.plot(x,np.exp(x)); plt.axis([0,3,-1,5]);
for k in [1,2,3]:
plt.plot(x,np.sin(k*x))
plt.grid()
def espiral(n):
t = np.linspace(0,n*2*np.pi,1000)
r = 3 * t
x = r * np.cos(t)
y = r * np.sin(t)
plt.plot(x,y)
plt.axis('equal')
plt.axis('off')
espiral(4)
import numpy.random as rnd
def randwalk(n,s):
p = s*rnd.randn(n,2)
r = np.cumsum(p,axis=0)
x = r[:,0]
y = r[:,1]
plt.plot(x,y)
plt.axis('equal');
import matplotlib.pylab as pylab
print(pylab.rcParams['figure.figsize'])
pylab.rcParams['figure.figsize'] = 5, 5
[6.0, 4.0]
randwalk(1000,1)
pylab.rcParams['figure.figsize'] = 8, 8
x = np.linspace(0,6*np.pi,100);
plt.subplot(2,2,1)
plt.plot(x,np.sin(x),'r')
plt.subplot(2,2,2)
plt.plot(x,np.cos(x))
plt.subplot(2,2,3)
plt.plot(x,np.sin(2*x))
plt.subplot(2,2,4)
plt.plot(x,np.cos(2*x),'g');
pylab.rcParams['figure.figsize'] = 6, 6
x,y = np.mgrid[-3:3:0.2,-3:3:0.2]
plt.contour(x,y, x**2-y**2-1, 0, colors=['k']);
plt.axis('equal');
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x,y,x**2-y**2, cmap=cm.coolwarm, linewidth=0.5, rstride=2, cstride=2);
pylab.rcParams['figure.figsize'] = 6, 4
import numpy.random as rnd
rnd.randn(2,3)
array([[ 1.10109629, -0.72138553, -1.2202826 ], [-0.53595346, -0.44037795, -0.98122434]])
r = 2+3*rnd.randn(100)
plt.hist(r,10);
np.mean(r)
2.2170177852015147
np.std(r)
2.938161931731035
np.median(r)
2.44634684308995
plt.plot(np.sort(r),range(1,101));
np.matrix('1 2; 3 5')
matrix([[1, 2], [3, 5]])
m = np.matrix([[1, 2, 3, 4]
,[5, 6, 7, 8]
,[9,10,11,12]])
m
matrix([[ 1, 2, 3, 4], [ 5, 6, 7, 8], [ 9, 10, 11, 12]])
m.transpose()
matrix([[ 1, 5, 9], [ 2, 6, 10], [ 3, 7, 11], [ 4, 8, 12]])
m.T
matrix([[ 1, 5, 9], [ 2, 6, 10], [ 3, 7, 11], [ 4, 8, 12]])
En los arrays de un índice (vectores) no se distinguen filas y columnas:
np.array([1,2,3]).T
array([1, 2, 3])
pero en los de dos índices (matrices), sí:
np.array([[1,2,3]]).T
array([[1], [2], [3]])
Matrices por bloques
np.bmat([[m, np.zeros([3,3])]
,[np.eye(3), np.ones([3,4]) ]])
matrix([[ 1., 2., 3., 4., 0., 0., 0.], [ 5., 6., 7., 8., 0., 0., 0.], [ 9., 10., 11., 12., 0., 0., 0.], [ 1., 0., 0., 1., 1., 1., 1.], [ 0., 1., 0., 1., 1., 1., 1.], [ 0., 0., 1., 1., 1., 1., 1.]])
"Automatic broadcasting"
m + np.array([[10],[20],[30]])
matrix([[11, 12, 13, 14], [25, 26, 27, 28], [39, 40, 41, 42]])
Extracción de elementos y submatrices
m
matrix([[ 1, 2, 3, 4], [ 5, 6, 7, 8], [ 9, 10, 11, 12]])
m[1,2]
7
m[2,1:4]
matrix([[10, 11, 12]])
m[:2, 2:]
matrix([[3, 4], [7, 8]])
m[[1,0,0,2,1],:]
matrix([[ 5, 6, 7, 8], [ 1, 2, 3, 4], [ 1, 2, 3, 4], [ 9, 10, 11, 12], [ 5, 6, 7, 8]])
Extracción de elementos que cumplen una condición:
n = np.array(range(1,11))
n
array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
n < 5
array([ True, True, True, True, False, False, False, False, False, False], dtype=bool)
n[n<5]
array([1, 2, 3, 4])
k = np.linspace(1,100,100)
(k ** 2)[np.logical_and(k>10 , k**3 < 2000)]
array([ 121., 144.])
Más ejemplos:
s = np.array([1,2,3,4])
type(s)
numpy.ndarray
np.ndim(s)
1
np.shape(s)
(4,)
sum(s)
10
x = np.array([k for k in range(1,11)])
print(x)
[ 1 2 3 4 5 6 7 8 9 10]
np.append(x,[13,14])
array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 13, 14])
np.size(x)
10
for x in np.linspace(0,1,10):
print(x)
0.0 0.111111111111 0.222222222222 0.333333333333 0.444444444444 0.555555555556 0.666666666667 0.777777777778 0.888888888889 1.0
np.sqrt(4)
2.0
import numpy as np
import scipy.linalg as la
Resuelve
$$ \begin{align*} x + 2y &= 3\\ 3x+4y &= 5 \end{align*} $$m = np.matrix([[1,2],[3,4]])
m
matrix([[1, 2], [3, 4]])
np.shape(m)
(2, 2)
la.inv(m)
array([[-2. , 1. ], [ 1.5, -0.5]])
x = la.solve(m,[[3],[5]])
x
array([[-1.], [ 2.]])
np.dot(m,x)
matrix([[ 3.], [ 5.]])
m @ x
matrix([[ 3.], [ 5.]])
np.shape(x)
(2, 1)
np.ndim(x)
2
Como ejemplo de uso de las herramientas de álgebra lineal realizaremos el ajuste de un modelo polinomial a unas observaciones ficticias. Encontraremos la solución de mínimo error cuadrático a un sistema de ecuaciones sobredeterminado.
En primer lugar generamos unos datos de prueba artificiales que simulan observaciones contaminadas con ruido de una función no lineal.
x = np.linspace(0,2,30)
y = np.sin(x) + 0.05*rnd.randn(x.size)
plt.plot(x,y,'.');
Vamos a ajustar un modelo del tipo $y = ax^2 + bx + c$. Los coeficientes desconocidos $a$, $b$ y $c$ se pueden obtener resolviendo un sistema de ecuaciones lineales.
La matriz de coeficientes tiene potencias de $x$ hasta el grado que nos interesa.
A = np.vstack([x**2, x, np.ones(x.size)]).transpose()
A
array([[ 0. , 0. , 1. ], [ 0.00475624, 0.06896552, 1. ], [ 0.01902497, 0.13793103, 1. ], [ 0.04280618, 0.20689655, 1. ], [ 0.07609988, 0.27586207, 1. ], [ 0.11890606, 0.34482759, 1. ], [ 0.17122473, 0.4137931 , 1. ], [ 0.23305589, 0.48275862, 1. ], [ 0.30439952, 0.55172414, 1. ], [ 0.38525565, 0.62068966, 1. ], [ 0.47562426, 0.68965517, 1. ], [ 0.57550535, 0.75862069, 1. ], [ 0.68489893, 0.82758621, 1. ], [ 0.80380499, 0.89655172, 1. ], [ 0.93222354, 0.96551724, 1. ], [ 1.07015458, 1.03448276, 1. ], [ 1.2175981 , 1.10344828, 1. ], [ 1.3745541 , 1.17241379, 1. ], [ 1.54102259, 1.24137931, 1. ], [ 1.71700357, 1.31034483, 1. ], [ 1.90249703, 1.37931034, 1. ], [ 2.09750297, 1.44827586, 1. ], [ 2.3020214 , 1.51724138, 1. ], [ 2.51605232, 1.5862069 , 1. ], [ 2.73959572, 1.65517241, 1. ], [ 2.97265161, 1.72413793, 1. ], [ 3.21521998, 1.79310345, 1. ], [ 3.46730083, 1.86206897, 1. ], [ 3.72889417, 1.93103448, 1. ], [ 4. , 2. , 1. ]])
En realidad es una matriz de Vandermonde:
np.vander(x,3)
array([[ 0. , 0. , 1. ], [ 0.00475624, 0.06896552, 1. ], [ 0.01902497, 0.13793103, 1. ], [ 0.04280618, 0.20689655, 1. ], [ 0.07609988, 0.27586207, 1. ], [ 0.11890606, 0.34482759, 1. ], [ 0.17122473, 0.4137931 , 1. ], [ 0.23305589, 0.48275862, 1. ], [ 0.30439952, 0.55172414, 1. ], [ 0.38525565, 0.62068966, 1. ], [ 0.47562426, 0.68965517, 1. ], [ 0.57550535, 0.75862069, 1. ], [ 0.68489893, 0.82758621, 1. ], [ 0.80380499, 0.89655172, 1. ], [ 0.93222354, 0.96551724, 1. ], [ 1.07015458, 1.03448276, 1. ], [ 1.2175981 , 1.10344828, 1. ], [ 1.3745541 , 1.17241379, 1. ], [ 1.54102259, 1.24137931, 1. ], [ 1.71700357, 1.31034483, 1. ], [ 1.90249703, 1.37931034, 1. ], [ 2.09750297, 1.44827586, 1. ], [ 2.3020214 , 1.51724138, 1. ], [ 2.51605232, 1.5862069 , 1. ], [ 2.73959572, 1.65517241, 1. ], [ 2.97265161, 1.72413793, 1. ], [ 3.21521998, 1.79310345, 1. ], [ 3.46730083, 1.86206897, 1. ], [ 3.72889417, 1.93103448, 1. ], [ 4. , 2. , 1. ]])
El lado derecho del sistema es directamente el vector con los valores de $y$, la variable independiente del modelo.
B = np.matrix(y).T
B
matrix([[-0.03017562], [ 0.12102307], [ 0.20472809], [ 0.22336368], [ 0.29015562], [ 0.28786901], [ 0.39513457], [ 0.46121004], [ 0.49673906], [ 0.55230101], [ 0.65127628], [ 0.65792914], [ 0.77057348], [ 0.79241806], [ 0.76122737], [ 0.94307161], [ 0.8812106 ], [ 0.89471824], [ 0.96952483], [ 0.95853218], [ 0.99921169], [ 0.96480888], [ 1.00485747], [ 0.94585026], [ 1.07448525], [ 0.94431881], [ 0.96184194], [ 0.96520362], [ 0.96336003], [ 0.93290306]])
El sistema que hay que resolver está sobredeterminado: tiene solo tres incógnitas y tantas ecuaciones como observaciones de la función.
$$A \begin{bmatrix}a\\b\\c\end{bmatrix}= B$$La solución de mínimo error cuadrático para los coeficientes del modelo se obtiene de manera directa con el operador de división matricial.
sol = la.lstsq(A,B)[0]
sol
array([[-0.36807329], [ 1.21964859], [-0.01910038]])
#ye = np.dot(A,sol)
ye = A @ sol
plt.plot(x,y,'.',x,ye,'r');
Se puede experimentar con polinomios de mayor o menor grado.
Otra forma de llegar a esta solución es resolver las "ecuaciones normales" del sistema $Ax=b$. Como está sobredeterminado, vamos a minimizar el error cuadrático $E=||Ax-b||^2$, que es matemáticamente tratable. Hay que resolver
$$\frac{\partial E}{\partial x_k} =0 $$Se puede comprobar que la solución es
$$x = A^+ b$$donde $A^+$ es la "pseudoinversa" de $A$:
$$A^+ = (A^TA)^{-1}A^T$$la.inv(A.T @ A) @ A.T @ B
matrix([[-0.36807329], [ 1.21964859], [-0.01910038]])
la.pinv(A) @ B
matrix([[-0.36807329], [ 1.21964859], [-0.01910038]])
Resuelve
$$x^4=16$$import scipy as sci
sci.roots([1,0,0,0,-16])
array([ -2.00000000e+00+0.j, 1.66533454e-16+2.j, 1.66533454e-16-2.j, 2.00000000e+00+0.j])
Resuelve
$$sin(x)+cos(2x)=0$$import scipy.optimize as opt
opt.fsolve(lambda x: sci.sin(x) + sci.cos(2*x), 0)
array([-0.52359878])
Resuelve
$$ \begin{align*} x^2 - 3y &= 10\\ sin(x)+y &= 5 \end{align*} $$def fun(z):
x,y = z
return [ x**2 - 3*y - 10
, sci.sin(x) + y - 5]
opt.fsolve(fun,[0.1,-0.1])
array([ 5.2511881 , 5.85832548])
Encuentra $(x,y)$ que minimiza $(x-1)^2 + (y-2)^2-x+3y$
def fun(z):
x,y = z
return (x-1)**2 + (y-2)**2 - x + 3*y
opt.minimize(fun,[0.1,-0.1])
fun: 2.500000000000014 hess_inv: array([[ 0.57758622, -0.18103452], [-0.18103452, 0.92241375]]) jac: array([ 0.00000000e+00, -2.38418579e-07]) message: 'Optimization terminated successfully.' nfev: 12 nit: 2 njev: 3 status: 0 success: True x: array([ 1.49999999, 0.49999988])
r = np.linspace(0,2,30)
x,y = np.meshgrid(r,r)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
z = fun((x,y))
ax.plot_surface(x,y,z, cmap=cm.coolwarm, linewidth=0.5, rstride=2, cstride=2);
Calcula una aproximación numérica para $f'(2)$ cuando $f(x) = \sin(2x)*\exp(\cos(x))$
from scipy.misc import derivative
derivative(lambda x: sci.sin(2*x)*sci.exp(sci.cos(x)),2,1E-6)
-0.40836700757052036
(lambda x: (-np.sin(x)*np.sin(2*x) + 2*np.cos(2*x))*np.exp(np.cos(x)))(2)
-0.40836700756782335
Calcula una aproximación numérica a la integral definida
$$\int_0^1 \frac{4}{1+x^2}dx$$from scipy.integrate import quad
quad(lambda x: 4/(1+x**2),0,1)
(3.1415926535897936, 3.4878684980086326e-14)
Resuelve
$$\ddot{x}+0.95x+0.1\dot{x}=0$$para $x(0)=10$, $\dot{x}(0)=0, t\in[0,20]$
from scipy.integrate import odeint
def xdot(z,t):
x,v = z
return [v,-0.95*x-0.1*v]
t = np.linspace(0,20,1000)
r = odeint(xdot,[10,0],t)
# plt.plot(r);
plt.plot(t,r[:,0],t,r[:,1]);
plt.plot(r[:,0],r[:,1]);
Podemos "tunearlo":
from IPython.display import HTML
HTML(open('../nb1.css').read())