lunes, 10 de septiembre de 2012

Expression templates y lazy evaluation

Estos días estoy dedicandome a escribir mi programa (simulación de transporte de sedimentos por oleaje, para los interesados) usando técnicas de template metaprogramming. Más concretamente, los expression templates. Ambas cosas tienen fama de complicadas (y hay que reconocer que lo son), pero también son muy interesantes. En este artículo voy a intentar explicar un poco cómo funcionan y por qué merecen la pena.

El artículo me ha salido largo y un tanto árido; pero quiero insistir en que creo que merece la pena ;-)

Templates

Lo primero que voy a hacer es hablar un poco de los templates según c++. Como casi todo en c++, tienen una sintaxis compleja y su funcionamiento y reglas lo son más, de hecho hay muchas cosas que aún no entiendo.

Otros lenguajes han optado por sistemas más sencillos, que en este artículo voy a englobar con el nombre de generics.

El problema es que los generics sirven para lo que hace años, la mayoría pensábamos que servían los templates; es decir, para cosas como:

vector<double> v;

O dicho de otra forma: para especificar clases que dependían de otros tipos con la intención de optimizar el código. Esto es lo que se llama instanciación explícita de los templates.

Por otro lado, c++ hablaba de la instanciación implícita o automática. El tipo de cosas de c++ que uno suele pasar por alto, porque las explicaciones de lo que significa suelen tener la luminosidad poética del BOE. Afortunadamente, otros no lo pasaron por alto. La idea básica es que si uno tiene una función template del tipo:

template<class T> sqr(T val) { return val*val; }

Cuando el código hagamos sqr(a), el compilador obtendrá el tipo de a y instanciará la versión adecuada de sqr, sin que tengamos que decirle el tipo explícitamente. Las reglas sobre cómo se deduce el tipo y cómo se escoge la versión adecuada del template son complicadas y no voy a entrar en ellas (por que no me las sé y lo que hago es procurar evitar ambigüedades).

¿Cuál es la ventaja de esto? Es difícil de explicar. Para hacerlo voy a usar el ejemplo que conozco mejor, y para ello abro otro tema.

La sobrecarga de operadores: ayer y hoy

Una de las novedades que trajo el c++ al mundo de la programación (al menos para el gran público) fue la sobrecarga de operadores. Por fin uno podía definir una clase de números complejos y operar con ella como si fueran números reales. Lo mismo podía decirse de vectores, matrices, etc.

Pero había un problema, cada vez que se realiza una operación hay una llamada a una función y se devuelve un valor. Algo que no es muy importante en el caso de objetos pequeños como números complejos, pero que es un problema grande en el caso de clases de matrices y vectores (pensemos en vectores de cientos de miles de elementos). Supongamos la siguiente operación:

\begin{aligned} a\times(x+2\times(y+z)) \end{aligned},

Donde todas las variables son vectores y se supone que las operaciones ocurren elemento a elemento. En este caso, la suma y+z genera un nuevo vector, la multiplicación por 2 otro más, la suma con x otro y, finalmente, el producto por a, otro más. Ni siquiera es seguro que esos temporales no se destruyan hasta el final de la operación, con lo cual tendremos la memoria ocupada con temporales innecesarios.

Otro problema asociado es que este tipo de expresiones no son adecuadas para la optimización. En una expresión como la anterior pero para números reales, el compilador puede encontrar todo tipo de optimizaciones. Por ejemplo, el producto por dos puede ponerse como una suma del valor consigo mismo. En el caso de vectores, no es posible.

El efecto de todo esto es que en c++, si queríamos hacer un código "bonito" teníamos que renunciar a la optimización, y no hablamos de porcentajes pequeños. En las pruebas que he hecho he visto reducciones de velocidad de menos de un 50%.

Así estaba la situación hasta que a alguien se le ocurrieron los expression templates. La idea es incluir en nuestras clases de vectores un método double Eval(index). Este método lo único que hace es devolver el valor en la posición index del vector.

El siguiente paso es crear un operador template como este:

