JSM_Programacion

1550 days ago by trinidad

Programación en Sage

El lenguaje de programación de Sage por defecto es Python. Según Wikipedia:

"Python es un lenguaje de programación interpretado cuya filosofía hace hincapié en una sintaxis muy limpia y que favorezca un código legible.

Se trata de un lenguaje de programación multiparadigma, ya que soporta orientación a objetos, programación imperativa y, en menor medida, programación funcional. Es un lenguaje interpretado, usa tipado dinámico y es multiplataforma.

Es administrado por la Python Software Foundation. Posee una licencia de código abierto, denominada Python Software Foundation License,1 que es compatible con la Licencia pública general de GNU a partir de la versión 2.1.1, e incompatible en ciertas versiones anteriores."

Variables, variables, funciones y funciones

Como en Sage usaremos Python en un entorno matemático, tendremos conceptos distintos con el mismo nombre. Para evitar confusiones, vamos a comenzar definiendo esos conceptos y viendo algunos ejemplos.

Variables

Una variable en Python será un identificador en el que se guarda un "objeto". A diferencia de otros lenguajes de programación, no es necesario decir a priori qué es lo que vamos a guardar en la variable.

a=5;a 
       

Además, Sage maneja las variables de modo dinámico, por lo que, si queremos podemos guardar otra cosa en la misma variable:

a='Hola';a 
       

Una vez tenemos guardado algo, podemos utilizarlo de la misma forma que utilizaríamos el objeto original:

a=a+" y adios";a 
       

Recordemos que para saber qué operaciones podemos hacer con la variable, basta escribir el nombre, seguido de un punto y pulsar <tab>.

a.count("a") 
       

Si queremos borrar una variable, podemos hacerlo con del:

del a;a 
       

También podemos asignar varias variables al mismo tiempo:

a,b=1,2; a,b 
       
a,b=b,a+b; a,b 
       

Variables (en sentido matemático)

En Sage también tenemos "variables" en sentido matemático. Por defecto la única que cumple esto es la letra "x":

(x+3)^3+x+1 
       

Si necesito más variables, puedo declararlas con var:

var("y") 
       
x+y-x^2 
       

Hemos visto antes que en una variable se pueden guardar valores. También se pueden guardar expresiones matemáticas:

p=x^3+x^2+1;p 
       
       

¿Qué ocurre si asignamos un valor a la variable?

x=3 
       
p+x 
       

Podemos recuperar la incógnita usando restore

restore('x');x 
       

Funciones (matemáticas)

Otro concepto matemático que podemos usar en Sage es el de función.

f(x)=x^3+3;f 
       

Cuando definimos una función, Sage entiende que las variables que ponemos son simbólicas:

f(x,y,z)=x*y*z+1;f 
       

Y las funciones admiten las operaciones normales para funciones matemáticas:

f.diff(x) 
       

Funciones

Para definir una función (programación) en Sage/Python, usamos la palabra clave def, seguida del nombre de la función, las variables de entrada entre paréntesis y dos puntos. En las siguientes líneas se escriben los comandos que queremos que evalúe la función.La identación es muy importante en Python: todo lo que esté con la misma identación será el cuerpo de la función, para terminar la función basta escribir sin identación. Por último, el valor que devuelve la función se determina por la palabra clave return.

def g(x,y): x=x+1; y=y+2; return x*y 
       

Podemos usar la función igual que usamos cualquier función matemática.

g(1,3) 
       

Listas, bucles y y bucles en listas

Un tipo de datos muy importante en Sage son las listas. Vamos a ver cómo crear, modificar y seleccionar elementos de ellas.

Una lista no es más que una colección ordenada de elementos. Cada elemento de la lista puede ser cualquier tipo de objeto. Para crear una lista se pone entre corchetes la lista de los elementos, separados por comas.

L = [2,3,5,7,11,13,17,2] 
       

