Introducción

El PDL (Perl Data Language) es un conjunto de módulos que le dan a Perl la facultad de realizar operaciones matemáticas en conjuntos matriciales de datos de cualquier dimensión. Y, según los recursos disponibles, de cualquier tamaño. Se podría decir que es una librería matemática para el lenguaje Perl.

Por qué usar PDL

Ya sabemos que nuestro querido Perl es un lenguaje interpretado y, por tanto, tiene esa mala fama de ser lento. Me refiero al momento de la ejecución del programa. Sabemos que si tenemos que procesar miles, millones o incluso miles de millones de datos, el hacerlo con un programa interpretado será desperdiciar recursos, frente a una solución compilada, como puede ser C o C++.

Pero Perl tiene algo a lo que no podemos resistirnos: su rapidez en la producción de una solución. Perl nos permite tener un programa listo para correr en muy poco tiempo. No hay que compilar nada y su lenguaje de muy alto nivel nos permite dar una solución en muy pocas líneas de código, según nuestro nivel de experiencia con él. Así que tendemos a resolver todos los problemas con Perl. Y con los ordenadores modernos, una solución interpretada es válida en muchísimos casos.

Entonces, si Perl es interpretado, no puede ayudarnos a resolver un problema que requiere una solución compilada para la máquina en la que se va a correr. Como mucho, Perl podría llamar a un programa de procesado, un programa externo a nuestro script. Pero eso se traduce también en una solución lenta, ya que se pierde un poco de tiempo, en el sistema, para llamar a un nuevo programa, pasarle los argumentos y esperar su resultado. Lo ideal es que Perl pueda tener su propio librería compilada de procesado de datos y ser él mismo el que haga todos los procesos.

PDL viene a llenar este hueco. Nos permite manejar millones y millones de datos con la sintaxis de Perl.

Ejemplo: Nubosidad

Se describe a continuación un ejemplo de aplicación de Perl con la librería PDL. El ejemplo es de grado medio, medio-alto. Aunque el lenguaje Perl utilizado es sencillo, la parte de PDL es muy especial, pero, con sólo leerlo, puede dar una idea de lo que se puede hacer con PDL. No se trata de un ejemplo para empezar a trabajar con PDL.

El ejemplo propuesto trata de resolver el siguiente problema:
Trabajamos en un departamento que se dedica a la investigación del clima y la vegetación.

Una de las rutinas normales del trabajo es la del procesado de imágenes de satélites de baja órbita. Estos satélites se caracterizan por una alta resolución (250m/pixel). Uno de las productos que se pueden sacar de estas imágenes es el nivel de nubosidad que existía en el momento en que el satélite hizo la pasada por encima de una determinada zona a analizar.

Definimos nubosidad como el porcentaje de píxeles de una comarca que están despejados de nubes con respecto al total de píxeles que cubre esa comarca. Saber la nubosidad es importante porque luego podremos deducir cálculos de índice de vegetación y riesgo de incendios.

La misión es: tenemos una serie de imágenes, todas ellas de la misma zona del planeta, de un tamaño de 6144 por 4096 píxeles, y de dos bytes por pixel cada uno, donde esos bytes pueden ser 0 si en ese pixel había una nube y un valor mayor que 0 si estaba despejado. Hay que calcular el índice de nubosidad para todas las comarcas de esa zona, por cada una de las imágenes.

Naturalmente, para saber dónde están las comarcas necesitamos un mapa. El mapa es otra imagen del mismo tamaño que las anteriores, donde cada pixel es de un byte, con un valor que representa el número de comarca. Si hubiese un 0, es que ese pixel pertenece a otro país o al mar.

Problemas...

Bueno, como hay que procesar una serie de imágenes y sacar un informe para cada una de ellas, enseguida pensamos en Perl. Pero... ¡cada imagen tiene más de 24 millones de píxeles!. Aunque no hay que hacer ninguna operación matemática, el simple hecho de recorrer 24 millones de puntos en un bucle en Perl... puede hacerse eterno (desde el punto de vista del ordenador). Y el enunciado del problema no comenta el número de imágenes a procesar... pueden ser dos o doscientas... Y si hay otra persona o programa esperando los resultados... mal vamos. Hay que buscar una solución más rápida.

PDL al rescate

