Comprobando presencia de subconjuntos dentro de conjunto de forma eficiente

Hoy estaba trabajando con fechas y tenía que encontrar subconjuntos de fechas que estaban presentes o no en otros conjuntos. Numpy nos ofrece varias funciones de ayuda para comprobar este tipo de cosas. Ahora mismo se me vienen a la cabeza numpy.in1d y numpy.setdiff1d. El tema es que necesitaba hacer muchas operaciones de este tipo y he estado buscando la forma de poder optimizar el código un poco y resulta que, dependiendo de la operación, numpy puede resultar hasta 10 veces más lento (según mis pruebas de hoy) que usar puro CPython3. Veamos un ejemplo sencillo:

1) Imaginemos que mis datos son arrays de numpy y quiero conocer los datos de un array a que no están presentes en b:

import numpy as np
a = np.arange(1000000)
b = np.arange(10, 1000010)
%timeit np.setdiff1d(a, b)

Usando IPython y en mi máquina el anterior código me da el siguiente resultado:

10 loops, best of 3: 148 ms per loop

2) Si ahora hago lo mismo usando conjuntos (sets) obtendría el siguiente resultado:

a = set(np.arange(1000000))
b = set(np.arange(10, 1000010))
%timeit a.difference(b)

De la misma forma, usando IPython y en mi máquina el anterior código me da el siguiente resultado:

10 loops, best of 3: 30.2 ms per loop

Para este caso concreto numpy es ¡¡¡hasta 5 veces más lento!!!

Esto es una micro-optimización que me ha servido a mi en un caso concreto y si aplica espero que os resulte útil pero ya sabéis que la optimización prematura es la fuente de todos los males 🙂

Saludos

Actualización:

Como bien comentan tanto Jaime como Chema más abajo, la comparación es tramposa ya que, normalmente, lo que mido no es lo que se hace en un programa completo, se debería de medir el proceso completo. La pretensión de este post era baja y solo mostrar la diferencia de implementaciones entre lo que hace numpy y lo que hace CPython. Os recomiendo que leáis los comentarios para aprender más.

Gracias, Jaime y Chema.

PFS: Multiplicación de Karatsuba

Vamos a seguir con nuestra serie sobre 'python, fuente de sabiduria' hablando sobre la multiplicación.

¿Te has preguntado alguna vez cómo se multiplican dos números en un ordenador o calculadora? Normalmente, los circuitos de los procesadores incluyen operaciones básicas como la multiplicación de enteros de tamaño pequeño. Sobre estas operaciones básicas se construye todo lo demás.

El caso que nos ocupa hoy, multiplicar números enteros muy grandes, se consigue mediante el uso de algoritmos vía software. Una forma sencilla de calcular enteros es usando la misma metodología que usábamos en la escuela, con la tabla de multiplicar del uno al diez y con sumas:

    23958233
        5830 ×
------------
    00000000 ( =      23,958,233 ×     0)
   71874699  ( =      23,958,233 ×    30)
 191665864   ( =      23,958,233 ×   800)
119791165    ( =      23,958,233 × 5,000)
------------
139676498390 ( = 139,676,498,390        )

(ejemplo extraído de wikipedia)

El algoritmo que hace uso de esa metodología se suele llamar long, grade-school, naïve,..., multiplication. En este caso vamos a hacer una implementación sucia de este algoritmo descomponiendo el 'multiplicador' en dígitos que multiplicaremos al 'multiplicando':

def naive(x, y):
    total = 0
    yy = [int(i) for i in str(y)]
    for i, yyy in enumerate(yy[::-1]):
        total += x * yyy * (10 ** i)
    return total

Para el ejemplo anterior (multiplicando = 23958233 y multiplicador = 5830) lo que estariamos haciendo sería lo siguiente:

  • el multiplicador lo descomponemos en una lista de los dígitos que lo componen, [5, 8, 3, 0]
  • el 0 lo multiplicamos por el multiplicando
  • el 3 lo multiplicamos por 10 y por el multiplicando y lo sumamos al resultado del paso anterior
  • el 8 lo multiplicamos por 100 y por el multiplicando y lo sumamos al resultado del paso anterior
  • el 5 lo multiplicamos por 1000 y por el multiplicando y lo sumamos al resultado del paso anterior
  • devolvemos el resultado final