Para acceder a un elemento se pone el nombre de la lista, seguido de corchetes y entre corchetes el número de elemento al que queremos acceder. El primer elemento se numera con 0.

L[0] 
       
L[0]=3;L 
       

También podemos comenzar contando desde el último elmento usando índices negativos. Así el último elemento está numerado con -1, el penúltimo con -2, etc

L[-2] 
       

Podemos seleccionar un segmento de la lista con el operador :, siguiendo la convención nombreDeLista[primerElemento:ultimoElemento].

L[2:5] 
       

Si el último elemento que queremos elegir es el último elemento de la lista, no es necesario ponerlo.

L[2:] 
       
L[:-2] 
       

Recordemos que los elementos de una lista pueden ser de distinto tipo, aunque no es usual emplearlas de este modo.

Lista_rara=[1,"Hola",[1,2,3]];Lista_rara 
       

Podemos añadir un elemento al final de una lista con append, añadir una lista con extend o símplemente unir varias listas con la operación +.

L.append(4) 
       
L.extend([10,11,12]) 
       
[1,3,5]+[2,4,6]+[100] 
       

Creación de listas

Python permite crear de modo sencillo una lista de elementos equiespaciados. En Sage esta utilidad está extendida también para trabajar de modo simbólico.

[1..7] 
       
[2,4..10] 
       
[pi,4*pi..32] 
       
show([x^2,3*x^2/2..2*pi*x^2]) 
       

Operando con listas

Sage soporta muchas operaciones directamente sobre listas. Como ejemplo, vamos a sumar los elementos de una lista, vamos a agrupar dos listas en una uniendo el elemento k de la primera con el k de la segunda, y por último map que aplica una función a cada elemento de una lista.

sum([1..100]) 
       
zip([1,2,3,4],['a','b','c','d'] ) 
       
map( cos, [0, pi/4, pi/2, 3*pi/4, pi] ) 
       

Estructuras de control

Las estructuras de control seleccionan cuál es la siguiente sentencia a ejecutar (selección -- if) o hacen que una sentencia se ejecute repetidamente (iteración -- while o for).

 

En Python, la selección tiene la siguiente sintaxis:

n=12343; if n%2 == 0: print n/2 else: print (n+1)/2 
       

Notese que la indentación es esencial: todas las sentencias que queremos que se ejecuten han de tener la misma identación.

Para repetir una sentencia un número indeterminado de veces, podemos usar un bucle while, con sintaxis:

i=0; while i < 5: i=i+1; print i^2 
       

Finalmente, el bucle for repite una instrucción para cada elemento de una lista:

for i in [0..4]: print i^2 
       
for str in ["juan","mario","eva","laura"]: print str.capitalize() 
       

Listas por comprensión

Python permite definir listas por comprensión. Es una manera rápida de crear una lista a partir de otra. La sintasis es [ expresión for indice in lista]. Devuelve una lista resultado de evaluar la expresión (que dependerá de la variable denominada índice) aplicada a cada uno de los elementos de la lista.

[ 2*k for k in [0..10] ] 
       
[ cos(x) for x in [pi/4, pi/2..2*pi]] 
       

Además se puede filtrar el resultado de la salida:

U = [ k for k in [1..19] if gcd(k,20) == 1]; U 
       

Y se pueden anidar bucles:

[ (a,b) for a in U for b in U ] 
       

Por último, la expresión inicial también puede contener condicionales:

[ 'prime' if x.is_prime() else 'not prime' for x in U] 
       

Otros tipos similares son secuencias y tuplas. En las tuplas no se pueden eliminar ni añadir elementos (son equivalentes al "array" en C). La secuencia es similar a la lista, pero puede elegirse que sea inmutable.


Sage es orientado a objetos

Anteriormente hemos definido el concepto de variables. Sin embargo, en Python las variables son algo más, pues conocen qué se puede hacer con ellas. Esto se denomina objetos.

