Este artículo estuvo en mi lista de borradores (draft posts) por mucho tiempo ya. Evidente mente no cumplí con la consigna release early del desarrollo tipo bazar. Enough!

Este artículo va a hablar un poco de simulación de sistemas discretos, va a presentar un módulo que escribí como prueba de que se podía hacer y me va a servir para contarles algunas cosas uq e aprendí en el camino.

Sigo Aprendiendo Python y esta es una clase más.

Programación temporal de eventos

En el primer cuatrimestre de este año cursé una materia llamada Simulación. Estudiamos principalmente como simular sistemas discretos (como un grupo de pasajeros subiendo a un avión. Suben de a uno, en un momento hay 7 pasajeros a bordo y en otro momento hay 8, nunca hay 7.5) en contraposición a los sistemas continuos (como el tanque de nafta del avión dónde su nivel baja instante a instante cuando el avión está volando).

Básicamente existen dos tipos de técnicas que se pueden utilizar para simular sistemas discretos, las de orientación al intervalo y las de orientación al evento. En las del primer tipo, cada cierto tiempo t se verifica de alguna forma si ha ocurrido un evento y se realiza el procesamiento asociado a este. Esto lleva a impresiciones y en ocasiones a un gasto innecesesario de cpu. En las técnicas de orientación al evento estos problemas no existen.

Una de las técnicas más primitivas de orientación al intervalo (no se si se usa en la actualidad) es la de Programación Temporal de Eventos. Básicamente se tiene un reloj que cuenta el tiempo de simulación y una lista de eventos ordenada en el tiempo. Uno debe escribir una función o procedimiento por cada uno de los eventos principales del sistema (por ejemplo la llegada de una pieza a una máquina) y una función principal que contenga el main loop de la simulación. Esta irá tomando eventos de la lista, avanzará el tiempo de simulación según corresponda y los ejecutará. Siempre se toma el primero y se lo ejecuta, si no hay más elementos en la lista termina la siumulación.

Los eventos son agregados a la lista con una primitiva llamada por ejemplo PPE (programar próximo evento), dónde además de indicar cual será el próximo evento, se debe indicar dentro de cuanto tiempo a partir del actual se debe ejecutar.

Mientras cursaba, estudiaba y rendía pensaba que sería muy fácil implementar en Python un módulo que provea las primitivas necesarias para realizar una simulción por programación temporal de eventos. Poner y sacar funciones de una lista seguro sería sencillo en un lenguaje muy dinámico y donde todo (incluidas las funciones) son objetos. La semana siguiente a rendir estuve escribiendo un poco de código Python. El resultado es un módulo llamado pte que cubre las funcionalidades mínimas que yo pretendía.

Resultado

No pretendía realizar un imponente desarrollo ni pienso que alguien use esto alguna vez en serío. Es más bien una prueba de concepto (Proof of concept). Lo había avandonado y ni siquiera funcionaba sin tirar errores. Objetivo: que funcione.

Este es el módulo: pte.py (versión orginal)

Una versión con el consejo de Daniel: pte-d.py (gracias Daniel!)

Y este un ejemplo en el que lo uso: pte_ej2.py

Ejemplo de ejecución:

juanjo@sarge:~/python$ ./pte_ej2.py



500 piezas fueron aceptadas en 1247.55 unidades de tiempo.



Máximo tamaño de la cola para procesamiento 5.

Qué aprendí

Módulos

Más de lo que pueda decir sobre módulos: http://docs.python.org/tut/node8.html

reload(modulo)

Reload a previously imported module. The argument must be a module object, so it must have been successfully imported before. This is useful if you have edited the module source file using an external editor and want to try out the new version without leaving the Python interpreter. The return value is the module object (the same as the module argument).

Entonces lo que hay que hacer cuando se está trabajando en la consola de Python[1] es, primer importar el módulo que se quiere probar:

>>> import modulo

Hacer algunas pruebas:

>>> modulo.var1

1

>>> modulo.cuad(3)

9

Si luego editamos el texto de modulo y queremos probar la función que agregarmos:

>>> modulo = reload(modulo)

Funciones como objetos

En Python todo es un objeto, incluso las funciones. Sabiendo esto y teniendo en mente lo que quería hacer (en la definición de una función, explicitar que esta se llame en el futuro, es decir almacenarla en una lista de eventos) hice una pequeña prueba en el intérprete:

>>> l = []

>>> def f():

    l.append(f)





>>> l

[]

>>> f()

>>> l

[<function f at 0x418bce9c>]