La operación a realizar en todos los píxeles es sencilla: saber cuántos píxeles de una zona son mayores que 0. Y luego dividirlo por el número de píxeles del total de la zona.

Lo vamos a expresar bien: queremos saber lo despejada que es una comarca, es decir, cuántos píxeles de la imagen que pertenecen a esa comarca son mayores que 0. Y esto, repetirlo para todas las comarcas y para todas las imágenes.

Empezamos

El código que resuelve el problema está al final del artículo. Lo vamos a comentar línea a línea.

Al principio, en la línea 1, activamos el sistema de avisos (warnings) con la opción -w.
Se podría haber colocado un use warnings; en el código, pero poniéndolo arriba del todo, nos avisará de cualquier posible fallo, incluído en todos los módulos que usemos.

En las líneas 8 y 9 incorporamos a nuestro programa toda la sabiduría de los módulos PDL y PDL::IO::FastRaw. El primero es obligatorio si vamos a usar alguna de las funciones de la librería PDL, y el segundo lo necesitamos para poder leer ficheros de forma muy rápida.

La línea 10 nos obliga a programar de forma cuidadosa...

De la línea 12 a la 15 inicializamos algunas constantes. Se podrían haber definido con use constant ... pero el programa es tan pequeño que no merecía la pena tener tanto cuidado. No vamos a cambiarlas a lo largo del programa.
Aparte del $ancho y $alto de las imágenes y del $directorio donde están almacenadas, tenemos puesto el nombre del fichero que guarda la posición de cada una de las comarcas.

En la línea 17 nos colocamos en nuestro directorio de trabajo... o morimos en el intento.
En la 19, quitamos el buffer de salida. Todos los mensajes de los print saldrán inmediatamente a pantalla (y fichero).

De la 20 a la 22 tenemos las variables que vamos a usar:

@comarcas es un array que guarda el código de cada comarca. Por ejemplo, $comarca[32] puede tener el código 452, indicando que ese es el número de código para la comarca 31ª. @indices_comarcas es un array que guarda apuntadores a cada uno de los píxeles que forman las comarcas. Esto hay que explicarlo un poco: necesitamos saber dónde está (qué píxeles forman) la comarca 31ª. Pues en $indices_comarcas[32] guardaremos un objeto que tendrá todos esos píxeles. Pueden ser un píxel o un millón. Nos da igual. Realmente, lo que estamos guardando es la posición (offset) desde el principio de la imagen de todos esos píxeles. El primer pixel de la imagen será el 0, luego el 1, etc. etc. hasta el 24 millones. Nosotros guardamos esas posiciones. Una comarca que tuviera sólo 3 píxeles, podrían ser por ejemplo, el 653711, el 653712 y el 659855 (por estar en la línea de abajo de los otros dos). @pixeles_comarcas es un array para ahorrar un poquito de tiempo... va a guardar cuántos píxeles hay por cada comarca. Podríamos habernos ahorrado esta variable, pero ya digo que es por ahorrar un poco (realmente muy poco) tiempo de procesado.

Comenzamos con los preparativos

En la línea 25 nos leemos el fichero del mapa de comarcas.
Utilizamos la función mapfraw, sacada del módulo PDL::IO::FastRaw. Esta función NO lee un fichero. Lo que hace realmente es mapear el fichero a la memoria, es decir, accederemos a la memoria reservada como si estuviera ocupada por el contenido del fichero, pero en realidad, lo que estará haciendo el sistema operativo es acceder al disco duro. Esto nos permite acceder a un fichero muy grande (incluso a veces mucho más grande que la propia memoria del ordenador) sin tener que leerlo enteramente a memoria (porque a veces es imposible, claro).
Para 'leer' el fichero, indicamos las dimensiones de la imagen ($ancho*$alto), el formato de cada pixel (ushort quiere decir un entero de dos bytes, sin signo) y, además, indicamos que no vamos a modificarlo (ReadOnly=>1).
Notar algo especial... hasta ahora hemos hablado de imágenes, zonas, comarcas... por lo que estamos pensando en 2 dimensiones (estamos trabajando en superficies)... pero... en nuestro problema, las operaciones que vamos a realizar son a nivel de píxel: lo que ocurre en un pixel no influye en el de al lado... ergo no es necesario tratar los datos como si estuvieran en dos dimensiones. En vez de leer una imagen de dimensiones $ancho y $alto, tenemos en su lugar una matriz de una dimensión, de tamaño $ancho*$alto. Y así trabajaremos con todas las imágenes. Esto acelera mucho más los cálculos (aunque con PDL, se nota poco la diferencia).
Ya se acabó la línea 25... una sola línea que da mucho juego...

