miércoles, 26 de diciembre de 2012

BREVES: Bombas autosimilares

Cuando pensamos en figuras geométricas autosimilares solemos tener en mente objetos fractales complicados. El copo de Koch que mostré hace unos días es un ejemplo claro, al igual que el archiconocido conjunto de Mandelbrot.

Para los incrédulos, esta es una parte autosimilar del conjunto de Mandelbrot. 
Sin embargo, no tenemos que irnos tan lejos, de hecho una línea recta también es autosimilar ante cambios de escala y cualquier polígono regular lo es ante rotaciones. La propiedad de autosimilaridad es muy utilizada en aparatos mecánicos, el ejemplo más sencillo es la rueda. El hecho de que sea autosimilar es la propiedad que la hace útil, al igual que los engranajes funcionan porque son autosimilares ante ciertas rotaciones o traslaciones.

Esta propiedad también se ha utilizado desde antiguo para mover agua, generalmente de un punto a otro más elevado. El ejemplo más conocido es el tornillo de Arquímedes que, como su nombre indica, no fue inventado por él:

En ocasiones también se usa para elevar bolas rojas

Una versión menos conocida de este sistema requiere una dimensión menos:

Sí, ya lo sé, esta no es exactamente autosimilar

Siguiendo con las espirales, también tenemos este ingenioso compresor que se suele utilizar en aparatos de aire acondicionado:
Se recomienda hacer un silbido periódico para una mejor experiencia.

Os dejo pensando en otros métodos para usar la autosimilitud o las espirales y su aplicación al movimiento de fluidos.


domingo, 23 de diciembre de 2012

BREVES: Una de arañas

Hoy voy a salirme de los temas habituales y voy a hablar de biología.

Me he encontrado con la siguiente noticia:

DESCUBIERTA UNA ARAÑA QUE CONSTRUYE SU PROPIO DOBLE


Para la foto han elegido una que tiene ocho patas, pero buscando he visto otras fotos con menos:

Ésta, sin ir más lejos, ha hecho una escultura de un hombre


A simple vista podría parecer la típica noticia del tipo: "naturaleza asombrosa", e incluso citan a expertos de la universidad de Cornell (que por lo que tengo entendido, no es mala) para confirmar la noticia.

Según dicen, la razón más probable para ese compartamiento es que sirva para protegerse de (o confundir a) depredadores despistados.

A mi me resulta extraño que una araña, cuyo método de caza es, precisamente, hacer telas invisibles para otros insectos y luego esconderse, se ponga luego a construir una réplica de sí misma. Me parece que eso desvirtúa la intención original de las telas.

Además, habría que creer que esas "esculturas" engañan a depredadores (pájaros, supongo) que en general tienen muy buena vista y un cerebro adaptado a reconocer este tipo de presas.

Y sin embargo, al parecer esas arañas existen y es cierto que construyen esas formaciones en la tela. Pero tengo la fuerte impresión de que el motivo tienen que ser otro.

¿Ideas?

sábado, 22 de diciembre de 2012

El motor de agua

Desde siempre he oído a gente hablar del motor de agua. Supuestamente se trataría de un motor que funciona con agua en lugar de gasolina y que las grandes petroleras están conspirando para que no se haga público y les quite el negocio.

Existen muchas razones lógicas para no creerse algo así. Por circunscribirnos a las razones físicas, debería estar claro para cualquiera que piense un poco en ello que algunas de las razones que hacen que el agua sea tan abundante y necesaria para la vida (estabilidad química), la hacen casi inútil como fuente de energía.

De hecho, si se pudiese sacar energía del agua, probablemente habría algún ser vivo que lo haría; pero hasta dónde yo sé, ninguno usa el agua más que para disolver otras sustancias. Y eso que hay organismos que se alimentan de casi cualquier fuente de energía.

Sin embargo, esto no quiere decir que no podamos usar el agua en un motor. De hecho, la mayoría de los coches usan agua en el circuito de refrigeración. Este circuito tiene una bomba que mueve el agua entre el motor, donde se genera calor, y el radiador, donde se traspasa ese calor al aire circundante. Así se ha hecho durante más o menos un siglo sin apenas cambios, pero ¿Es posible hacerlo mejor?

No tengo respuesta a esa pregunta, aunque sí una posibilidad de que la respuesta sea afirmativa. Os presento el motor de seis tiempos de Crower:

Los primeros cuatro tiempos son los mismos que los de un motor diesel convencional: admisión, compresión, explosión y escape.

En el quinto tiempo se introduce agua líquida en el piston ya comprimido y se cierra la válvula correspondiente. El agua líquida se transforma en vapor (robando calor al sistema en el proceso) y genera presión para ayudar a una nueva expansión del cilindro.

En el sexto tiempo el vapor de agua es expulsado hacia un condensador, donde el vapor vuelve a convertirse en agua.

En este caso, el agua no sólo sirve de refrigerante, sino que también produce potencia. Puesto que la expansión del vapor ayuda a mover el pistón, el vapor se enfría en el proceso con lo que, al parecer, no es necesario que el motor tenga un radiador.

Imagino que hay muchos problemas por resolver y quizás, en la práctica, el motor de Crower no sea mejor que un motor diesel convencional. Eso no quita que sea una idea interesante.

domingo, 16 de diciembre de 2012

BREVES: Se hace energía al andar

Últimamente parece que hay gente empeñada en sacar energía de todas partes. Mucha gente parece no darse cuenta de la enorme utilidad que tienen las leyes de conservación: permitirnos no tener en cuenta los detalles. No importa lo complicado que sea un aparato, de un sistema cerrado no vamos a sacar más energía de la que le metemos.

Mi casero está empeñado en montar un sistema para sacar energía mientras él pedalea en una bicicleta estática. Le he dicho que la energía que produzca (y más) la tiene que comer primero y que el precio de los alimentos es generalmente mayor que el de la gasolina (en cuanto a valor energético). Su respuesta es que ahora hay unos generadores nuevos que producen más electricidad.