>>> l[0]

<function f at 0x418bce9c>

>>> l

[<function f at 0x418bce9c>]

>>> l[0]()

>>> l

[<function f at 0x418bce9c>, <function f at 0x418bce9c>]

Esto es posible justamente gracias a que en Python las funciones son objetos. En ciencias de la computación esto se conoce como Funciones de Primera Clase, es decir que las funciones sean Objetos de Primera Clase, es decir que puedan ser creadas durante la ejecución de un programa, almacenadas en estructuras de datos, recuperadas más tarde para usarlas, pasadas como argumentos y retornadas en otras funciones. Más al respecto en wikipedia.

Ordenar según un criterio propio

Una lista puede ser ordenada según un criterio propio pasándole al método sort de una lista ese criterio en forma de función:

def _comparar_tiempos(a,b):

    """ Función auxiliar usada para ordenar los pares eventos,tiempos

    en la lista de eventos.

    """

    t1 = a[1]

    t2 = b[1]

    if (t1 > t2): return 1

    elif (t1 < t2): return -1

    else: return 0

y luego..

eventos.sort(_comparar_tiempos)

Espero les haya resultado interesante! [1] También conocida como REPL, algo muy piola cuando uno esta desarrollando.

Este es el material de la charla que dí el sábado pasado (19 de agosto) en el 1º Python Day del Grulic:

<li><a href="https://viejoblog.juanjoconti.com.ar/files/python/fu/charla-cordoba/index.html">Presentación (navegable y link al formato original)</a></li>

<li>Código fuente</li>

• colorFrame.py
• sky.py

También tengo un texto que escribí para ordenar mis ideas, pero está muy desprolijo como para publicarlo. SI alguien lo quiere (es buena ayuda para seguir la presentación) pídamelo por mail y se lo mando.

El próximo sábado 19 se va a llevar a cabo en la Facultad Regional del Córdoba de la UTN el primer PythonDay del Grulic. Con unos amigos de la facu (César, Mariano y Pablito y Bruno!*) vamos a viajar para participar.

Yo preparé una pequeña charla sobre Python-Fu. Espero que salga bien :)

Nos vemos en Córdoba!

• Actualización a pedido del público :)

En el wiki del LUGLi ya está todo el material de las jornadas, audio, vidoe, fotos, slides y más! : http://www.lugli.org.ar/mediawiki/index.php/1%C2%BA_Jornada_Python_en_Santa_Fe

También hay un par de relatos sobre las jornadas: uno en el blog de Facundo y otro de Nubis :) en la lista de correos de pyar.

Saludos!

En el evento Python Santa Fe 2006 que se hizo este sábado entregamos a los asistentes un trifolio que en el interior tenía una cartilla introductoria al lenguaje de programación python. Dejo aca el archivo original hecho con OpenOffice (hay dos versiones de igual contenido pero la diagramación se veía distinta en diferentes máquinas) y una versión en pdf:

<li><a href="http://firebirds.com.ar/~juanjo/wordpress/files/python/cartilla/cartilla-python.tar.bz2">cartilla-python.tar.bz2</a></li>

• cartilla-python-0.1.1-a.odt
• cartilla-python-0.1.1-b.odt
• cartilla-python-0.1.1.pdf

En el intérprete interactivo de python:

>>> import this

The Zen of Python, by Tim Peters

Beautiful is better than ugly. Explicit is better than implicit. Simple is better than complex. Complex is better than complicated. Flat is better than nested. Sparse is better than dense. Readability counts. Special cases aren't special enough to break the rules. Although practicality beats purity. Errors should never pass silently. Unless explicitly silenced. In the face of ambiguity, refuse the temptation to guess. There should be one-- and preferably only one --obvious way to do it. Although that way may not be obvious at first unless you're Dutch. Now is better than never. Although never is often better than *right* now. If the implementation is hard to explain, it's a bad idea. If the implementation is easy to explain, it may be a good idea. Namespaces are one honking great idea -- let's do more of those!

Logré bajar la complejidad de mi solución de un orden cúbico a un orden cuadrático, la diferencia en cuanto a tiempo de ejecución requerido es sorprendente:

Pero todavía no e suficiente, con un N igual a 10000 la respuesta se obtiene luego de 10 minutos (en mi computadora de 700 Mhz).

En esta solución ya no uso la función pythagoras(x,y,z). Solo hay 2 ciclos for anidados y el valor de z se calcula como z = sqrt( pow(x,2) + pow(y,2) ). Si z es entero, se tiene un triplete válido y solo falta verificar si x,y,z son primos relativos.