En la línea 28 vamos a averiguar todos los códigos de las comarcas. Vamos a usar dos funciones de PDL: uniq y list. uniq nos va a devolver un piddle (ahora mismo lo explico) que contiene todos los valores distintos (únicos) del piddle $imagen_comarca. list transforma ese piddle en una lista normal de Perl.

Ahora es cuando nos encontramos con un nuevo objeto llamado piddle... ¿Qué es un piddle? Pues un piddle es un objeto PDL. Es una 'entidad matemática' que contiene datos numéricos, en una serie de dimensiones. Podría ser un simple valor escalar (un dato en una dimensión), o un hipercubo de un millón de datos, en 6 dimensiones.

Al final, como resultado, tenemos en @comarcas todos los números distintos que hemos encontrado en el mapa de comarcas. Y además, ordenados, de menos a mayor.

En la línea 29 quitamos el primer valor del array @comarcas, pues siempre vale 0. Resulta que nos hemos enterado que en el mapa de comarcas sí aparecen valores 0, representando píxeles que son el mar. Y el mar no es ninguna de nuestras comarcas, así que lo quitamos.

De la línea 31 a la 40 vamos a sacar una serie de estadísticas, que nos ayudarán más tarde, comarca por comarca.
En la línea 34 averiguaremos dónde está esa comarca. Usaremos la función which, del PDL. Lo que hace esa función es buscar por toda la $imagen_comarca aquellos píxeles que tienen el código de la comarca que estamos analizando. El resultado es un piddle que guarda las posiciones a todos esos píxeles. Ese piddle lo guardamos dentro del array @indices_comarcas.

Después de saber dónde está la comarca, en la línea 37 sabremos cuántos píxeles la forman. Llamamos de nuevo a una función PDL, nelem, que nos devuelve esa cantidad, y la guardamos en el array @pixeles_comarcas.

La línea 37 sólo es informativa, por si nos interesa ver qué es lo que está haciendo nuestro programa.

Empieza la magia

El bucle que va desde la línea 42 hasta el final recorre todas las imágenes a procesar. Cada $fichero sale del operador diamante, que se utiliza en este caso como un patrón de nombres de ficheros con un carácter comodín. Este patrón se pasa al sistema operativo y este nos devolverá todos los nombres de ficheros que coincidan con ese patrón. Repito: el comodín (*) que usamos ahí no tiene significado en Perl, sino para el sistema operativo.

Después de informar al usuario de la imagen que estamos procesando en ese momento, en la línea 47 la leemos, de la misma forma que hemos hecho antes con la imágen de comarcas. El resultado se guarda en el piddle $imagen_nubosa. Este piddle será una matriz de una dimensión cuyos valores serán 0 para el caso de que el pixel tenga una nube encima o un valor superior si está despejado.

En las líneas 51 y 52 vamos a abrir el fichero de texto donde guardaremos el resultado del procesado del $fichero. El nombre del fichero de texto lo derivamos del propio nombre del $fichero original de la imagen, cambiando la extensión 'img' por 'txt'.

Bueno, ya tenemos cargada una imagen. Hay que procesarla. El proceso significa que recorremos todas las comarcas y calculamos su índice de nubosidad.

Pues eso hacemos... el bucle de las comarcas empieza en la línea 54.

Y aquí llegamos a la magia...

De la línea 56 a la 62 -que en realidad, sólo es UNA- averiguamos cúantos píxeles de esa comarca están despejados.

La información la tenemos, pero está repartida entre varios piddles y arrays. Sólo hay que acceder a ellos de forma correcta.

Primero leamos los comentarios que hay puestos a la derecha del código, pues todos juntos forman la frase que resuelve el problema, y cada renglón corresponde con el código Perl que tienen a la izquierda. Es un conjunto de llamadas a funciones PDL. Vamos a describir el funcionamiento en el mismo sentido en que Perl las ejecuta.