Esta es una forma de hacer multiplicaciones pero cuando los números involucrados en la multiplicación son muy grandes podemos encontrarnos con problemas de rendimiento y perder mucho tiempo haciendo un cálculo.

Como en el desarrollo de Python hay gente muy inteligente, si buceas por el código fuente de Python verás que para multiplicación de long's se usa el algoritmo de la multiplicación rápida de Karatsuba. ¿Y quién demonios es este japonés? Pues este japonés es un ruso al que, con 23 añitos (en 1960), se le ocurrió llevar la contraria al mismísimo Kolmogorov.

Lo mejor es que le echéis un ojo a la siguiente entrada (o en la Wikipedia en español) para entender más claramente el funcionamiento (DRY). Básicamente, de lo que se trata es de minimizar el número de multiplicaciones (aunque no siempre se consigue) usando una estrategia de 'divide y vencerás'.

A continuación os dejo una implementación muy compacta que he encontrado en project nayuki (si no entendéis algo del código preguntad en los comentarios):

## http://nayuki.eigenstate.org/res/karatsuba-multiplication/karatsuba.py
## http://nayuki.eigenstate.org/page/karatsuba-multiplication
_CUTOFF = 1024
def kara(x, y):
    if x.bit_length() <= _CUTOFF or y.bit_length() <= _CUTOFF:  # Base case
        #print('naive')
        return naive(x, y)
    else:
        #print('karatsuba')
        n = max(x.bit_length(), y.bit_length())
        half = (n + 32) // 64 * 32
        mask = (1 << half) - 1
        xlow = x & mask
        ylow = y & mask
        xhigh = x >> half
        yhigh = y >> half
        a = kara(xhigh, yhigh)
        b = kara(xlow + xhigh, ylow + yhigh)
        c = kara(xlow, ylow)
        d = b - a - c
        return (((a << half) + d) << half) + c

Bien, pues vamos a ver si esto es capar de optimizar mi función naive. Para ello vamos a multiplicar el siguiente número (tened cuidado si vuestro ordenador no es muy potente ya que el cálculo se puede llegar a alargar):

num = 1234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890*(10**14000)
print(num)

1234567890123456789012345678901234567890123456789012345678901234567890
1234567890123456789012345678900000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000

( Qué post más largo... 🙂 )

Usando la función naive.

a = naive(num, num)
%timeit naive(num, num)

1 loops, best of 3: 1.31 s per loop

Usando la función kara

b = kara(num, num)
%timeit kara(num, num)

10 loops, best of 3: 172 ms per loop

print(a == b)
print(1.31 / 0.173)

True
7.572254335260117

Vemos que hemos conseguido una ganancia de 7.5x!!!

Aunque si comprobamos como lo hace Python vemos que no hemos hecho una gran optimización (aunque no era lo importante de esta entrada):

%timeit num * num

1000 loops, best of 3: 455 µs per loop

Parece que solo es alrededor de 375 veces más rápido 😉

Este código lo he probado en Python 3.3, si veis alguna errata o alguna barbaridad que haya podido decir, por favor, avísadme.

Esta entrada la puedes descargar en formato notebook desde nuestro repositorio en github.

P.D.: Pretendía que esta entrada fuera más rigurosa y completa pero mi tiempo se ha visto disminuido de forma drástica últimamente. La familia manda 🙂

PFS: Exponenciación binaria

Hoy voy a iniciar una serie nueva de entradas que se llamará 'python, fuente de sabiduría' (PFS para abreviar). ¿Qué pretende esta serie? La idea básica es mostrar la influencia y realimentación entre las matemáticas y la programación. La idea e inspiración para estas entradas salió después de leer esta entrada.

Inauguraremos esta entrada con una forma de calcular la potencia entera de un número entero. Si tuviera que crear una función que calculase esto, la primera idea que se me ocurre es usar la fuerza bruta (soy así de primario)...

def pow_fuerzabruta(x,n):
    y = x
    for i in range(1,n):
        y = y * x
    return y

Con la anterior función, la cual no he pensado mucho, lo que estamos haciendo son n-1 multiplicaciones. A medida que aumente el exponente aumentará el tiempo de cálculo. Veámoslo en una gráfica:

from matplotlib import pyplot as plt
n = [5, 25, 100, 1000, 10000, 100000]       # exponente
t = [1.13, 3.93, 17.8, 330, 15000, 1260000] # t en microsegundos
plt.xscale('log')
plt.yscale('log')
plt.plot(n,t)
plt.xlabel('valor del exponente n')
plt.ylabel('t en microsegundos')

pow_bruto

[Nota: la escala de los ejes en la anterior figura es logarítmica]

Esto, cuando usamos números muy grandes, no es muy deseable. Por ello, vamos a mostrar el algortimo de hoy, la exponenciación binaria, también llamado algoritmo de izquierda a derecha binario, que nos permitirá reducir de forma drástica el número de operaciones a realizar. Su uso en Python lo podéis ver, por ejemplo, en la siguiente implementación de long_pow (ver a partir de línea 3376).

Podemos pensar que estas cosas son muy modernas y sofisticadas pero, en este caso, no es así. Los primeros orígenes de esta implementación para calcular potencias parece que sería anterior a 2200 años atrás. Según el libro de Donald Knuth titulado Seminumerical Algorithms, volume 2 of The Art of Computer Programming, página 441, podemos leer:

El método es bastante antiguo; apareció antes del año 200 A.C. en Pingala's Hindu classic Chandah-sutra [ver B. Datta y A.N. Singh, History of Hindu Mathematics 1, 1935]; sin embargo, parece que existen otras referencias a este método fuera de la India a lo largo de los 1000 años posteriores. Una discusión muy clara sobre como calcular 2n de forma eficiente para valores de n arbitrarios sería debida a al-Uqlidisi de Damscus en 952 D.C.; podéis ver The Arithmetic of al-Uglidisi por A.S. Saidan (1975), p. 341-342, donde las ideas generales están ilustradas para el caso de n=51. Ver también al-Biruni's Chronology of Ancient Nations (1879), p. 132-136; durante el siglo XI, los avances del mundo árabe tuvieron una gran influencia.

Más información histórica sobre los orígenes la podéis encontrar en este artículo.

Vamos a ver como funciona, más o menos, este algoritmo. La idea básica detrás del mismo es la siguiente (ver wikipedia para una explicación más completa):

Las potencias más rápidas de calcular son las que tienen la siguiente forma \({x^2}^n\). Estos es debido a que se pueden encontrar las potencias calculando el cuadrado del número de forma repetida.

Vamos a explicar el mismo ejemplo que al-Uqlidisi de Damscus (unos diez siglos más tarde...) considerando un número elevado a 51. Por tanto, puedo calcular \(5^2\) con solo una multiplicación, \(5^4={5^2}^2\) con solo dos multiplicaciones, \(5^8={{5^2}^2}^2\) en tres multiplicaciones,...

Creo que vais viendo por donde van los tiros. El algoritmo explicado de forma muy básica sería para el ejemplo \(5^{51}\):

  • Obtener el exponente como suma de potencias de 2

En este caso sería 51 = 32 + 16 + 2 + 1

  • Encontrar la base de potencias de dos que aparecen en el exponente

\(5^1=5\)

\(5^2 = 5 cdot 5 = 25\)

\(5^4= 25 cdot 25 = 625\)

\(5^8= 625 cdot 625 = 390625\)

\(5^{16} = 390625 cdot 390625 = 152587890625\)

\(5^{32} = 152587890625 cdot 152587890625 = 23283064365386962890625\)

  • Y, por último, se multiplicarían de forma conjunta las potencias ya encontradas

\(5^{51} = 5^{32} cdot 5^{16} cdot 5^2 cdot 5 = 444089209850062616169452667236328125\)

De esta forma, en lugar de multiplicar 5 por si mismo 50 veces solo hemos de realizar 8 multiplicaciones.

El que vamos a implementar es un algoritmo de exponenciación binaria pero de derecha a izquierda. Podéis ver los detalles en la wikipedia.

def my_pow(x,n):
    r = 1
    y = x
    while n > 1:
        if n%2:
            r = r * y
        n = int(n/2)
        y = y * y
    r = r*y
    return r