template<class L, class R>
ExprAdd<L,R> operator+(const L & op1, const R & op2)

{
  return
ExprAdd<L,R>(op1,op2);
}
Basicamente, lo que hace es devolver una clase que construye sobre la marcha basándose en los operadores. En el resto de la explicación voy a hablar sólo de sumas, pero es claramente extensible a otros operadores.

Claramente, este operador es aplicable a cualquier suma, tengan los operandos el tipo que tengan. Afortunadamente, si hay otros operadores definidos sin template, el compilador escogerá esas versiones.

La parte interesante consiste en que si hacemos algo como:

\begin{aligned} a+b+c+d+e \end{aligned},

El compilador coge (d+e) y crea una nueva clase ExprAdd<T,T> donde T es el tipo de las variables. (c+d+e) tendrá tipo: ExprAdd<T, ExprAdd<T,T> >...y así sucesivamente. La expresión final tendrá tipo:

ExprAdd<T, ExprAdd<T, ExprAdd<T, ExprAdd<T, ExprAdd<T,T>>>>>



Afortunadamente para nosotros, la instanciación implícita hace que no sea necesario que conozcamos el tipo de esa clase.

Lo que hemos hecho hasta ahora es construir clases, pero hasta ahora no se ha hecho ninguna operación. Para eso tenemos que ver cómo es la clase ExprAdd. Una versión sencilla sería algo como:

template<class L, class R>
struct ExprAdd
{
  const L &Lop;
  const R &Rop;

  Expr(L &ALop, R &ARop) : Lop(ALop), Rop(ARop) {};

  double Eval(int idx) {return Lop.Eval(idx)+Rop.Eval(idx); };
};

En pocas palabras, es una clase que guarda referencias a los operandos Lop y Rop y, lo más importante, tiene el mismo método Eval(index) que tenían los vectores. Este método llama a los Eval() de los operandos y suma el resultado.

El caso es que los operandos pueden ser vectores (con el Eval definido en la clase vector) o clases del tipo ExprAdd. La función Eval final de la expresión llamará a las funciones Eval de subexpresiones que a su vez llamarán a otras. Todas estas llamadas son conocidas en tiempo de compilación. Aquí viene otra de las ideas importantes, al ser llamadas conocidas en tiempo de compilación, el compilador las inlineará, de forma que, al final, la función Eval de la expresión contendrá un código similar a: a+b+c+d+e, pero donde las variables son ahora reales.

El detalle final es que los cálculos se realizan en el operador igual. Cuando tenemos una expresión del tipo: z = a+b+c+d+e, donde todas las variables son vectores del mismo tamaño,  en a+b+c+d+e se construye una clase que incluye un método Eval y en el operador = se lo llama y se realizan las operaciones, elemento a elemento.

Es decir, no se construyen vectores temporales y el compilador es capaz de hacer todo tipo de optimizaciones. Además, nos permite escribir código claro y bonito. Todos contentos.

Lazy evaluation

De hecho es aún mejor. Esta técnica nos permite una expresividad nueva, que con los métodos antiguos no hubiésemos usado. La técnica está basada en la lazy evaluation.

En pocas palabras, lazy evaluation significa que los cálculos no se hagan hasta que no hagan falta. Para ver cómo va, es mejor poner un ejemplo.

Supongamos que tenemos un vector v y que en ciertas partes del código necesitamos usar su cuadrado. Una posibilidad es hacer algo asi:

v2 = v*v;
a = v2+b;
c = v2-e;
etc...


Que requiere más memoria para v2, a lo que hay que sumar los consabidos problemas de rendimiento. Otra posibilidad es construir la clase:

struct Tsqr
{
  const Vector &V;

  Tsqr(Vector &AV) : V(AV) {};

  double Eval(int idx) { return V.Eval(idx)*V.Eval(idx); };
};

Y la función:

Tsqr sqr(Vector &v) { return Tsqr(v); }

La función sólo sirve para construir la clase Tsqr, que es la que realiza el trabajo. Puesto que el único requerimiento de una clase para poder ser usada, es que incluya el método Eval, la clase funcionará como cualquier otro vector. Ahora podemos usarla en expresiones como:

a = sqr(v)+b;
c = sqr(v)-e;
etc...

Donde, ahora no se guarda nada en memoria. Los cálculos del cuadrado de v se realizan de forma lazy, cuando son necesarios.

De hecho, usando la nueva instrucción auto  de c++11, podemos hacer cosas como:


auto X = a+(b*c);

Donde X es ahora la expressión a+(b*c), sin que se realice ningún cálculo en esa línea (auto permite definir variables deduciendo el tipo de la expresión tras el igual, lo cual añade un nuevo uso a la instanciación automática). Luego podemos usar X más adelante en otras expresiones y será entonces cuando se realizen los cálculos, de nuevo, de forma lazy.

Esto es todo. Como no he usado ninguna imagen hasta ahora, pondré una como premio a los que hayáis llegado hasta aquí.

Los 80 también tuvieron cosas buenas


SVG y asteroides (duros)



Hasta la semana pasada, para mí el SVG era un formato vectorial de gráficos que los navegadores de internet deberían aceptar (y que no todos lo hacían), además de ser el formato usado por Inkscape.

Mi opinión sobre él era que era el típico resultado de un diseño por comité: lo justo para contentar a todos pero insuficiente para gustar realmente a nadie.

El caso es que siempre me había llamado la atención el hecho de que Inkscape permite asignar eventos a los elementos del dibujo. Intentando averiguar qué era eso me encontré con que es posible hacer SVG dinámicos; aunque, cómo siempre, la posibilidad no está implementada en todos los navegadores.

Sin embargo me picó la curiosidad. Más que nada porque desconozco casi todo sobre la programación web, y cuando digo casi todo quiero decir que el otro día descubrí lo que era el DOM.

Para probarlo (y en el proceso, "probar" que soy un friki  que vivió en los 70), se me ocurrió hacer una versión del juego Asteroides más o menos fiel al original, con algún toque de color, pero sin convertirse en uno de esos engendros de remake que se ven por ahí.
Cuando digo "engendro" me refiero a algo así
La cosa resultó sorprendentemente fácil. En un par de horas tenía una nave girando y un asteroide. Después de tres días de baja dedicación ya tenía casi todo lo importante funcionando. El movimiento es bastante suave, incluso para tamaños grandes de pantalla. Eso sí, el código dejó de funcionarme en firefox cuando añadí la detección de choques entre disparos y asteroides; pero en Opera y Chrome sí funciona (aunque aún hay problemas con el zoom). Si queréis ver cómo va podéis probarlo aquí.

Y si no lo queréis probar, esta es una imagen
En este ejemplo, el SVG se abre directamente desde la URL, pero también es posible incrustarlo en una página HTML normal. Pienso que gran parte de las cosas que se hacen en Flash (salvo vídeo, que quizás sea, ahora, la más importante) pueden hacerse mejor en SVG.

Como bonus, aquí tenéis lo que, creo, es el diseño original de los asteroides que encontré en esta página. El que hayan puesto la raya en la O es lo que me convence de que es un documento original de la época, aunque ahora que lo pienso, la raya es precisamente para distinguir ceros de oes...


domingo, 9 de septiembre de 2012

El canal de Garonne

Viendo el puente siguiente, se me dió por pensar brevemente la siguiente pregunta:
¿Como calcular el peso máximo del barco/s que pueden pasar por el puente?.
Os lo dejo para pensar (poco pero algo).

lunes, 3 de septiembre de 2012

En el principio fue el código máquina (I)