La primera diferencia frente a una variable normal nos la encontramos cuando igualamos dos objetos. Python no copia el contenido de una variable en la otra, sino que asigna la nueva etiqueta al objeto que tenía la etiqueta. Veámoslo con un ejemplo:

a = [1,2,3] b = a # Asignamos el nombre b al objeto guardado en a. Ahora el objeto tiene dos nombres. b[1] = 0 # Modificamos uno de ellos a # y el otro se modifica. 
       

Un objeto en Python tiene asignados diversos métodos, es decir, operaciones que se pueden hacer con el objeto. Para ver las operaciones podemos escribir el objeto, seguido de un punto y pulsar <tab>:

a.append(3);a 
       

Clases

Vamos a construir una pequeña clase en Python.

class Glass(object): def __init__(self, size): # Cuando se crea el vaso se llama a este método self._size = size # Tamaño del vaso self._content = "nada" # Contenido def __repr__(self): # Este método controla cómo se muestra por pantalla if self._content == "nada": return "Un vaso vacío de tamaño %s"%(self._size) else: return "Un vaso de tamaño %s lleno de %s"%(self._size,self._content) def fill(self,content): # Método para llenar el vaso self._content = content def empty(self): # Método para vaciar el vaso self._content = "nada" 
       
myGlass = Glass("grande"); myGlass 
       
myGlass.fill("cerveza"); myGlass 
       
myGlass.empty(); myGlass 
       

Una convención en Python es que los métodos o atributos que comienzan con _ no deben ser utilizados, sino que se usan internamente.

En la primera línea hemos puesto entre los paréntesis "object". Toda clase "hereda" de otra. Cuando definimos una clase genérica, heredamos de la clase "object" (en Sage normalmente de SageObject). Pero podemos heredar de cualquier otra clase. De ese modo, tenemos todos los métodos y propiedades definidos para dicha clase.

class WineGlass(Glass): def fill(self): # Método para llenar el vaso self._content = "vino" 
       
vasoVino=WineGlass("mediano");vasoVino 
       
vasoVino.fill();vasoVino 
       

Elementos, clases y categorías

Todo objeto en Sage pertenece a una clase. Podemos conocer cuál es esa clase mediante la función type o el método parent().

Por ejemplo, el número 1 es un objeto de tipo número entero (un elemento del anillo de los números enteros).

type(1) 
       
1.parent() 
       

El anillo de los números enteros se denota en Sage como ZZ.

ZZ 
       

Tal y como definimos las clases anteriormente, cada clase puede heredar de una única clase. Sin embargo en muchos casos hay clases que se obtienen por la intersección de dos clases, es decir, que deberían heredar de dos de ellas. Para ello Python dispone de lo que se denomina "categorías". Como ejemplo, podemos ver las categorías a las que pertenece el anillo de los números enteros y mostrar en forma de grafo sus relaciones:

ZZ.categories() 
       
EuclideanDomains().category_graph().plot(talk = True) 
       

Incrementando la eficiencia

Estadísticas del desempeño de un programa

En Sage utilizamos muchos tipos complejos de datos y métodos asociados a ellos. Esto hace que en ocasiones sea complejo determinar la eficiencia de un algoritmo. Las estadísticas de la ejecución de los algoritmos nos darán información de qué puntos del programa son los que consumen más tiempo y, por tanto, dónde debemos mejorar la eficiencia.

Vamos a comenzar con un ejemplo de cálculo de la constante de Brun. La constante de Brun se define como:

\[\sum_{p,p+2\in\mathbb{P}} \left(\frac{1}{p}+\frac{1}{p+2}\right),\]

donde $\mathbb{P}$ es el conjunto de los números primos.