Si calculamos los tiempos que hemos obtenido con la nueva función (my_pow) y los comparamos con la anterior implementación (pow_fuerzabruta) veremos la ganancia que hemos obtenido:

from matplotlib import pyplot as plt
n = [5, 25, 100, 1000, 10000, 100000]       # exponente
t_bruta = [1.13, 3.93, 17.8, 330, 15000, 1260000] # t  para la primera función
t_bin = [1.34, 2.52, 3.68, 19.8, 765, 32800] # t para la segunda función
plt.plot(n,t_bruta, label='bruta')
plt.plot(n,t_bin, label='bin')
plt.legend()
plt.xlabel('valor del exponente n')
plt.ylabel('t en microsegundos')

pow_bruto_bin

[Nota: Ahora, las escalas de los ejes son lineales]

Vemos que la ganancia, dependiendo del valor del exponente, es cada vez más importante a medida que aumenta el valor del mismo.

for i, t1, t2 in zip(n, t_bruta, t_bin):
    print('Ganancia para n = {0:>6}, {1:>10.2f}'.format(i, t1/t2))


Ganancia para n = 5, 0.84
Ganancia para n = 25, 1.56
Ganancia para n = 100, 4.84
Ganancia para n = 1000, 16.67
Ganancia para n = 10000, 19.00
Ganancia para n = 100000, 38.00

Para números pequeños vemos que incluso puede llegar a ser peor pero para números grandes la ganancia es significativa llegando a valores en torno a 40 para un número n de orden \(10^5\)

Por último, solo comentar que las funciones `pow` y `math.pow` son más eficientes que lo que he hecho (programadas en C). La idea de la entrada solo es mostrar que, gracias al trabajo y saber colectivo, la inteligencia que se esconde en las tripas de Python (y el resto de lenguajes de programación) es mucha y en el día a día ni nos damos cuenta de ello.

P.D.: Las funciones implementadas no contemplan exponentes negativos ni 0. He intentado simplificar el código al máximo para que se vea el bosque en lugar de los árboles.

Esta entrada la puedes descargar en formato notebook desde nuestro repositorio en github.

Nuestro primer algoritmo heurístico

Normalmente, cuando se trata de optimizar 'algo' se pueden usar diferentes aproximaciones.

Vamos a proponer un problema sencillo y lo vamos a resolver de forma aproximada usando un algoritmo heurístico de 'búsqueda en escalada' o 'hill-climbing'.

Vamos a ello, la resolución del problema es obvia en este caso, pero solo se trata de que veáis el funcionamiento de todo esto. Tenemos un conjunto de n números que se hallan dentro de un subconjunto de los números reales. ¿Cuál sería el conjunto de estos n números que me permiten minimizar lo siguiente?

\(sum_{k=1}^N k^2\)

Para simplificar permitiremos que los n números se puedan repetir por lo que la respuesta obvia es que el conjunto de los números sea [0, 0, 0, ..., 0] por lo que la solución exacta sería 0. Si pensamos por un momento que la solución no es obvia, la búsqueda exhaustiva nunca sería óptima puesto que estamos buscando entre los números reales. Por tanto, ¿como resolveríamos esto con una búsqueda heurística? Primero de todo pondremos un poco de código, lo explicamos y luego haremos uso del mismo para ver los resultados que obtenemos:

import numpy as np
import matplotlib.pyplot as plt
plt.ion()
class Busqueda:
    u"""
    Clase (fea) para explicar el funcionamiento de
    un algoritmo heurístico
    """
    def __init__(self, ngeneraciones, ngenes, lim_sup, lim_inf):
        self.ngeneraciones = ngeneraciones
        self.ngenes = ngenes
        self.lim_sup = lim_sup
        self.lim_inf = lim_inf
    def padre(self):
        self.individuo0 = np.random.rand(self.ngenes) *
                          (self.lim_sup - self.lim_inf) + self.lim_inf
        return self.individuo0
    def fitness(self, individuo):
        funcion = individuo * individuo
        return np.sum(funcion)
    def hijo(self, individuo):
        nindividuo = individuo + np.random.normal(0,1, self.ngenes)
        nindividuo[nindividuo < self.lim_inf] = self.lim_inf
        nindividuo[nindividuo > self.lim_sup] = self.lim_sup
        return nindividuo
    def calculos(self):
        fit_acum = []
        individuo = self.individuo0
        poblacion = self.individuo0
        fit_acum = np.append(fit_acum, self.fitness(individuo))
        for generacion in xrange(self.ngeneraciones - 1):
            nindividuo = self.hijo(individuo)
            if self.fitness(nindividuo) <= fit_acum[-1]:
                individuo = nindividuo
            fit_acum = np.append(fit_acum, self.fitness(individuo))
            poblacion = np.append(poblacion, individuo)
        return poblacion.reshape(self.ngeneraciones, self.ngenes), fit_acum
    def representa_proceso(self, poblacion, fitness_acumulado):
        plt.subplot(311)
        plt.xlim(self.lim_inf, self.lim_sup)
        plt.ylim(0., self.lim_sup ** 2.)
        for gen in range(self.ngenes):
            plt.text(poblacion[0, gen],
                     self.fitness(poblacion[0, gen]),
                     gen+1)
        plt.subplot(312)
        plt.xlim(self.lim_inf, self.lim_sup)
        plt.ylim(0., self.lim_sup ** 2.)
        for generacion in range(self.ngeneraciones):
            for gen in range(self.ngenes):
                plt.text(poblacion[generacion, gen],
                         self.fitness(poblacion[generacion, gen]),
                         gen+1)
        plt.subplot(313)
        plt.xlim(0, self.ngeneraciones)
        plt.ylim(0., np.max(fitness_acumulado))
        plt.plot(fitness_acumulado)

Vamos a guardar el anterior código en un fichero que se llame heuristico.py y vamos a abrir ipython en el mismo lugar donde se encuentra el fichero heuristico.py.

Hemos creado una clase llamada Busqueda que hemos de inicializar con una serie de parámetros, el primero sería el número de iteraciones que vamos a usar en nuestra búsqueda, el segundo sería el número de miembros que vamos a usar, n, y el tercer y cuarto serían el límite superior e inferior del intervalo de números reales donde vamos a restringir la búsqueda. Dentro de ipython hacemos lo siguiente:

In [1]: from heuristico import Busqueda
In [2]: busqueda = Busqueda(1000, 3, 100, -100)

Vamos a hacer una simulación que contará con 1000 generaciones (o pasos, o iteraciones, como más os guste), la solución tendrá en cuenta un conjunto de 3 números que estarán restringidos entre [-100, 100].

Lo siguiente será crear la primera solución (la solución padre de la que derivarán todas las soluciones hijas siguientes) que no es más que una solución aleatoria para empezar la búsqueda. Lo hacemos de la siguiente forma:

In [3]: busqueda.padre()
Out[3]: array([-56.14070247,  96.07815303,  31.11068139])

Veis que el primer miembro ha tomado un valor en torno a -56, el segundo un valor en torno a 96 y el tercero un valor en torno a 31.

Ahora vamos a ver la chicha del algoritmo. El cálculo completo se realiza en el método calculos, que a su vez hace uso de los métodos fitness e hijo. El método fitness es el que calcula la solución del problema (en este caso se trata de minimizar el valor). El método hijo es el que crea una posible solución nueva cercana a la solución actual (solución padre). En el método calculos hacemos las iteraciones (1000 en este ejemplo) para comparar al padre con el hijo y ver si tiene mejor fitness y en caso de que sea así actualizamos al padre con el hijo, es decir, le pasamos los valores del hijo al padre.

In [4]: pob, fit_acum = busqueda.calculos()

Hemos calculado la población final, pob, el conjunto de 3 elementos que nos minimizan el resultado, y hemos obtenido el valor de la función fitness en cada iteración, fit_acum.

Por último, vamos a hacer una representación de los resultados de la siguiente manera:

In [5]: busqueda.representa_proceso(pob, fit_acum)

Que nos mostraría la siguiente ventana:

Donde el primer panel (arriba), nos muestra la posición de la primera solución aleatoria que hemos generado con el método padre. Como hemos visto anteriormente, el primer miembro ha tomado un valor en torno a -56, el segundo un valor en torno a 96 y el tercero un valor en torno a 31. En el segundo panel (en medio) vemos como va evolucionando la solución, cada uno de los miembros de la solución va descendiendo poco a poco hacia la solución final (cercana a cero). Por último, en el tercer panel (abajo) vemos el valor de nuestra solución y como se va a cercando a una solución cercana a 0, de hecho, vemos que a partir de las 300 iteraciones ya estaríamos cerca de 0. Aunque después de las 1000 iteraciones no hemos alcanzado la solución exacta, 0, puesto que estos algoritmos dan valores aproximados.

