Max/Min heap, priority queues… ¿Qué son y cómo lo hago? Comparativa Python, C/C++.


Introducción al algoritmo. ¿Porqué es útil?


Heapsort es el nombre que se le da a un tipo de algoritmo de ordenación basado en trees o árboles y que opera con complejidad O(n\log_{}n). El algoritmo realiza se realiza in place, es decir, sólo un número constante de elementos del array se almacenan fuera del mismo durante la ejecución del algoritmo.

Este algoritmo hace uso de un tipo de estructura de datos llamada “heap” para manejar la información, como veremos más abajo el uso de heaps permite implementar una cola de prioridad muy eficiente.


Como dijo Jack…


Aviso a navegantes: huelga decirlo, pero al igual que en el resto de posts, agradeceré cualquier mejora o corrección sobre la entrada.

Una binary-heap (me tomaré la licencia de no usar la traducción del término como montículo binario) es un tipo de estructura de datos que trata de representar un array en forma de árbol binario semicompleto (el término anglosajón es nearly complete binary tree).

Cada elemento del array constituye un nodo del árbol, que tendrá todos sus nodos completos excepto, opcionalmente, el último de los nodos.

Un array A que representa nuestro heap tiene dos parámetros, a saber:

  • A.length: longitud del array, es decir, el número de elementos que podemos almacenar.
  • A.size: número de elementos contenidos en el array.

Si la raíz (root) del árbol se encuentra en A[1], dado el índice i de un nodo podremos calcular el índice de su padre y sus hijos (derecho e izquierdo).

  • Parent = (\frac{2}{i})
  • Left = 2\cdot i
  • Right =2\cdot i+1

Paremos por un momento a reflexionar sobre el hecho de que la raíz del árbol se encuentre en A[1]. ¿Porqué sucede esto cuando todos los lenguajes de programación que conozco comienzan a indexar los arrays a partir del índice 0? Podemos alegar, como yo hice al leer el algoritmo, que son cosas del pseudocódigo pero hay una razón más sutil: la optimización. El hecho de colocar nuestra raíz en A[1] y no en A[0] simplifica las fórmulas anteriores de forma notable (puedes deducirlas rápidamente). Puede parecer insignificante, pero no lo es cuando buscamos eficiencia.

Según la bibliografía, las implementaciones eficientes de algoritmos heapsort codifican estas operaciones como macros o procedimientos inline y calculan los índices a nivel de bit. Como veremos más abajo este hecho introduce una mejora en el rendimiento.

Como adelanta la introducción, uno de los usos más comunes de los heaps es construir una cola de prioridad (priority queue). Las hay de dos tipos; max y min, o más comúnmente max-heap o min-heap.