Por petición popular (de una sóla persona) inicio aquí una serie de entradas sobre diferentes aspectos de la optimización de programas: ¿Qué cosas funcionan y cuáles no? ¿cuándo es necesario optimizar a mano y cuándo es mejor dejar que el compilador lo haga por tí? etc. Me da que no es este un tema popular, entre otras cosas porque la historia nos dice que, debido a los cambios en el hardware, lo que era una optimización ayer, ahora resulta que es más lento que un código menos trabajado. La sensación es que el esfuerzo es fútil. De hecho ese es uno de los temas que pienso tocar en el futuro: ¿Por qué el paralelismo y las jerarquías de memoria (registros->caches->RAM->Disco) no son meros caprichos de la tecnología actual, sino que hay leyes físicas fundamentales que garantizan que vamos a tenerlas en el futuro? El tema de hoy comenzó, precisamente, como una pregunta referente a la memoria: ¿cuándo es preferible guardar unos cálculos en memoria para no tener que rehacerlos? La respuesta es muy compleja y totalmente caso-dependiente. Existen reglas generales, pero cada caso concreto tiene sus diferentes problemas. Por lo tanto, vamos a considerar el caso concreto, relacionado con un programa de ajedrez:

const int UP = 8;

inline int Sign(int color)
{ return (1 - (color * 2)); }  // White sign is +1, Black sign is -1

inline int Ahead(int color)
{ return (Sign(color) * UP); }  // UP is -8, DOWN is +8


Supongamos que la función Ahead se llama muchas veces siendo color un valor que sólo admite los valores 0 y 1. ¿Es más rápido precalcular un array ahead de dos valores?

En este caso creo que la respuesta es clara. Los cálculos son extremadamente sencillos y pueden hacerse sobre registros, así que ni hablar de precálculos. Aunque para dos únicos valores es seguro que la cache se utilizaría eficientemente, acceder a registros es más rápido que acceder a cache.

Sin embargo aún quedan algunas preguntas por responder:
  • ¿Los inline hacen algo?
  • ¿Se conseguiría alguna ganancia optimizando a mano Ahead y sign?
Para responder a esas preguntas lo mejor es recurrir al código máquina.

La solución está aquí


Para ver el código que genera el compilador escribí el siguiente programa:

int main()
{
    int b=rand() & 1;  // Me aseguro de que b=0,1

    int a=Ahead(b) ;

    printf("Hello world %i !\n", a);
    return 0;
}


El rand() es necesario porque si ponemos una constante literal el compilador no se molesta en incluir las funciones y llama a printf con el valor calculado en tiempo de compilación. El printf es necesario porque si no usamos el valor de a, el compilador no se molesta en incluir nada.

El código máquina generado por gcc con la máxima optimización (-O3), pero sin más florituras es:

sub    rsp,0x8              
call   0x400470 <rand@plt>  // Llama a rand()

and    eax,0x1              // b = rand() & 1  

    
mov    esi,0x4006bc         // N.P.I.

mov    edi,0x1              // N.P.I.

neg    eax                  // eax = -b  
add    eax,eax              // eax = eax+eax = -2*b

lea    edx,[rax*8+0x8]      // edx = 8*eax+8 = 8 - 16*b

 
xor    eax,eax              // eax = 0
call   0x400460 <__printf_chk@plt> // Llama a printf

xor    eax,eax
add    rsp,0x8
ret    
                             


La parte del programa donde se hace la "llamada" a Ahead son simplemente tres instrucciones en código máquina que realizan la operación a = 8 - 16*b.

Se pueden hacer muchos comentarios respecto al código y la cantidad de trucos que usa:

1. Se da cuenta de que la llamada combinada a las dos funciones se puede simplificar como:

Ahead(color) = (1 - (color * 2)) * UP = 8 - 16*color

2. Aprovecha que teníamos color en eax y lo suma consigo mismo para obtener el doble de una forma rápida y sin necesidad de parámetros adicionales.

3. Usa la instrucción lea (load effective address) como un truco para obtener un cálculo aritmético. Esta función suele usarse para calcular punteros y moverse por un array (índice + offset), pero aquí se usa para hacer un cálculo. Desconozco si es posible escoger cualquier multiplicador o si hay sólo unos cuantos a elegir (lo que explicaría por qué no hace simplemente [rax*16+0x8]).
 
En resumen, reto a cualquiera a que escriba un código máquina más eficiente que el que genera el compilador para este caso.

Un detalle interesante es que, en este caso, si hubiéramos usado un array, los "cálculos" necesarios para acceder a los elementos son casi igual de complicados que el cálculo en sí.