def is_int(i):
return i == int(i)



def solv106f(N):
miembros = []

tripletes_rprim = 0

for x in range(1,N-1):

    for y in range(x+1,N):

        z = sqrt( pow(x,2) + pow(y,2) )

        if (is_int(z) and z <= N):

            for a in [x,y,z]:

                if a not in miembros:

                    miembros.append(a)

            if (rprim(x,y,z)):

                tripletes_rprim +=1



no_miembros = N - len(miembros)

return N, tripletes_rprim, no_miembros

Esta es mi solución al problema número 6 de una larga lista de problemas de programación. A pesar de que, como estudiante de Ingeniería, pasé los últimos 3 años y medio estudiando distintas formas de las matemáticas, no conocí el Último Teorema de Fermat hasta que mi amigo Joel, estudiante de Filosofía, me lo comentó.

El problema en realidad no tiene mucho que ver con Fermat, pero sí con Pitágoras: dado un valor N, hay que encontrar cuántos tripletes x,y,z satisfacen x² + y² = z² tales que

<li>x,y,z sean menores o iguales a N</li>

<li>x,y,z sean primos relativos (es decir que 1 es su único divisor común)</li>

<li>x &lt; y &lt; z</li>

Además hay que encontrar cuántos números mayores que 0 y menores o iguales a N no pertenecen a ninguno de los tripletes (no solo a los tripletes de primos relativos).

El enunciado original aquí.

Una posibilidad sería explorar exasutivamente todas las combinaciones de números x,y,z menores o iguales a N, y para cada uno verificar si cumplen las características pedidas. Pero esto sería de altísima complejidad, sería de orden cúbico. Hay que pensar una mejor solución.

Empecemos con una función que dado 3 valores x,y,z me diga si satisfacen el Teorema de Pitágoras:

from math import pow



def pythagoras(x,y,z):

    return (pow(x,2)+pow(y,2) == pow(z,2))

>>> pythagoras(3,4,5)

True

>>> pythagoras(1,1,1)

False

>>> pythagoras()

True

También necesito una función, que dado 3 valores x,y,z me diga sin son primos relativos entre sí:

def mcd(a,b): # algoritmo de euclides

    if (a == b): return a

    elif (a < b): a,b = b,a # a es el más grande

    while(b != 0):

        r = a % b

        a = b

        b = r

        return a



def rprim(x,y,z):

    return (mcd(mcd(x,y),z) == 1)

La solución que explora todas las combinaciones de 3 números constaría de 3 ciclos for anidados y en el curerpo del más profundo se haría las preguntas pertitenes a la resolución del problema: ¿ese triplete satisface el teorema de pitagoras? ¿son primos relativos? ¿x<y<z?

for(x=1;x=<N;x++)

    for(y=1;y=<N;y++)

        for(z=1;z=<N;z++)

            # ¿...?

Una de las condiciones es que x<y<z, esto permite una redefinición tal que NO se exploren todas las combinaciones de 3 números menores que N sino solo las combinaciones que cumplen con esta restricción:

for(x=1;x=<N-2;x++)

    for(y=x+1;y=<N-1;y++)

        for(z=y+1;z=<N;z++)

            # ¿..?

¿Esta solución es buena o no es mucho mejor a una en la que todos los ciclos for vayan de 1 a N y por cada uno se pregunte si x,y,z?

En el archivo pythagoras.py están definidas dos funciones: solv106s(N), basada en la primera aproximación y solv106(N) que espera ser mejor.

Este programita sirve para medir el tiempo de las dos ejecuciones y compararlos.

import timing



def test(N=1000):

    timing.start()

    solv106s(N)

    timing.finish()

    print "la primer solución tardo " + str(timing.seconds()) + " segundos"

    timing.start()

    solv106(N)

    timing.finish()

    print "la segunda solución tardo " + str(timing.seconds()) + " segundos"

Hago una prueba con un valor relativamente grande, N=1000 el valor por defecto:

>>> test()

la primer solución tardo  1879 segundos

la segunda solución tardo 1561 segundos

31 minutos contra 26 minutos.. no me parece una diferencia muy grande. Nota: mi máquina es un AMD K6 II de 700 Mhz.

Esta gráfica compara las dos soluciones luego de haberlas probado con diferentes valores de N, obtenido unos puntos e interpolar una curva representativa:

mmm, hay que pensar una mejor solución.

ACM (Association for Computing Machinery) organiza unas competencias de programación llamadas ICPC (International Collegiate Programming Contest). En esta dirección hay muchos problemas de las competencias: http://acm.uva.es/problemset/.

