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.

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...

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).

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á)

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.

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.