Respecto a la otra pregunta: pues no, el gcc genera el mismo código. Da igual si hemos puesto la directiva inline o no.

(Continuará)

sábado, 25 de agosto de 2012

Un céntimo por tus pensamientos

El otro día leí una noticia que me llamó la atención: una técnica para ver lo que una persona está viendo reconstruyéndolo a partir de una resonancia magnética del cerebro. Concretamente el precórtex visual. En el siguiente vídeo lo explican:



Para los que no pilléis bien el inglés hablado, podéis saltar a 1:47, donde se muestran las imagenes. A la izquierda está el vídeo que el sujeto está viendo y a la derecha la reconstrucción.

La técnica usada parece sorprendentemente simple: muestran vídeos al sujeto durante unas horas mientras graban la actividad en el precórtex visual. Después realizan un análisis de correlación entre ambos. Finalmente muestran al sujeto otro vídeo (supuestamente diferente a los usados para obtener la correlación) y son capaces de reconstruir imágenes similares a las que el sujeto está viendo. El investigador explica en el vídeo que el sistema funciona mejor cuando lo que el sujeto está viendo son cosas simples, como caras y objetos fijos.

Esta investigación ha sido publicada en Nature, con lo cual no es probable que se trate de un hoax. Sin embargo, cuando uno se pone a mirar las reconstrucciones, asaltan las dudas.

La primera curiosidad es que aparecen textos en las imágenes:



Supongo que esto es debido a que los vídeos iniciales que se le mostraron durante el proceso de entrenamiento del sistema incluían texto; pero me da la impresión de que hay demasiadas. No sé que relevancia tiene esto, pero me llama la atención.

La segunda cosa que me llama la atención es que muchas imágenes parecen estar compuestas por las superposición de unas pocas imágenes:



No sé vosotros, pero yo ahí veo un plato en un lavavajillas, superpuesto a un par de imágenes más. Imagino que parte de lo que hacen es buscar las imágenes del proceso de entrenamiento que mejor correlacionan con las medidas obtenidas en la resonancia y combinarlas entre sí, quizás con un peso dado por lo bien que se correlacionan. No está mal, pero posiblemente esto haga que los resultados parezcan más espectaculares de lo que son. El caso anterior se beneficia de que es más común encuadrar objetos centrados en la imagen. Pero éste no es el caso más claro, los hay peores:


En este caso, está claro que la reconstrucción está basada en un 83% (siempre uso 83% cuando me invento porcentajes, lo cual ocurre menos de un 83% de las veces) en imágenes de un tío con camiseta negra. El hecho de que esté borroso es debido a que es la suma de varias imágenes en diferentes posiciones, y no a que la reconstrucción es imperfecta (bueno, aunque también puede verse así).

Es decir, las imágenes son tan buenas porque para la reconstrucción usan imágenes previas que ya tienen gran parte de los detalles necesarios.

En resumen, el sistema es interesante; pero probablemente es por el momento menos espectacular de lo que puede parecer por el vídeo. De todas formas, démosle unos años para mejorar, supongamos que podemos construir escáneres de resonancia magnética portátiles y se me ocurren un par de aplicaciones interesantes para este sistema (de las cuales sólo dos están relacionadas con el porno).

Para los que tengáis tiempo y ganas, podéis intentar experimentar con ello. Dicen que sólo dan los datos a los que investigan en neurociencia, pero desde una universidad casi seguro que es posible obtener los datos.

viernes, 10 de agosto de 2012

Python embebido

Iba a empezar esta entrada diciendo que "embebido" no significa lo que creemos que significa; pero me equivocaba. Significa lo que creemos que significa (es decir, no me equivocaba ¿o sí?)

Al parecer embebido no sólo significa estar absorbido (o absorto), sino que también tiene el significado de estar incrustado en algo. Es curioso porque "embedded" y "embebido" parecen tener orígenes bastante distintos, pero acaban significando cosas parecidas.

En fin, que voy a hablar de embeber Python en un programa, concretamente, un modelo numérico.

