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