Este es el enunciado del primer problema de la guía: http://acm.uva.es/p/v1/100.html

El fin de semana hice esta solución:

(nota: los lenguajes soportados por las competencias son C y Java, mi solución está escrita en Python por que me resulta más ágil. De todas formas importa más el algoritmo que el lenguaje.)

def odd(n):

    return not (n%2 == 0)

def tnu(n=1):

    print n

    if (n == 1): return

    elif (odd(n)): n = 3*n+1

    else: n = n/2

    tnu(n) # la recursividad por la cola es evitable

>>> tnu(22)

22

11

34

17

52

26

13

40

20

10

5

16

8

4

2

1

La longitud del ciclo para un n dado es la cantidad de valores que la función tnu() develve (incluido el 1 final)

El problema pide leer un archivo con pares i,j y para cada par devolver la mayor longitud de ciclo para todos los enteres enter i y j incluyendo a i y a j.

Solución Exaustiva

Una solución exaustiva podría consistir en calcular para cada par i,j todos las longitudes de ciclo asociadas y devolver la mayor.

Redefino tnu para que en lugar de imprimir los sucesivos valores de n devulva una lista con los valores, el segundo argumento de la función es una lista vacía.

def tnu(n,l):

    #print n

    l.append(n)

    if (n == 1): return l

    elif (odd(n)): n = 3*n+1

    else: n = n/2

    return tnu(n,l) # la recursividad por la cola es evitable

Suponiendo que ya se leyeron los pares del archivo correspondiente y se tienen en una lista llamada pairs:

def solv1(pairs):

    for i,j in pairs:

        maxc = 0

        for n in range(i,j+1):

            maxc = max(len(tnu(n,[])),maxc)

        print i, j, maxc

>>> pairs = [(1,10),(88,123),(7,1),(55,455),(10000,10003)]

>>> solv1(pairs)

1 10 20

88 123 119

7 1 0

55 455 144

10000 10003 180

>>> pairs = [(1,10),(100,200),(201,210),(900,1000)]

>>> solv1(pairs)

1 10 20

100 200 125

201 210 89

900 1000 174

Bien, esta solución ANDA, pero no es óptima. De hecho probablemente sea la solución menos eficiente :-) la de la fuerza bruta.

Solución Dinámica

La idea es almacenar los resultados ya calculados para no tener que hacerlos una y otra vez. Si se pide la máxima longitud de ciclo (mlc) enter 100 y 200 y obtengo 125 y luego se pide la máxima longitud de ciclo enter 1 y 300, voy a calcular la mlc para el par 1,100, luego lo haré para 200,300 y luego compararé los resultados con 125 quedándome con el mayor de los 3.

Pero.. que pasa si ya he calculado la mlc para los valores entre 2 y 99?

def solv1din(pairs,dic):

    for i,j in pairs:

        maxc = f_maxc(i,j,dic)

        print i, j, maxc

        dic[i,j]=maxc

def f_maxc(i,j,dic):

    wide = 0

    for a,b in dic.keys():

        if a == i and b == j: return dic[a,b]

        if (i<=a and b<=j and b-a > wide):

            f,g = a,b

            wide = b-a

    if (wide == 0): #ningún rango del diccionario estaba dentro de i,j

        maxc = 0

        for n in range(i,j+1):

            maxc = max(len(tnu(n,[])),maxc)

        return maxc

    else:

        return max(f_maxc(i,f,dic),dic[f,g],f_maxc(g,j,dic))

>>> solv1din(pairs,{})

1 10 20

100 200 125

201 210 89

900 1000 174

Una versión no recursiva de tnu podría ser algo como esto:

def tnu(n):

    l = []

    while(n != 1):

        l.append(n)

        if (odd(n)): n = 3*n+1

        else: n = n/2

    l.append(1)

    return l

Download

• Una implementación completa que resuelve el problema: 3n+1.py
• Un archivo de entrada de ejemplo: input100.txt

Ejemplo de uso

juanjo@sarge:~/problemas$ ./3n+1.py input100.txt

1 10 20

100 200 125

201 210 89

900 1000 174

88 123 119

7 1 0

55 455 144

10000 10003 180

En el LUGLi estamos organizando una jornada de charlas sobre el lenguaje de programación Python para el 3 de junio. Nos estamos organizando por medio del wiki del lugli y hasta tenemos un sitio web oficial! ¿Ya te anotaste para venir? Te esperamos!