Veamos un segundo caso:

In [1]: from heuristico import Busqueda # Solo necesario si habéis cerrado ipython
In [2]: busqueda = Busqueda(1000, 20, 200, 0)
In [3]: busqueda.padre()
Out[3]:
array([  71.65433435,   13.38201722,   57.30044929,   74.69591418,
197.39440847,   42.3185271 ,   69.08520774,    7.40961601,
79.87846973,  161.14232418,  135.47599187,  183.17622057,
30.69728341,   21.85297791,  157.85448115,  139.02125074,
98.8888481 ,   49.02073966,   46.10618954,  116.69194813])
In [4]: pob, fit_acum = busqueda.calculos()
In [5]: busqueda.representa_proceso(pob, fit_acum)

Que nos daría el siguiente resultado:

En este caso vemos hemos usado un conjunto de 20 números y observamos que quizá serían necesarias más iteraciones para llegar a un número más aproximado.

Normalmente, el algoritmo se para cuando se han superado las iteraciones o cuando la solución no ha mejorado en un número determinado de pasos (esto segundo no lo he implementado por lo que lo tenéis como ejercicio).

Saludos.

P.D.: Como siempre, se agradece cualquier corrección, crítica constructiva, propuesta de mejora, debate,..., en los comentarios.

¿Cómo encontrar el mínimo de una función usando scipy?

NOTICE: Esta entrada es una traducción libre (algunas cosas no serán traducidas por no existir una traducción sencilla en castellano) de un artículo publicado en The Glowing Python con permiso de su autor.

Para la siguiente entrada se ha usado python 2.7.2, numpy 1.6.1, scipy 0.9.0 y matplotlib 1.1.0

En este ejemplo veremos como usar la función fmin para minimizar una función. La función fmin se encuentra en el módulo optimize de la librería scipy. La función fmin usa el algoritmo downhill simplex para encontrar el mínimo de la función objetivo empezando por un punto inicial dado por el usuario. En el ejemplo emezaremos a partir de dos puntos iniciales diferentes para comparar los resultados.

import numpy
import matplotlib.pyplot as plt
from scipy.optimize import fmin
# Función objetivo
rsinc = lambda x: -1 * numpy.sin(x)/x
# Empezamos a partir de x = -5
x0 = -5
xmin0 = fmin(rsinc,x0)
# Empezamos a partir de x = -4
x1 = -4
xmin1 = fmin(rsinc,x1)
# Dibujamos la función
x = numpy.linspace(-15,15,100)
y = rsinc(x)
plt.plot(x,y)
# Dibujo de x0 y el mínimo encontrado empezando en x0
plt.plot(x0,rsinc(x0),'bd',xmin0,rsinc(xmin0),'bo')
# Dibujo de x1 y el mínimo encontrado empezando en x1
plt.plot(x1,rsinc(x1),'rd',xmin1,rsinc(xmin1),'ro')
plt.axis([-15,15,-1.3,0.3])
plt.show()

La función fmin escribirá algunos detalles sobre el proceso iterativo llevado a cabo (dejamos la salida en inglés, que es lo que encontraréis cuando corráis el ejemplo):

Optimization terminated successfully.
         Current function value: -0.128375
         Iterations: 18
         Function evaluations: 36
Optimization terminated successfully.
         Current function value: -1.000000
         Iterations: 19
         Function evaluations: 38

Y la siguiente gráfica aparecerá en la pantalla:

Uso de la función fmin en scipy.optimize

El punto azul es el mínimo encontrado empezando a partir del diamante azul (x=-5) y el punto rojo es el mínimo encontrado a partir del diamante rojo (x=-4). En este caso, cuando empezamos a partir de x=-5 fmin se 'atasca' en un mínimo local mientras que si empezamos a partir de x=-4 fmin alcanza el mínimo global.

[Thanks to @JustGlowing to allow us to translate their articles]

Saludos.