Claro que quizás le venga bien ese ejercicio y se pueda contabilizar todo el proceso como win-win. De todas formas, lo que es seguro es que no va a sacar electricidad suficiente como para vendérsela a la compañía eléctrica, que es lo que él pretende.

EPO: ahora con más Kw/h

Otros quieren obtener (aunque robar es probablemente una palabra mejor) energía de los coches que pasan.

Y también hay varios grupos investigando cómo obtener energía de una persona caminando. De nuevo, hay que tener en cuenta que esta energía debe provenir de un mayor esfuerzo al caminar (a no ser que los humanos seamos muy ineficientes al andar, cosa que dudo muchísimo). Quizás existan nichos en los que algo así pueda ser útil, pero he hecho unos cálculos rápidos para una persona normal:


Energía que consume una persona al caminar 1600 m: 480 kj. (http://www.ncbi.nlm.nih.gov/pubmed/15570150)
Energía que contiene una sola pila AA alcalina: 10 kj (http://www.allaboutbatteries.com/Energy-tables.html)
Lo que anda una persona normal al día es, precisamente, 1600 m (http://walking.about.com/cs/pedometers/a/2000steps.htm).


Pero claro, una cosa es la energía que consumimos al caminar y otra la que realmente se necesita, parte de la energía se va a calentar el cuerpo, sudar y otras cosas relacionadas con no morirse. Además, de la energía usada al caminar sólo un pequeño porcentaje debería emplearse en la recarga, ya que sino nos cansaríamos en exceso. También hay que tener en cuenta las pérdidas en el proceso de generación de energía y recarga.

Creo que siendo optimistas, sólo un 1% de la energía consumida al caminar podría recuperarse, es decir, en un día normal sacamos suficiente para recargar sólo media pila AA.

Por otro lado, la densidad de energía de la gasolina es unas 30 veces mayor que la de una pila alcalina (http://en.wikipedia.org/wiki/Energy_density), así que la gasolina equivalente a lo que generamos en el proceso ocupa un volumen 30 veces menor que el de una pila AA, es decir, aproximadamente una gota.

¿Alguien cree que merece la pena la investigación, el desarrollo y, finalmente, el uso de un aparato así para ahorrarse una gota de gasolina al día?

Cabe la posibilidad de que mis cálculos sean incorrectos, así que os animo a revisarlos; aunque no creo que la cosa cambie por mucho.

sábado, 15 de diciembre de 2012

BREVES: Huellas "digitales"

En un artículo de ayer de El Mundo dice (si no lo editan posteriormente):

El asesino no aparece en su anuario de la escuela secundaria('Camera Shy', dice el anuario) y dejó pocas huellas digitales, ya que al parecer ni siquiera tenía perfil en Facebook. En 2006 su hermano terminó el instituto y se marchó a estudiar a la Universidad de Quinnipiac en Connecticut.

No parece que el autor esté haciendo un juego de palabras consciente. En ese caso, lo propio habría sido entrecomillar (como he hecho yo en el título para hacerlo más interesante) o poner en cursiva la palabra digital, o quizás, redactar el párrafo de otra forma. Pero no, la frase no da a entender nada de eso y mi opinión es que el periodista no se dio cuenta de que las huellas digitales solían ser otra cosa.

Siempre se puede contar con Google Images para ilustrar un artículo de forma poco original

En una época en la que las series policiacas han dejado de lado las huellas digitales por otros métodos más modernos como el reconocimiento de huesos, fluídos corporales y un conocimiento enciclopédico de todas las sustancias químicas y quién las manufactura, no es extraño que asociar "digital" con "dedos" resulte arcaico.

¿Puede ser éste el primer uso de un nuevo modismo? Aunque hay un nicho de significado que llenar, las huellas digitales tradicionales siguén usándose y me parece que la cosa llevaría a múltiples equívocos.

Además, tengo la impresión de que casi todos los neologismos españoles son adaptaciones de sus equivalentes yankis y, en este caso, los fingerprints ingleses no aceptan el juego de palabras.

De hecho, existe en inglés el término digital fingerprint, que define sistemas para marcar la autoría de obras digitales. ¿Cómo debemos traducirlo? ¿Huella digital digital? ¿Digital y tal?


viernes, 14 de diciembre de 2012

BREVES: ¿Es toda la energía, energía cinética?

Inauguro aquí las entradas "breves", que como su nombre indica serán más cortas (y menos pensadas) de lo habitual. También serán, espero, más frecuentes.

La primera es simplemente, esa pregunta. ¿Es toda la energía, energía cinética?

Desconozco la respuesta. Espero que aquellos que hayan estudiado el modelo estándar (si es que hay alguien) me respondan. Mi impresión es que cuando nos vamos a la teoría cuántica de campos, tres de las fuerzas fundamentales se reducen a intercambio de partículas virtuales y ese contexto no parece que tenga cabida la energía potencial. Evidentemente, otros tipos de energía como la calorífica, química, elástica, etc. pueden ponerse como combinación de fuerzas fundamentales y agitaciones (energía cinética).

Y aún en el caso de que esto fuera cierto, queda la energía gravitatoria, que salvo que alguien encuentre una teoría compatible con gravitones, no parece que pueda ponerse en forma de energía cinética.

Energía en estado puro


Voy a terminar con una cita del libro de física de Feynman:


Lo cual está muy bien. Hay que desconfíar de aquellos libros que eluden las preguntas difíciles para no reconocer que hay cosas que no se saben. Estaría muy bien que aunque no sepamos qué es la energía, al menos pudiéramos decir que sólo hay un tipo.


lunes, 10 de diciembre de 2012

Concurso: planetas fotorealistas

En esta entrada voy a plantear un concurso. Los premios serán el consabido café y un regalo sorpresa por definir (cuya característica más probable será que se pueda encontrar en una dollar store).

Antes de ir a los detalles, una breve introducción.

Siempre me ha parecido curioso, que siendo los planetas uno de los objetos aparentemente más sencillos de pintar con un ordenador, la mayoría de las imágenes 3D de planetas que he visto, no se parecen gran mucho a la realidad. Empezando por el programa Celestia, que lo hace (bastante bien) en tiempo real:


Y terminando por las más elaboradas imágenes generadas con programas de renderizado 3D. Estas requieren más tiempo de proceso y suelen incluir complicados efectos atmosféricos y varias texturas superpuestas).


Los resultados, aunque atractivos, se parecen poco a la realidad:


En general, parece que la tendencia suele ser hacer la atmósfera más grande de lo que es para que resulte más llamativa, hacer las nubes más transparentes para que se vean mejor los continentes, dar un poco de luz ambiente para mostrar la cara oculta y, algo que se ve mucho últimamente, hacer visible la luz artificial en las partes oscuras.

Dicho esto, vamos a por el concurso. Se trata de mejorar un programa base en C++ que yo os proporciono. Este programa incluye una rutina de renderizado sencilla, pero fácilmente extendible. No usa OpenGL ni requiere grandes conocimientos de C++ o programación gráfica. También incluye capacidades de movimiento de la cámara y tres cuerpos celestes (Sol, Tierra y Luna).

¿En qué pueden consistir las mejoras? he pensado en tres categorías:

1. Mejorar la rutina de pintado de forma que los planetas sean lo más realistas posibles.

Cuando digo "realistas" me refiero a que parezcan realistas. No es necesario que simulen exactamente la costa terrestre o los cráteres de la Luna. Deben ser lo suficientemente parecidos para que los planetas sean reconocibles a simple vista, pero dejo espacio por si alguien quiere meterse en cosas de generación algorítmica de texturas y cosas similares. Como este es un tema subjetivo, al final habrá una votación para decidir los resultados.

2. Implementar la física de una nave espacial.

En el código actual he incluído, de forma muy sencilla, el movimiento de la cámara/nave. Mejorar esta parte es bastante sencillo, pero hacerlo bien no. El código incluye una pequeña librería de operaciones entre vectores y matrices 3D, así que casi todo el trabajo es de modelado físico.

3. Aceleración del código.

El código, tal y como está, va a unos 12 FPS en una situación típica. Esto empeorará con rutinas de pintado más complicadas, pero también hay mucho espacio para la optimización. Inicialmente me concentrE en que el programa fuese fácil de entender, pero a partir de ahora voy a trabajar en la optimización y daré actualizaciones al que las quiera; pero cosas como usar la GPU no entran de momento en mis planes. También es posible mejorar los algoritmos en formas que no se me ocurran a mi.

Ejemplo del programa tal y como está. Para manejarlo usar cursores, control izquierda, barra espaciadora y enter.
Teniendo en cuenta las estadísticas del foro y lo ocupados que solemos estar, no me hago muchas ilusiones de que alguien quiera participar en el concurso, pero eso no me detiene. Como fecha tope voy a poner el 28F del 2013.

El proyecto se puede descargar desde aquí. Incluye un archivo de proyecto para code:blocks además de makefiles para Linux y Windows. Está pensado para ser compilado con gcc, pero no dudo de vuestras habilidades para convertirlo a cualquier otro compilador (o, incluso, otro lenguaje). Eso sí, es necesario tener la librería SFML 1 instalada (en Ubuntu se puede hacer desde el Software center).

Una breve descripción de las rutinas

El programa tiene un par de módulos auxiliares donde se definen las clases usadas por el programa: TPlanet, TCamera, TRay y TSolarSystem. Las tres primeras son más o menos evidentes, la cuarta es simplemente un agregado de planetas con (de momento) un sólo Sol.

La parte interesante del programa se encuentra en main.cpp, que es un programa de sólo 180 líneas. Consta de sólo dos funciones, la rutina de pintado:


void Render(sf::Image &img, TCamera const &Cam, TSolarSystem const &SolarSystem)
{
  const int W=img.GetWidth ();
  const int H=img.GetHeight();

  static const uint32_t clBlack = ColorToUint(sf::Color(0,0,0));

  uint32_t *ptr = (uint32_t *)img.GetPixelsPtr();
  for (int j=0;j<H;j++)
  {
    for (int i=0;i<W;i++)
    {
      // For every pixel (i,j) in the image...

      uint32_t col=clBlack, c;
      Vector3d n;
      double dist, lastdist=1e30;

      // Computes the ray associated to each pixel.
      const TRay ray=Cam.Ray(i,j);

      for (int k=0;k<SolarSystem.Count;k++)
      // For every planet...
      {
        const TPlanet &planet=*(SolarSystem.Planets[k]);

        // If planet k intersects ray...
        if (planet.Intersects(ray, n, dist))
        {
          c=planet.Color(n, SolarSystem.Sun);
          if (dist<lastdist)
          // If this intersection is closer than the last...
          {
            col=c;
            lastdist=dist;
          }
        }
      }

      // Deals with the Sun in a slightly different way.
      if (SolarSystem.Sun.Intersects(ray, n, dist))
      // If Sun intersects ray...
      {
        c=ColorToUint(sf::Color(255,255,240));
        if (dist<lastdist)
        // If this intersection is closer than the last...
        {
          col=c;
          lastdist=dist;
        }
      }

      // Sets the collor of pixel (i,j).
      *ptr=col;

      ptr++; // Next pixel.

    }
  }
}


Esta función hace algo muy sencillo: Existe una cámara asociada a la pantalla y por cada píxel de esa imagen crea un rayo (TRay) de luz y comprueba si intersecta algún planeta. En el caso de que haya varios planetas que intersecten el mismo rayo, se queda con el más cercano y obtiene el color de ese píxel. Si el rayo no intersecta ningún planeta el píxel será negro.

La clase TPlanet implementa la rutina Intersects(ray, n, dist) que devuelve la distancia al punto de intersección más cercano a la cámara y el vector normal a la superficie del planeta en ese punto. TPlanet también implementa una versión básica de la función Color que devuelve el color del píxel en el punto de intersección.

El siguiente gráfico muestra la situación de forma más clara:


La parte del manejo de la nave también se encuentra en main.cpp y es la siguiente:

    // Process events
  sf::Event Event;
  while (App.GetEvent(Event))
  {
    rot = Vector3d(0.0,0.0,0.0);

    if (Event.Type == sf::Event::KeyPressed)
    {
       if (Event.Key.Code==277) acc=-10000000.0;    // Acceleration
       if (Event.Key.Code==257) acc= 10000000.0;    // Deceleration
       if (Event.Key.Code==293) rot += 0.02*Cam.xv; // Up
       if (Event.Key.Code==294) rot -= 0.02*Cam.xv; // Down
       if (Event.Key.Code==291) rot +=  0.2*Cam.zv; // Left
       if (Event.Key.Code==292) rot -=  0.2*Cam.zv; // Right
       if (Event.Key.Code==278) vel = Vector3d(0.0,0.0,0.0); // Stops movement.
    }

    if (Event.Type == sf::Event::KeyReleased)
    {
       if (Event.Key.Code==277) acc= 0.0;
       if (Event.Key.Code==257) acc= 0.0;
    }

    // Close window : exit
    if (Event.Type == sf::Event::Closed) App.Close();
  }

  // Camera/ship movement.
  vel=vel+(acc*dt)*Cam.zv;
  Cam.Pos+=dt*vel;
  if (norm(rot)>0) Cam.Rotate(unitary(rot), dt*norm(rot));


Que como se ve, no tiene mucha complicación y es bastante cutre.

miércoles, 5 de diciembre de 2012

El número de Reynolds (III): Copos de Koch y superfluidez

En esta tercera entrega sobre el número de Reynolds voy a hablar de los efectos de escala.

Antes de empezar, echad un vistazo a la siguiente animación:

Mira a la imagen fijamente durante un minuto...

Ahora mandame todo tu dinero.


Como muchos habréis reconocido, se trata de un zoom del famoso copo de nieve de Koch. La idea del zoom no es mía, si buscáis veréis animaciones similares a esta en otros sitios. Mi aportación en este caso se limita a ponerle motion blur (y, accidentalmente, ese movimiento a trompicones).

Ahora imaginad que ese copo de Koch representa una costa. En cada una de las infinitas bahías que van apareciendo podemos imaginar un remolino girando en su interior. Remolinos cada vez más pequeños según vamos acercándonos.

¿Qué es lo que hace que esta situación sea irreal?

Bueno, en primer lugar está el hecho de que la materia no es continua. Llegado un momento nos toparemos con los átomos y la situación cambiará completamente. Una situación así sólo puede darse en el caso simplificado de que supongamos que la materia es continua.

El caso es que esa es, precisamente, una de las primeras simplificaciones que hacemos al derivar las ecuaciones de Navier-Stokes, que los fluídos son continuos: las ecuaciones nos permiten, en principio, realizar zooms indefinidos.

...Y sin embargo, las ecuaciones de Navier-Stokes predicen que los remolinos irán desapareciendo a medida que nos vamos acercando.

De nuevo, las ecuaciones de Navier-Stokes

Aquí están, sin más preámbulos, en su forma no dimensional:



Si definimos la vorticidad como el rotacional de la velocidad, \(\boldsymbol\omega = \nabla\times\boldsymbol u  \). Suponiendo que el fluído sea incompresible y la densidad constante podemos tomar el rotacional de la ecuación anterior y obtener:

\begin{aligned} \frac{\partial\boldsymbol\omega}{\partial t} =- \boldsymbol u\cdot\nabla\boldsymbol\omega +\frac{1}{Re}\nabla({\boldsymbol\omega}^2) \end{aligned}

Esta es la llamada ecuación de la vorticidad. La razón para usar ésta en lugar de la anterior es que en este caso las únicas variables que tenemos son la velocidad o sus derivadas (la vorticidad), ya que la gravedad y la presión se nos han ido al tomar rotacionales.

La ecuación está en forma no dimensional, lo cual es una forma un tanto extraña de decir que está en unas unidades que no son fijas, sino dependientes del problema. En un caso como el del zoom al copo de Koch, las únicas varas de medir longitudes que tenemos deben estar relacionadas con las dimensiones de la imagen. Es importante darse cuenta de que no tenemos forma de saber en qué punto del zoom estamos, así que el tamaño de la imagen en un instante dado es realmente, la única cosa que podemos usar para medir.

Echando un vistazo a las ecuaciones, queda claro que casi todos los términos son invariantes bajo cualquier elección de unidades; la única cosa que puede variar es el número de Reynolds. Echémosle un vistazo:

\begin{aligned} \frac{1}{Re} =\frac{\nu}{UL} \end{aligned}

U suele definirse como una velocidad característica del flujo; pero eso no es una propiedad de las ecuaciones, sino de una solución de éstas. En lugar de eso, voy a definir U=L/T, donde L está relacionada con el tamaño de la imagen en un momento dado (pongamos que es su diagonal) y T es un tiempo característico por definir.

Como ya he dicho, para que la ecuación de la vorticidad sea invariante bajo cambios de escala tiene que darse que, según bajemos en el nivel de zoom, el número de Reynolds, se mantenga constante. Para que esto ocurra, es preciso que:

\begin{aligned} \frac{1}{Re} =\frac{\nu}{\frac{L^2}{T}} = const ~~\Rightarrow ~~ T=kL^2 \end{aligned}

Donde k es una constante. Es decir, las ecuaciones se mantienen invariantes ante cambios de escala de la forma:

\begin{aligned}x'=\lambda x ~~,~ t'=\lambda^2 t \end{aligned}

Esto es interesante, especialmente teniendo en cuenta que empecé este artículo con la intención de mostrar que las ecuaciones de Navier-Stokes no son invariantes ante cambios de escala y resulta que sí lo son, aunque no sea un cambio de escala lineal.

Dicho de otra forma, lo que quiere decir todo esto es que usando el cambio de escala propuesto, y usando las mismas condiciones iniciales y de contorno, obtendremos la "misma" solución.

He entrecomillado el "misma" porque la solución es numéricamente la misma, pero expresa resultados en distintas unidades. Supongamos que tenemos dos velocidades v numéricamente iguales, pero una para un dominio de longitud L y otra para un dominio de la mitad de longitud.

Si comparamos ambas velocidades, tendremos:

\begin{aligned} \frac{vLT^{-1}}{v0.5L (0.25T)^-1} = \frac{1}{2} \end{aligned}

Es decir, si usáramos las mismas unidades para medir ambas velocidades resultaría que aquellas medidas en el dominio pequeño serían el doble de grandes que las medidas en el dominio grande.

Para el caso del copo de Koch, esto significa que es imposible encontrar una solución real de flujo que sea autosimilar. Cualquier vórtice en una de las "bahías", tendría velocidades cada vez mayores en bahías cada vez más pequeñas, en una ascensión exponencial.

Lo anterior es válido para viscosidad \(\nu\) constante. Es interesante observar lo que ocurre para la viscosidad equivalente debida a la turbulencia, \(\nu_T\). En este caso, usando argumentos dimensionales tenemos que la viscosidad de remolino escala como:

\begin{aligned}  \nu_T = \alpha \frac{L^2}{T}  \end{aligned}

Es decir, que la viscosidad de remolino sí que escala de la forma adecuada, es únicamente la viscosidad molecular (cuya longitud característica está asociada a la agitación térmica) la que nos impide tener una solución autosimilar del flujo.

Afortunadamente, existe un fluído sin viscosidad. En ciertas condiciones el helio muestra superfluidez, así que la posibilidad de una solución autosimilar existe, al menos hasta llegar al nivel de los átomos.



martes, 13 de noviembre de 2012

Juegos de palabras

Al principio fue el verbo.

Hay tantas interpretaciones del significado del comienzo del evangelio de Juan como teólogos (es decir, pocas). La más habitual es suponer que "el verbo" (o "the word" o "logos", según el idioma) se refiere a Jesús.

Otra interpretación, que es la que yo prefiero, es suponer que el verbo es exactamente eso, una orden verbal. El pensamiento mágico asigna poder real a las palabras. En ese contexto, Dios sería el über mago, creando el universo con un simple "hágase", en un big-bang abracadabrante.

Ni que decir tiene que yo no creo en el poder mágico de las palabras, pero sí me parece que tienen mucha importancia en el modo en el que nuestro cerebro es capaz de razonar. El hecho mismo de asignar una palabra a un concepto nos permite pensar sobre ello con más facilidad. La mente (al menos la mía) necesita palabras para representar internamente categorías.

Casualmente, ayer leí la siguiente frase de Gödel:
The human mind is incapable of formulating (or mechanizing) all its mathematical intuitions, i.e., if it has succeeded in formulating some of them, this very fact yields new intuitive knowledge [...]
Dicho de otra forma, cada nueva formulación matemática nos sirve como un clavo al que aferrarnos para generar nuevo conocimiento y nuevas intuiciones. De forma similar, me parece a mi, nuevas palabras cubriendo nuevos conceptos nos permiten generar nuevas intuiciones y aumentar el espacio de lo cognoscible. De ser así, esto podría ser una explicación del effecto Flynn.

No sólo las palabras, la estructura del lenguaje en sí podría tener efecto en la forma en que pensamos y nuestra capacidad intelectual. Hay incluso lenguas artificiales como Ithkuil, pensadas específicamente para permitir pensar sobre cosas más profundas de lo que las lenguas naturales permiten.

¿Y todo este rollo para qué? os preguntaréis. La verdad es que quería hablar de la herramienta de n-gramas de Google y el prólogo se me fue de las manos.




Pero vamos a ello.

Google ngrams

Esta página de Google permite estudiar la frecuencia de uso de palabras en inglés. La fuente son los libros escaneados por Google y, por tanto, hay que tener en cuenta que no es perfecta. En particular tenemos los siguientes problemas:
  • El OCR no es perfecto y puede haber errores en la interpretación de las palabras.
  • El número de obras escaneadas decrece cuanto más nos movemos atrás en el tiempo. Puesto que los resultados aparecen en forma de porcentaje, esto no resulta claramente visible, pero hay que tener en cuenta que tales porcentajes están más sujetos a error cuanto más antiguos son. Para evitar esto no voy a incluir resultados muy antiguos, donde el ruido es más claro.
  • Hay palabras que pueden haber tenido usos en cierta época que no son obvios actualmente: una marca comercial como General Electric puede aumentar artificialmente el uso de "electric" y "general" en determinada época, lo mismo es aplicable a apellidos.
Teniendo en cuenta estas cosas, vamos a divertirnos echando un vistazo a ciertas gráficas (no soy el primero y quizás hasta me repita).



Con este gráfico intentaba ver cómo habían cambiado los temas en la ciencia ficción (y junto a ellos, en la población en general). Al principio, las historias hablaban de cohetes, más tarde de robots y últimamente de realidad virtual. Hay un hueco alrededor del 1975 que yo identifico con la ciencia ficción sobre sociedades humanas (distopías y similares). Como curiosidad, la palabra "society" tiene un pico claro sobre esos años, aunque al ser una palabra común se nos sale de la escala.


No hay una fina línea entre "love" y "hate"; hay un salto considerable, piensen lo que piensen los pesimistas. Además ¿quién habría dicho que los años 60 no sólo están cerca de ser el minimo absoluto del uso de "love", sino que además la tendencia es descendiente?

Al parecer el porno se inventó más o menos a la vez que internet. Entre el 1970 y el 2000 parece que el porno es la derivada de internet (yo pensaba que sería al revés).

Hay muchas más posibilidades y no quiero privaros de la diversión de buscarlas vosotros. Sólo voy a incluir una más que no sé cómo interpretar.


Me parece curioso que una palabra tan común como "table" cambie de frecuencia de uso tan clara y consistentemente. Sus dos significados principales son "mesa" y "tabla" (en su sentido gráfico, no de madera). Ninguno de ellos me parece que debería variar significativamente, pero aquí vemos cambios importantes. ¿A alguien se le ocurre alguna razón? ¿Por que los americanos dejaron de hablar de mesas o tablas en los 70?

Pongo "table" como ejemplo, pero prácticamente todas las palabras muestran cambios similares. Es realmente difícil encontrar una que no esté sujeta a ellos. La más constante que he visto es "in", y eso que parece que su uso ha empezado a decaer en los últimos años (y sí,coincide con un aumento de "out", aunque ésta última tiene mucha más variabilidad).

miércoles, 26 de septiembre de 2012

El número de Reynolds (II): difusión superlumínica

En la entrada anterior hablé de los procesos de advección y difusión de un contaminante, la diferente importancia de cada uno y de cómo el número de Péclet (Reynolds) muestra la importancia de cada uno.

Al final había comenzado a hablar de cómo las ecuaciones de Navier-Stokes (NS) no son más que ecuaciones de conservación del momento con sus términos de advección y difusión, cuya importancia relativa está dada por el número de Reynolds.

En esta entrada voy a volver un momento atrás y comentar algo que se me pasó por alto en la anterior. No es que no me acordase de contarlo sino que no se me ocurrió hasta tiempo después.

Difusión superlumínica

Voy a empezar con un gráfico similar al de la entrada anterior, pero en esta ocasión he añadido una línea que representa un punto que se aleja del punto de interés a la velocidad de la luz.


Un físico no se encuentra cómodo hasta que no mete la velocidad de la luz en algún sitio.

Como dije anteriormente, para cualquier velocidad, existe un radio en las cercanías de cualquier punto, dentro del cual la difusión es la forma más rápida de "transportar información". Evidentemente, esto no puede ser correcto físicamente, es sólo un resultado de nuestro modelo de la difusión/viscosidad.

Esto también explica algo que afirmé en la anterior entrada, pero que no me molesté en argumentar: que en escalas pequeñas la difusión deja de ser una aproximación válida.

¿Cómo de pequeñas? Veámoslo. Vamos a suponer que hablamos de agua y de viscosidad. ¿A que distancia máxima la difusión es más rápida que la velocidad de la luz?

Una forma de hacerlo es calcular la distancia L que hace que el número de Reynolds sea 1 para la velocidad de la luz:

\begin{aligned} Re = 1 = \frac {Lc}{\nu} \end{aligned}

Usando valores estándar:

\begin{aligned} L = \frac {\nu}{c} = \frac {10^{-6} (m^2/s)}{3~10^5 m/s}\approx 10^{-11} m \end{aligned}

Es decir, más o menos la décima parte de un átomo de Helio. Parece que esto no nos va a dar problemas, al menos para el caso del agua. La brea tiene una viscosidad que es casualmente \(10^{11}\) veces mayor que la del agua, así que en ese caso, el radio es del orden de un metro. ¿Cómo tendrán en cuenta eso en las ecuaciones?

Me habrán dado un Ig-Nobel por parsimonioso, pero la difusión va que se las pela

En realidad, la difusión no puede ser más rápida que el proceso físico que la crea, que no es otro que la agitación térmica de las moléculas del agua. Nos quedamos sin ver radiación Cherenkov asociada a la difusión.

Sin embargo, la razón por la que decidí escribir esta entrada no es sólo porque me pareció curioso. Hay un motivo más fundamental.

Una de las cosas que he oído repetir como un mantra desde que llegué a Cornell es: "una vez que escojas un modelo de la realidad, olvídate de la realidad". La idea es que una vez que has hecho una serie de hipótesis simplificadoras no es lógico usar argumentos que caen fuera de esas simplificaciones.

En este caso, una de las hipótesis simplificadoras es la hipótesis del continuo: todas las magnitudes consideradas son continuas a todos los efectos. Por lo tanto, debemos olvidarnos de moléculas, agitación térmica y demás.

Es decir, la hipótesis del continuo implica la transferencia de información a velocidades mayores que la de la luz. Es parte de nuestro modelo y tenemos que vivir con ello. Al fin y al cabo, el suponer que un fluído es incompresible también tiene la misma implicación, y nadie se queja. Qué le vamos a hacer.







jueves, 20 de septiembre de 2012

El número de Reynolds (I): Introducción.

Aunque tengo por ahí empezada una serie sobre las ecuaciones de Navier-Stokes, quiero hablar un poco sobre el número de Reynolds sin tener que entrar a fondo en las ecuaciones, así que empiezo una serie nueva.

En esta primera entrada voy a hablar de conceptos básicos para luego ir entrando en materia. De hecho, antes de hablar del número de Reynolds en sí, voy a tratar un tema similar pero más sencillo.

Difusión y advección de sustancias


Supongamos que tenemos dos situaciones diferentes, en una tenemos agua inmovil, en otra agua moviéndose a velocidad constante \(u\). En el instante t=0 ponemos una gota de sustancia en un punto y esperamos a que pase un cierto tiempo. Supongamos además que la sustancia del primer caso tiene un coeficiente de difusión en el agua de \(\nu\), mientras que la del segundo tiene un coeficiente nulo (no se difunde).

Para los físicos, advección = convección. Para los ingenieros, sustancia = contaminante. Para todos 2D=3D.

La figura muesta el comportamiento de la sustancia en ambos casos. En el primero se diluye debido a la agitación térmica de las moléculas de agua, en el segundo avanza sin cambiar de forma a la misma velocidad que el agua. La ecuación que rige ambos casos (aunque no es necesaria para entender lo que sigue) es:

\begin{aligned} \frac {\partial c}{\partial t} + \nabla \cdot (c\mathbb{u}) = {\nu}\nabla^2c \end{aligned}

Una forma de comparar ambas situaciones es estudiar cuanto avanza la sustancia en un tiempo \(t\) desde el punto en el que fue depositada. En el caso de la difusión es necesario hacer una suposición adicional, así que definiré el avance como la distancia a la que la concentración es igual a la mitad de la concentración en el centro.

En el caso de la advección, claramente la sustancia avanza a una velocidad  \(u\) o, en términos de la distancia recorrida \(L_{adv}=ut\).

No voy a entrar en detalles de cómo se calcula, pero en el caso de la difusión la distancia recorrida es \(L_{dif}=\sqrt{4\nu t}\). Intuitivamente, es fácil ver que la difusión es menor cuanto más aumenta el tamaño de la mancha y disminuyen los gradientes, así que una dependencia con la raíz cuadrada del tiempo parece razonable.

La siguiente gráfica muestra la comparación entre ambos avances:



Como puede verse, una consecuencia interesante de que la difusión avance con la raíz cuadrada del tiempo es que existe un radio \(L_0(u,\nu)\) en el cual la difusión es más rápida que la advección. Esto es cierto independientemente de que velocidades o coeficientes de difusión tengamos; siempre habrá un valor de \(L_0\) dentro del cual la difusión avanza más rápido que la advección. Si pensamos en términos de información, todo punto dentro de ese radio recibirá noticia de la sustancia depositada más rapido debido a la difusión que a la advección.

Para tener en cuenta este efecto, se suele usar el número de Peclet;

\begin{aligned} Pe_L=\frac {LU} {\nu} \end{aligned}

Que no es más que un número de Reynolds aplicado a la difusión de sustancias en lugar de la difusión del momento (viscosidad). Como yo he usado la misma letra \(\nu\) que se suele usar para la viscosidad, ambos números tendrán la misma expresión.

También lo he hecho porque nunca me ha gustado dar nombres diferentes a cosas similares (¿series de Taylor y McLaurin?). Pienso que no hace más que confundir, el número de Peclet no es más que una versión del número de Reynolds.

Siguiendo con lo que estábamos, en la expresión anterior \(L\) es una longitud, digamos que entre el punto donde se ha vertido una sustancia y un punto de interés, \(U\) una velocidad característica del flujo de agua y \(\nu\), como ya he dicho, el coeficiente de difusión. Para entenderlo mejor podemos usar las expresiones anteriores para calcular el tiempo que tarda en avanzar la concentración desde el punto inicial a un punto a una distancia \(L\):

\begin{aligned} t_1 = \frac {L}{U},~t_2 = \frac{L^2}{4\nu} \end{aligned}

Dividiendo \(t_2\) entre \(t_1\):

\begin{aligned} \frac {t_2}{t_1} = \frac {LU}{4\nu} \end{aligned}
Es decir, el número de Peclet es una relación entre los tiempos que tarda en llegar la información. Si es mucho menor que uno significa que en nuestro problema, la difusión es más importante que la advección y su es mucho mayor que uno, todo lo contrario.

Si utilizamos \(L\) y \(U\) para adimensionalizar la ecuación de advección difusión obtendremos:

\begin{aligned} \frac {\partial c}{\partial t} + \nabla \cdot (c\mathbb{u}) = \frac{1}{Pe_L}\nabla^2c \end{aligned}

Es decir, el número de Peclet indica la importancia del término difusivo respecto al resto de términos.

Antes de pasar a otra cosa, quiero hacer un último comentario. La definición del número de Peclet (y del número de Reynolds) se basa en una cantidad objetiva \(\nu\), otra más o menos objetiva \(U\) y finalmente una que es subjetiva (en el sentido de que tenemos libertad para elegirla) \(L\). Esto tiene utilidad desde el punto de vista práctico, ya que en la realidad solemos estar interesados en calcular la difusión y advección entre el punto de vertido y la zona de interés.

Pero a nivel de procesos físicos, no existe una \(L\) única. La sustancia está, en cada punto y a cada instante, sujeta a procesos de advección y difusión locales, sin que la distancia \(L\) tenga ningún significado físico.

Más aún, podríamos argumentar que en cada punto, la difusión es infinitamente más importante que la advección (ver la gráfica anterior). ¿Cómo es posible que exista advección si localmente la difusión es infinitamente más importante?

Parte de la respuesta es que la difusión es siempre una aproximación (no es válida en un entorno infinitesimal), otra parte es que la ecuación es lineal (así que podemos ver ambos procesos como independientes) y luego está la parte de la respuesta que desconozco.

Pero no hablemos de cosas tristes. Quedémonos con la idea de que en cada punto existe un radio \(L_0\) dentro del cual la difusión es más importante que la advección.

Las ecuaciones de Navier-Stokes

La ecuación de conservación del momento, en forma no dimensional, puede escribirse así:

\begin{aligned} \frac {\partial \mathbb{u}}{\partial t} + \nabla \cdot(\mathbb{u}\mathbb{u}^T) = -\frac{1}{\rho}\nabla p +  \frac{1}{Re} \nabla^2(\mathbb{u}) + \mathbb{g} \end{aligned}

Donde \(Re=\frac {LU} {\nu}\)

Como es sabido, esta ecuación no es más que una ecuación de advección-difusión para el momento, más un par de términos adicionales: presión y gravedad. En esta ocasión, el número de Reynolds sustituye al número de Peclet, que como ya había dicho, son más o menos la misma cosa.

Continuara...








viernes, 14 de septiembre de 2012

En el principio fue el código máquina (II) : números y constantes

Como ya dije en un artículo anterior, me encuentro actualizando un código para incluir en él transporte de sedimentos. Se trata de un programa que originalmente "sólo" era capaz de simular multiples fluídos inmiscibles, incluyendo el efecto de la temperatura, electromagnetismo y otras cosas que nunca he tenido tiempo de mirar con detalle. Afortunadamente, el código está bastante bien compartimentado.
El programa inicialmente se desarrolló en Los Álamos. Un  insituto especializado en modelado numérico y que a la mayoría os sonará porque se creó para desarrollar la bomba atómica. Si en algún sitio deberían tener buenos programadores -pensaba yo- sería allí.

Stanisław Ulam, Richard Feynman y John Von Neumann en Los Alamos
Sin embargo, la realidad es que el programa es un caos, como todos los modelos numéricos que he tenido la desgracia de tener que modificar. Hay que decir que el diseño inicial no es del todo malo (aunque claramente no pensaban en eficiencia cuando lo hicieron); pero el resultado final es horrible: Trozos de código repetidos, funciones que hacen lo mismo que otras, sin explicación de por qué. Códigos que cláramente no hacen nada o hacen cosas erróneas... La lista sigue. El hecho de que el programa funcione aún así, da idea de la cantidad de redundancia que tiene.

En cualquier caso, no empecé este artículo para quejarme del programa en general (tal vez en otro artículo), sino para comentar un detalle que me parece curioso. Echemos un vistazo al módulo de constantes:
! Real whole numbers
  real(real_kind), parameter :: zero          =   0.0
  real(real_kind), parameter :: one           =   1.0
  real(real_kind), parameter :: two           =   2.0
[...]
  ! Real exponential numbers
  real(real_kind), parameter :: ten_tominus14 = 1.0e-14
  real(real_kind), parameter :: ten_tominus12 = 1.0e-12

[...]
  ! Real fractions
  real(real_kind), parameter :: one_168th     = one/168.0
  real(real_kind), parameter :: one_105th     = one/105.0
[...]
  ! Real constants
  real(real_kind), parameter :: pi  = 3.1415926535

Los [...] son míos, el código define todas las constantes que nos podamos imaginar. La última que he puesto, es la única que parece tener algo de sentido. La razón por la que la puse es porque me recordó el famoso comentario en un manual de Fortran de Xerox:

The primary purpose of the DATA statement is to give names to constants; instead of referring to pi as 3.141592653589793 at every appearance, the variable PI can be given that value with a DATA statement and used instead of the longer form of the constant. This also simplifies modifying the program, should the value of pi change.
Supongo que los programadores de Los Alamos, aún más cautelosos, tuvieron en cuenta la posibilidad de que números cómo 1 ó 2 cambien un día de estos. Nunca se es lo suficientemente precavido.

Lo cierto es que no se me ocurre ninguna razón válida para hacer algo así. Quizás la razón sea que así se aseguran de que one es un double (real(real_kind)) y no un float (real(4)), aunque hay pocas situaciones en las que eso sea un problema y siempre es posible escribir: 1.0d0.

En cualquier caso, una cosa me llevó a pensar otra y de ahí a una idea para un artículo. ¿Cómo manejan los compiladores las constantes reales?

A diferencia de las constantes enteras, que se pueden incorporar directamente en el código máquina; las reales se suelen guardar como datos (aunque en el segmento de código). Es decir, cada vez que en nuestro código hacemos algo como:


x=2.3*y;

El compilador crea una entrada para esa constante. Si le echamos un vistazo al ensamblador que genera el compilador, veremos una lista con todas los números (sin repetir) usados en el programa. Lo mismo pasa con los parameters listados antes.

Es decir, los números reales no se incluyen en el programa cada vez que usan (como sucede con los enteros), sino se guardan en una única dirección de memoria... memoria que se puede modificar.

Dicho de otra forma, es posible escribir código que modifique todos los 1.0 de nuestro programa (o todos los ones) y los cambie por otro número.

Aprovecho el impacto que la frase anterior os habrá provocado, para cambiar momentáneamente de tema.

Parte de lo que he hecho con el programa es escribir todo el código nuevo, y reescribir parte del viejo, en c++. Esto tiene sus ventajas e inconvenientes. Las ventajas las voy a dar por sabidas, las desventajas están relacionadas con la interacción entre ambos lenguajes, ya que no existe un estándar. Lo que significa que hay que tener en cuenta diferentes compiladores.

En el caso de gcc+gfortran, las funciones Fortran siempre pasan los parámetros por referencia. Incluso cuando el parámetro es un número, la función pasa el puntero a la dirección en memoria donde ese número está guardado. Si uno no se sale de Fortran no pasa nada, pero cuando una función c++ recibe el parámetro, recibe el puntero.

Afortunadamente, como dije antes, esos datos están en el segmento de código, que es de sólo lectura, así que normalmente el mayor peligro que tenemos es que obtengamos un segmentation fault que no sabemos de donde viene.

Pero claro, nadie dijo que no se podían cambiar las propiedades del segmento de código. Eso me trae de vuelta al tema original. Veamos el siguiente pedazo de código:

program hello
    interface
      subroutine func(d)
          real(8), intent(in)   :: d
      end subroutine func
    end interface


    print *, '2 + 2 = ',4.0D0
    call func(4.0D0);
    print *, '2 + 2 = ',4.0D0

end program

El código en rojo define la interfaz a una función en c++ hecha por mi. En particular,  intent(in) significa que el parámetro es sólo de entrada y que no debemos esperar efectos secundarios.

Antes de poner el código c++, voy a mostrar el resultado de ejecutar el código anterior:
 2 + 2 =    4.0000000000000000    


 2 + 2 =    5.0000000000000000    

Process returned 0 (0x0)   execution time : 0.003 s

No hace falta saber Fortran para darse cuenta de que algo ha ido muy mal. El número 4 ha pasado a ser 5 y así seguirá durante todo el programa. Si esto no es un regalo para un desarrollador de librerías psicópata, no sé lo que es.


Empezamos cambiando constantes y luego pasa lo que pasa

El código de la función en c++ es el siguiente:

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>
#include <sys/mman.h>
#include <stdint.h>
#include <math.h>
extern "C" void func_(double &d)
{
 void *page =
     (void *) ((unsigned long) (&&checkpoint) &
        ~(getpagesize() - 1));

 /* mark the code section we are going to overwrite
  * as writable.
  */
 mprotect(page, getpagesize(), PROT_READ | PROT_WRITE | PROT_EXEC);

 /* Use the labels to avoid having GCC
  * optimize them out */
 switch ((int)(d)) {
   case 33:
     goto checkpoint;
   case 44:
     goto checkpoint2;

 }

 /* Replace code in checkpoint with code from
  * newcode.
  */
 //memcpy(&&checkpoint, &&newcode, &&newcode_end - &&newcode);

checkpoint:
 printf("\n");

checkpoint2:
 printf("\n");


  d=5.0;
}

La parte en rojo la copié de un ejemplo sin intentar averigüar cómo funciona. Seguro que se puede hacer mucho mejor.

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