A diferencia de otro tipo de programas, los modelos numéricos suelen tener el problema de que en depuración es difícil obtener los valores que necesitamos. Si estamos resolviendo en una malla triangular, por ejemplo, necesitariamos hacer algún tipo de interpolación antes de poder ver los datos, y eso no es algo que el debugger pueda hacer.

La solución típica suele ser ir grabando resultados cada paso temporal y verlos desde otro programa. Los inconvenientes de esto son varios:

1. Generalmente sólo se graban ciertas variables importantes, así que no podemos acceder a ciertos valores.
2. Los archivos generados suelen ser muy grandes, pero poca de la informacion que se guarda se usa realmente.
3. A veces conviene ver los valores en momentos distintos al punto en el que se graba.

Hace tiempo tuve la idea de incluir un lenguaje interpretado junto con el programa, de forma que pudiésemos acceder a las variables en ciertos puntos de la ejecución. Bueno, la idea no es realmente mía, pero nunca la he visto implementada realmente.

Empecé a intentarlo con Python, que es un lenguaje bastante adecuado para ello (parece que cada vez se postula más como la alternativa de Matlab); pero la falta de documentación y ejemplos hizo que lo dejara... hasta ayer.

Ayer me puse con ganas renovadas y ya lo tengo funcionando. En pocas palabras, tengo un código que se compila con mi programa y que se enlaza con Python. En cada paso temporal tengo una llamada a una función que detecta si se ha pulsado una tecla y en caso afirmativo pone la simulación en pausa y permite escribir comandos Python. También es posible hacer que esto ocurra en un punto concreto del programa.

La siguiente figura muestra la línea de comandos, en la que se ve cómo accedo a ciertas variables (la presión y la componente X de la velocidad).




Y también como ejecuto un código externo (unas 6 líneas) que realizan el mallado e interpolación de la presión y la muestran en forma gráfica.

Es una gráfica aburrida, pero es la que tenía.


También es posible modificar las variables sobre la marcha, así que es posible hacer cambios en el programa sin recompilarlo. Hay muchas más posibilidades, aún estoy empezando a explorar el invento.

Supongo que es algo muy específico que no motivará a todo el mundo, pero a mí estas cosas me ponen.

A los interesados, puedo pasaros el código.

lunes, 23 de julio de 2012

El espagueti, ese misterio de la ciencia

Hace años leí que la costumbre que tienen los espaguetis secos de partirse en más de dos trozos cuando se doblan, había sido objeto de asombro y estudio por científicos del calibre de Richard Feynman. Yo, que soy de los que parten los espaguetis antes de meterlos en la olla (sí, ya sé que es un tremendo faux pas), había observado el fenómeno en muchas ocasiones, pero jamás se me ocurrió pensar en ello. Se vé que eso es lo que diferencia a los grandes científicos de los que simplemente escribimos un blog.

En YouTube hay varios vídeos de espaguetis rompiendose, ninguno de ellos especialmente bueno (al menos entre los que he encontrado). Aquí hay uno:

 

¿La explicación? No parece trivial, he visto algunas que no me convencen y otras que me convencen algo más. En este artículo se habla de todo esto y hay varias referencias.

Pero aún hay más. Vamos ahora a los espaguetis cocidos. Algo que todos hemos hecho (y algunos seguimos haciendo) es aspirar un espaguetti con la boca. ¿Pero cuántos se han preguntado cuál es la fuerza que hace que el espaguetti se introduzca? Yo no, pero me encontré la pregunta con su respuesta hace poco y es cierto que la respuesta no es sencilla.

La cuestión es: cuando aspiramos hacemos un vacío en la boca, lo que quiere decir que hay menos moléculas de aire golpeando la superficie del espagueti en el interior de la boca comparados con los que hay en la parte interior. Lo que quiere decir que en realidad, el espaguetti es empujado desde la parte exterior y no "tirado" desde dentro, como parece... ¿o sí?

Aquí tenéis una respuesta con discusión.

Está visto que la pasta oculta muchos misterios.
¿Estará aquí la llave a una nueva teoría unificada?