Ambos tipos de heaps satisfacen una propiedad llamada heap property (imaginativo, ¿verdad?) que queda resumida de la siguiente manera:

  • Max-heap: A[Parent(i) <= A[i]
  • Min-heap: A[Parent(i) >= A[i]

Esta regla no es más que un criterio de ordenación que determina si los valores se ordenan desde la raíz a las hojas de mayor a menor, max-heap, o viceversa, min-heap.

Ciertos funciones, o métodos, son necesarios para crear y mantener nuestra estructura. De ahora en adelante trataremos exclusivamente con min-heaps, pero tanto los procedimientos como las conclusiones serán análogas a las max-heaps. Vamos con ellos:

  • Min-Heapify: Una función que se encarga de mantener la propiedad de min-heap y que opera en tiempo O(log_{}n). Es el alma del algoritmo.
  • Build-Heap: en tiempo lineal genera el min-heap a partir de un array de entrada. Podemos pensar en ella como un constructor que recibe el array a ordenar como parámetro.
  • Heap-Sort: ordena el array in place en O(nlog_{}n)
  • Heap-insert, Heap-extract-min, Heap-increase-key: rutinas de complejidad O(log_{}n) que nos permiten interaccionar con la cola; extraer, insertar y modificar elementos.

Así contado puede no tener demasiado sentido, se comprende mejor con unos esquemas o un vídeo, podéis ver este:



Ahora ya queda claro que el principio básico del heap es mantener la propiedad de los heaps, es decir, que el nodo inferior siempre sea mayor que el padre en una cola de prioridad mínima y al revés en el caso de máxima prioridad.


¿Vamos ya con el código?

No he encontrado en internet mucho sobre la implementación de una priority queue en C, seguramente porque los lenguajes más modernos las incorporan en sus librerías estándar. Así que he tenido que implementarla casi desde cero mirando algunos ejemplos.

Las sorpresas las dejo para el final, así que primero vamos con C++ y con Python.

Como vamos a ver ahora, implementar los heaps en C++ es bastante sencillo, dejo aquí la implementación más rápida entre las que he probado, adaptada desde el código de cplusplus.com.

#include <queue>
#include <functional>
#include <ctime>
#include <iostream>
#include <math.h>

using namespace std;

int main(int argc, char* argv[])
{
    srand(time(0));

    int elements;
    int *values = (int*)malloc(sizeof(int)*elements);

    for (int j = 0; j < 9; ++j)
    {

        elements = pow(10,j);
        int *values = (int*)malloc(sizeof(int)*elements);
        // Generamos un array de elementos enteros
        for( int i = 0; i != elements; ++i ) {
            values[i] = rand();
        }

        // Construimos la cola a partir del array
        std::priority_queue<int, std::vector, std::greater > mypq (values, values+elements);

        // Extraemos los elementos uno a uno
        for (int i = 0; i < elements; ++i)
        {
            mypq.pop();
        }
    }
    return 0;
}


Y la de Python es aún más simple, como nos tiene acostumbrados:

import numpy as np
import heapq
import sys
import time

time_heapify = np.zeros(8)
time_extract = np.zeros(8)

for i in range(8):
    elements = pow(10,i)
    # Creamos una lista de enteros
    a = list(np.random.randint(low=sys.maxint, size=elements))
    # Transformamos el array de enteros en una heapq
    heapq.heapify(a)

    for j in range(elements):
        # Extraemos el elemento minimo
        heapq.heappop(a)


Resultados, una de carreras.


Los resultados son interesantes, no me esperaba que el programa en C++ distara tanto del programa en C. C marca la diferencia incluso sin optimizar.

3 lang. insert()

Se repite el mismo escenario que el caso anterior, C marca la diferencia…

3lengExtract


Optimizaciones implementadas sobre el programa en C.


He implementado algunas mejoras básicas sobre el programa en C y los resultados los podéis ver resumidos aquí:

extractc
insert() min-heap in c

  • El programa básico ya parte con su raíz en A[1], por lo que no podremos apreciar la reducción de tiempo que introduce esa mejora.
  • El siguiente paso ha sido sustituir las expresiones de Parent, Left y Right por sus operaciones a nivel bit (operadores >>, << y | ) e implementarlas como macros.
  • El tercer paso consiste en emplear registros para algunas variables que se usan a menudo y hacer algunas funciones clave inline.
  • El paso más significativo es la optimización agresiva del compilador, sólo para ver que marca la diferencia la he incluido aquí.


Notas finales.


Puede que los tiempos no reflejen el rendimiento que la aplicación pueda ofrecer en condiciones óptimas ya que las he realizado en mi ordenador personal, si todo va bien en un par de meses probaré la aplicación completa en el centro de computación y realizaré unos benchmarks más adecuados.

La versión de C definitiva la dejaré en Github en breves, tengo que terminar alguna prueba más! (;

Los datos sobre las versiones del compilador g++ empleado tanto para C como para C++:

i686-apple-darwin11-llvm-g++-4.2

Y la versión de Python; Python

Python 2.7.5 :: Anaconda 1.7.0 (x86_64)

Publicado en Algoritmos, C/C++, python | Etiquetado , , , , , , , , | Deja un comentario

Uso de numpy.optimize.leastsq() o cómo ahorrar trabajo.

Como sospechaba en la entrada de ayer, es posible ‘engañar’ a la función optimize.leastsq de numpy para trabajar con arrays n-dimensionales.
El resultado es mejor (por poco😉 ) al obtenido mediante el algoritmo de Gauss-Newton ayer, aunque no me cabe la menor duda de que será mucho más eficiente.

A trabajar.

Para que esta función haga su trabajo necesitaremos darle algunas cosas hechas:
  • Una función que calcule la imagen (matemática) de nuestra función. La llamaré
  • Una función que calcule el residual. La llamaré residual.
  • Un arras (no es obligatorio pero sí más cómodo) que contenga las variables a optimizar. Lo llamaré beta.
Con esto podremos poner a trabajar a nuestra función de la siguiente forma:

def GaussianFunct( beta, x, y):
 H,B,sigma,x0,y0 = beta

 exponente = -( (x-x0)**2 +(y-y0)**2) / (2*sigma*sigma)
 return H*np.e**(exponente) + B

def residual(beta, x,y,z):
 return (z-GaussianFunct(beta,x,y))**2

Ahora vamos con la optimización:

beta,cov,infodict,mesg,ier = optimize.leastsq( residual, beta, args=(meshX, meshY, Z), full_output=True)

Al tratarse mi caso de una función tridimensional (x,y,z) (no en el ámbito matemático, pero así nos entendemos) es cómodo evaluarla mediante ‘mallas’ obtenidas mediante meshgrid(), pero optimize.leastsq() no los acepta. No hay problema, las aplanamos primero mediante ravel() o reshape():

z = star[0:24,0:24]
x = np.arange(24)
y = np.arange(24)

meshx,meshy = np.meshgrid(x,y)

Z = np.ravel( z )
meshX = np.ravel( meshx )
meshY = np.ravel( meshy )

beta = np.array([ 100 ,6, 1.6,12, 12])

Resultados

El resultado es muy similar al de la entrada anterior en términos de error máximo porcentual, pero mirando con un poco más de detenimiento la distribución del error en la superficie no es la misma…

errorNumpy

Dejo aquí el notebook para que se vea más claro el flujo del programa😉.

Publicado en matemáticas, python | Etiquetado , , , , | 10 comentarios

Gauss-Newton en Python. Mínimos cuadrados no-lineales n-dimensionales.

El problema de optimización no lineal mediante mínimos cuadrados (no sujeto a restricciones) se enuncia formalmente como:
\displaystyle \underset{x}{minimizar}\, f(x)=\underset{i=1}{\overset{m}{\sum}}f_{i}(x)^{2}\equiv\frac{1}{2}F(x)^{T}F(x)
donde la función objetivo está definida en términos de funciones auxiliares {\left\{ f_{i}\right\} } y {F(x)} es una función vectorizada:
\displaystyle F(x)=(f_{1}(x)\: f_{2}(x)\:\cdots\: f_{m}(x))^{T}
En el contexto del ajuste de datos, las funciones auxiliares {\left\{ f_{i}\right\} } no son funciones arbitrarias, sino que corresponden los residuales del ajuste.
Las componentes del gradiente {\nabla f(x)} pueden calcularse como sigue:
\displaystyle \nabla f(x)=\nabla F(x)F(x)
De la misma forma el valor de la matriz Hessiana:
\displaystyle \nabla^{2}f(x)=\nabla F(x)\nabla F(x)^{T}+\underset{i=1}{\overset{m}{\sum}}f_{i}(x)\nabla^{2}f_{i}(x)
La matriz Hessiana es en muchas ocasiones complicada, cuando no imposible, de obtener. Tendremos que realizar {n^{2}} derivadas parciales dobles, por lo que para un problema como el que sigue de 5 variables tendremos que calcular 25 de estas operaciones. En cualquier caso su cálculo es pesado tanto si lo realizamos analíticamente como si aproximamos mediante diferencias finitas, aquí la potencia del método:
Llamemos {x_{*}} a la solución del problema, que causaría un ajuste perfecto de los datos, o lo que es lo mismo:
\displaystyle F(x)\approx0\: para\, x\approx x_{*}
Lo que nos lleva a:
\displaystyle \nabla^{2}f(x)=\nabla F(x)\nabla F(x)^{T}+\underset{i=1}{\overset{m}{\sum}}f_{i}(x)\nabla^{2}f_{i}(x)\approx\nabla^{2}f(x)=\nabla F(x)\nabla F(x)^{T}
\displaystyle
Por lo tanto este método sólo requiere el cálculo de las primeras derivadas para el cálculo del gradiente. Ahora sólo tendremos que resolver el sistema de ecuaciones alegraicas:
\displaystyle \nabla^{2}f(x)p=-\nabla f(x)
Implementación del algoritmo.
Mi implementación es un borrador para resolver un problema de optimización que consiste en lo siguiente, sólo para poneros en contexto. debo encontrar los valores para las constantes {H,\sigma,B,x_{o},y_{o}} que minimizen la distancia de la superficie definida por una campana doble de Gauss:
\displaystyle G=He^{-((x-x_{0})^{2}+(y-y_{0})^{2})\frac{1}{\sigma^{2}}}+B
respecto de los datos de brillo de una estrella, los cuales obtengo de una imagen directamente (podremos considerarlos a todos los efectos datos del problema).
Quien juegue un poco con numpy en el ámbito de la optimización puede venirle a la cabeza algunas funciones que proporciona, como la de mínimos cuadrados. No hay ningún problema con ellas, aunque no se si pueden trabajar con arrays n-dimensionales o mallas, pero el caso es que necesitaba tener control sobre el algoritmo y poder cambiar algunas partes.
Los pasos para resolver el algoritmo son los siguientes:
  • Hacer una estimación inicial de la semilla a la que llamaré {\beta}.
  • Calcular el gradiente, nabla, de la función G definida anteoriormente. Este calculo se hará analíniticamente porque siempre necesito aplicar la misma función.
  • Calcular la matriz Hessiana llamada “Hessian”.
  • Resolver el sistema de ecuaciones algebraicas para obtener unos nuevos pesos.
  • Sumar el vector de pesos a beta y realizar una nueva iteración hasta satisfacer el criterio de parada.
Vamos ahora con el código:
def Gaussian_distrib_3d(beta,x,y):
H,B,sigma,x0,y0 = beta
exponente = -( (x-x0)**2 +(y-y0)**2) / (2*sigma*sigma)
return H*np.e**(exponente) + B
Ahora vamos con el cuerpo del algoritmo:
gradientF = np.zeros([Nvariables,Ndata])

for step in range(0,N):
 H,B,sigma,x0,y0 = beta
 # Calculo de algunos parametros que se repiten
 p1 = (meshx-x0)**2 +(meshy-y0)**2
 p2 = 2*sigma*sigma
 p3 = sigma**3
 p4 = sigma*sigma
 exponente = -p1/p2

 Fx = np.ravel( H*np.e**(exponente) + B -z )
 Fx = Fx.transpose()

 # Calculo del gradiente
 gradientF[0,:] = np.ravel(np.e**(exponente))
 gradientF[1,:] = 1
 gradientF[2,:] = np.ravel((p1/p3)*H*np.e**(exponente))
 gradientF[3,:] = np.ravel(((meshx-x0)/p4) * H * np.e**(exponente))
 gradientF[4,:] = np.ravel(((meshy-y0)/p4) * H * np.e**(exponente))

 gradientFTransp = gradientF.transpose()

 nabla = np.dot(gradientF,Fx).reshape(Nvariables,1)

 #Calculamos la aproximacion de la matriz Hessiana
 Hessian = np.dot(gradientF,gradientFTransp)

 #Resolvemos el sistema algebraico
 p = np.linalg.solve(Hessian, -nabla)
 beta += p
El algoritmo se inicia tomando unos valores iniciales para la semilla, beta en nuestro caso:
beta = np.array([100,6,1.2,12,10.5]).reshape(Nvariables,1)
Por último los valores de x e y, puesto que se trata de una función ‘tridimensional’ (no en términos matemáticos) los generaremos mediante la función meshgrid()
z = data #datos que buscamos aproximar.
x = np.arange(24)
y = np.arange(24)
meshx,meshy = np.meshgrid(x,y)

Este algoritmo no es tan muy a la aproximación inicial, pero sí al conjunto de datos que buscamos ajustar. Si los residuales son muy altos o hay datos fuera de rango, ‘outers’, el algoritmo será inestable y seguramente lleve a una matriz singular que no podremos resolver directamente y deberemos resolverla por mínimos cuadrados. Esto último no será en la mayoría de los casos una opción viable.

El algoritmo toma pocas iteraciones para ajustar con una tolerancia muy baja. Un resultado para un subcaso bidimensional del problema anterior que empleé para validar el algoritmo:

fit

En el caso tridimensional el algoritmo funciona igualmente bien;

3adjust

Finalmente el mapa de errores, como se puede ser el ajuste es bastante bueno.

error_map

Conclusiones.

No estoy seguro de si las herramientas de las que dispone numpy pueden manejar los datos de la forma en que mi problema los dispone, la función numpy.optimize.lsqr() es versátil y seguramente pueda manejar arrays generados por la función meshgrid() si los pasamos propiamente aplanados (mediante ravel() o reshape() ).
En cualquier caso es un algoritmo fácil de implementar para problemas de optimización complejos y de varias variables en los que tenemos más o menos claro cuál es el tipo de función que va a ajustar bien con los datos.
Puesto que el factor sensible del algoritmo son los datos fuera de rango es recomendable preprocesar los mismos mediante alguna herramienta estadística (en breve pondré algo relacionado).
Publicado en matemáticas, python | Etiquetado , , , , , , , | 3 comentarios

Qué mejor despedida.

El libro ganado en la PUMPS summer school del Barcelona Supercomputing Center firmado por los autores.

IMG_20130731_203137 IMG_20130731_203216

Minientrada | Publicado el por | Deja un comentario

Here is our family picture for PUMPS 2013!

pumps_family_2013

PUMPS 2013 on BSC, Barcelona (Spain)

Imagen | Publicado el por | Etiquetado , , , | Deja un comentario

PUMPS 2013 at bsc, check in.

Just the beginning, almost a hundred people in the room!

PUMPS 2013

PUMPS 2013

Publicado en Sin categoría | Deja un comentario

Scripts, road to PUMPS 2013

Purpose.

I want to left here a couple of Python scripts for connecting to a remote server or synchronizing folders between host and server without using an “smart client”, all from terminal.

Connecting with the server.

This small script connect your pc via SSH to the server. Pexpect carries out the login process, you should install it -very easy step- before continuing.

import pexpect
import os as os

'''
 Login.py scrip.

 Change only #fields
 @ #server_direction: something@something.es
 @ #password : your password
'''

def screen_logo ():
 print " "
 print "--------------------------------------------------------------------"
 print " Login message... write something hacker' style here! "
 print "--------------------------------------------------------------------"
 print " "

screen_logo()
child = pexpect.spawn ('ssh #server_direction -p 2222')
#this will be the server' answer, so could be totally different!
child.expect ("#server_direction's password:")
print "... Login, please wait a few seconds."
child.sendline ("#password")
child.interact()

Synchronization part.

This one is sightly more complicated due to the huge range of options  that rsync offers, but the minimal script is the following;
import pexpect
import os as os

'''
 Sync.py scrip.

 Complete these fields bellow
 @ server_direction: something@something.es
 @ password : your password
 @ server_path : destiny path
'''

def screen_logo ():

 print " "

 print "----------------------------------------------------------------------"
 print " Synchronizing files "
 print "----------------------------------------------------------------------"

 print " "

def close_transfer():
 print " "

 print "----------------------------------------------------------------------"
 print " Files synchronized! "
 print "----------------------------------------------------------------------"

 print " "

screen_logo()

child = pexpect.spawn ('rsync -rvuh -e "ssh -p 2222"
 --stats --progress --exclude "/sync_folder/."
 --exclude "/#optional_argument/.."
 /#host_path/ #server_direction:/#path')

#this will be the server' answer, so could be totally different!
child.expect ("#server_direction's password:")

print "... synchronizing, please wait a few seconds."

child.sendline ("password")
child.interact()

close_transfer()

You probably don’t need more complicated options out of these ones (excluding folders and sync directories recursively).

Cuda snippets for Vim (cooking right now!).

I’m working on that. This evening I’ll have five free hours in the train so this will be finished. Tonight I’ll post for all old school code-writers or modern nostalgics. ,)
Answer me here as usual!
Publicado en Sin categoría | Deja un comentario