#Suma de los inversos de los primos gemelos #1: Encuentra primos def criba(ls): '''Se queda con los elementos irreducibles de una lista de enteros''' primos = [] while ls: p = ls[0] primos.append(p) ls = [k for k in ls if k%p] return primos def lista_primos(K): 'genera los numeros primos menores que K' return criba(range(2,K)) #2: Selecciona los gemelos (Nos quedamos con el menor de cada par) def criba_gemelos(ls): '''recibe una lista de primos, y devuelve los numeros p tales que p+2 tambien esta en la lista''' return [p for p in ls if p+2 in ls] #3: Sumamos los inversos para aproximar la constante de Brun def brun(K): '''Devuelve la suma de los inversos de los primos gemelos menores que K''' primos = lista_primos(K) gemelos = criba_gemelos(primos) return sum( (1.0/p + 1.0/(p+2)) for p in gemelos) 
       
%time print brun(1e4) 
       
#importamos los modulos cProfile y pstats para ver las estadisticas de cuanto tiempo se pasa en cada parte del codigo import cProfile, pstats #No necesitamos entender la siguiente linea: tomalo como una version avanzada de timeit cProfile.runctx("brun(10000)", globals(), locals(), DATA + "Profile.prof") s = pstats.Stats(DATA + "Profile.prof") #Imprimimos las estadisticas, ordenadas por el tiempo total s.strip_dirs().sort_stats("time").print_stats() 
       

Vemos que la mayor parte del tiempo los dedica a criba gemelos y criba. Es sencillo ver que se puede mejorar la criba de gemelos pues dos números primos gemelos aparecerán juntos en la lista.

#2: Selecciona los gemelos. Nos quedamos con el menor de cada par def criba_gemelos(ls): return [ls[j] for j in xrange(len(ls)-1) if ls[j+1]==ls[j]+2] 
       
%time print brun(1e4) 
       

Ahora podemos ver que todo el tiempo lo dedica a la criba, que sería el siguiente punto a optimizar (si estás muy interesado, en las referencias está la optimización completa)

import cProfile, pstats cProfile.runctx("brun(50000)", globals(), locals(), DATA + "Profile.prof") s = pstats.Stats(DATA + "Profile.prof") s.strip_dirs().sort_stats("time").print_stats() 
       

Compilando el código en cyton

Una segunda manera (complementaria) de obtener algoritmos más rápidos es compilar el código. Comenzaremos por un código en Python para el cálculo de una integral por el método del trapecio compuesto.

def f(x): return sin(x**2) def integral(a, b, N): dx = (b-a)/N s = 0 for i in range(N): s += f(a+dx*i) return s * dx 
       
time integral(0.0, 1.0, 500000) 
       

Para compilar el código, hemos de escribirlo en Cython. Para ello escribimos %cyton al principio de la celda. Sage entiende que todo lo que escribamos en la celda será Cython. Además hemos de importar la función seno y cambiar ^ por **

%cython #Tenemos que importar la funcion seno from math import sin def f(x): return sin(x**2) def integral_cy1(a, b, N): dx = (b-a)/N s = 0 for i in range(N): s += f(a+dx*i) return s * dx 
       
time integral_cy1(0.0, 1.0, 500000) 
       

Se puede ver que el tiempo de ejecución se ha reducido a la mitad. Sin embargo, podemos mejorarlo sustancialmente indicando el tipo de datos.

%cython from math import sin def f(double x): return sin(x**2) def integral_cy2(double a, double b, int N): cdef double dx = (b-a)/N cdef int i cdef double s = 0 for i in range(N): s += f(a+dx*i) return s * dx 
       
time integral_cy2(0.0, 1.0, 500000) 
       

Incluso podemos hacerlo más rápido llamando directamente a librerías en C.

Por último, en las referencias se muestra cómo aplicar el mismo proceso para calcular el conjunto de Mandelbrot.