En la línea 60, $indices_comarcas[$num_comarca] es un piddle que nos devuelve las posiciones de los píxeles de la comarca que estamos analizando.

En la línea 59, con la función index, 'sacamos' los píxeles de la $imagen_nubosa que acabamos de leer que corresponden a esas posiciones. Por si alguien no se ha dado cuenta: las posiciones de los píxeles las hemos sacado del fichero de comarcas. Ahora las usamos en OTRA imágen, la nubosa que hemos leído. Si suponemos que los dos ficheros tienen las mismas dimensiones, lo que estamos haciendo es mirar en la imágen nubosa los píxeles que coinciden en posición con los mismos píxeles de la comarca. El resultado es otro piddle, que en este caso contiene esos píxeles (puede parecer un poco vago decir esto, pero en realidad, eso es lo que ocurre: 'tenemos' esos píxeles).

Ahora, en la línea 58, preguntamos cuáles de esos píxeles son mayores que 0 (que está en la línea 62, sólo hay que seguir la lista de paréntesis). El resultado es un piddle que contiene los píxeles que están despejados (son mayores de 0, en nuestro problema). Sólo falta que en la línea 57, nelem nos devuelva el número de elementos de ese piddle.

Ya está... de los 24 millones de puntos, sólo miramos aquellos que corresponden a la comarca. Y de esos, cuáles y cuántos son los que son mayores de 0.

Naturalmente hay otras formas de hacerlo. En C o C++ lo que haríamos sería recorrer toda la imagen e ir mirando pixel por pixel a qué comarca pertenece (del fichero de mapas) y si está despejado (del fichero de imagen nubosa). Se puede hacer todo a la vez y al final, sacar los resultados. Incluso hasta podría ser más rápido... pero no tan bonito como aquí.

Probé a realizarlo en C y la mejora de velocidad no es significativa. Lo que sí fué notable fue el tiempo que tardé en hacerlo en C.

En la línea 64 calculamos el índice de nubosidad. Una simple división entre los píxeles que están despejados con respecto al total de píxeles que forman la comarca.

Finálmente, sacamos los valores calculados a pantalla (línea 67) y al fichero de texto (línea 71). El resultado que nos piden es que a de ser un fichero de texto con dos columnas. En el primero, el código de comarca y en la segunda, el porcentaje de nubosidad, multiplicado por 100 y con un decimal de resolución.

Espero que con este humilde programa os haya picado la curiosidad de esta fantástica librería.

Ahora comentaros los aspectos negativos...

La curva de aprendizaje de PDL es muy fuerte. Pero mucho. Muchísimo más que el propio Perl.

Me he encontrado muchas veces completamente perdido. La documentación es escasa (aparte del manual y un libro), a veces poco actualizada. La mejor manera de solucionar un problema es participar en la lista de correo de PDL. Hay gente que, muy amablemente, te echarán una mano y te asombrarás de las soluciones que te ofrecen.

Realmente, lo más difícil de esta librería es saber qué función debes usar en cada momento. Con el tiempo aprender a usar unas cuantas de ellas, pero las distintas formas de combinarse entre ellas y otro aspecto importante, la forma de acceder a los datos, es lo más complicado. En cierto sentido, si usas piddles en 3 dimensiones, tiendes a pensar en 3 dimensiones para resolver el problema.

Y no digamos si necesitas más dimensiones. Me encontré con un problema que requería 8 dimensiones. Con la calculadora en la mano, me salían tiempos de procesado de más de una semana, con un programa en C. Al final de muchos días de pruebas y programación, con Perl y PDL, el programa tarda 8 horas en procesar 521.838.526.464.000 datos. Asombroso.

Nuestro trabajo figura como caso de éxito en la página de PDL, entre otros maravillosos ejemplos.

Desde hace unos meses que ya no sigo allí y que no necesito usar PDL, pero algunas cosas no se olvidan fácilmente. PDL se olvida muy rápido -por lo complicado que es-, pero nó sus efectos.

Cortesía de Joaquín Ferrero

Joaquín Ferrero es programador y moderador del foro de Perl en Español, y gentilmente ha ofrecido algunos de sus tutoriales para su publicación en este sitio. Puedes visitar la página personal de Joaquín Ferrero con el siguiente link: http://www.joaquinferrero.com/ O puedes ver su perfil de "Perl en Español".