%cython import numpy cimport numpy #para declarar los tipos de los arrays tb tenemos que usar cimport def mandelbrot_cy5(float x0,float y0,float side,int N=200, int L=50, float R=3): '''returns an array NxN to be plotted with matrix_plot ''' cdef double complex c, z, I cdef float delta cdef int h, j, k cdef numpy.ndarray[numpy.uint16_t, ndim=2] m m = numpy.zeros((N,N), dtype=numpy.uint16) I = complex(0,1) delta = side/N for j in range(N): for k in range(N): c = (x0+j*delta)+ I*(y0+k*delta) z=0 h=0 while (h<L and z.real**2 + z.imag**2 < R*R): z=z*z+c h+=1 m[j,k]=h return m 
       
time m=mandelbrot_cy5(-0.59375, 0.46875, 0.046875,600,160) matrix_plot(m).show(figsize=(8,8)) 
       

Sage y C

En el apartado anterior, hemos comentado que podemos llamar a librerías en C. Vamos a ilustrar esto con un ejemplo.

Para cargar una función (o un "script") almacenado en fichero independiente podemos usar load o attach.

load DATA+"example.sage" 
       
attach DATA+"example.sage" 
       

Podemos cargar directamente una función en Sage que haga referencia a otra en C:

attach DATA+"test.spyx" 
       
test(10) 
       

Un poco de programación funcional

La programación funcional consiste en partir de una lista base e ir aplicando funciones (transformaciones) a esa lista hasta obtener la solución al problema. Hemos visto las listas por comprensión que sería una primera aproximación al problema. Sin embargo esa solución es muy ineficiente. Veamos una aproximación de Python a la programación funcional.

Vamos a crear una función que otiene el segundo número primo entre 10000 y 100000.

%time ls = srange(10000,100000) primos = [t for t in ls if is_prime(t) ] print primos[1] 
       

El problema de este algoritmo es que hemos generado todos los números primos entre 10000 y 100000 para luego sólo usar 38 de ellos. Para evitar generar todos esos primos, sustituimos srange por xsrange.

srange(100) 
       
xsrange(100) 
       

Podemos ver que xsrange(100) no devuelve ningún número. Lo que nos devuelve es un método para ir generando dichos números.

xlista = xsrange(100) 
       
print xlista.next() 
       
print xlista.next() 
       

Veamos un ejemplo de uso de xsrange:

acumulador = 0 for j in xsrange(101): acumulador += j acumulador 
       

También podemos generar listas por comprensión con xsrange. La diferencia es que tendremos que utilizar paréntesis en lugar de corchetes y que en lugar de darnos la lista, lo que nos da es un objeto denominado "iterador", que va generando los elementos de la lista según se van necesitando.

genera_cuadrados = (x^2 for x in xsrange(10)) for c in genera_cuadrados: print c 
       

Podemos ahora volver al programa original y hacer el cambio:

%time ls = xrange(10000,100000) primos = (t for t in ls if is_prime(t) ) primos.next() print primos.next() 
       

De modo alternativo:

%time contador = 2 for j in xsrange(10000,100000): if is_prime(j): contador -= 1 if contador == 0: break print j 
       

Generación de iteradores

Aunque xsrange nos ofrece una manera de generar listas "bajo demanda", en programación funcional es muy común necesitar listas "infinitas", pues a priori no se conoce el número de elementos que habrá que recorrer. Para generar iteradores infinitos, tenemos la palabra clave yield. Esta instrucción es similar al return con la diferencia de que no interrumpe la evaluación de la función: devuelve un elemento y se queda en ese punto de la función esperando a que le pidan otro elemento. En ese caso, continúa la evaluación por donde estaba.

def fib(): a, b = 0, 1 while 1: yield b a, b = b, a+b 
       
fibo=fib() 
       
for i in range(10): print fibo.next() 
       

Sage y Lisp

Como no podía ser de otra manera, si nos vemos necesitados de un lenguaje puro de programación funcional, Sage nos proporciona una interfaz para varios de ellos. El siguiente ejemplo muestra cómo ejecutar código Lisp dentro de Sage.

lisp.eval('(defun fib (n) (if (< n 2) n (+ (fib (- n 1)) (fib (- n 2)))))') 
       
lisp(10).fib()