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 a la velocidad de un programa compilado.
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 de una comarca 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 dos bytes, 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 y por lo tanto, no nos interesa.
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 complicada, 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 está 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 de la imagen y para todas las imágenes restantes.
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. Lo quitaremos cuando estemos seguros de que no hay fallos en nuestro programa.
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... lo mismo de antes, lo quitaremos en cuanto estémos seguros de que hemos terminado un programa correcto.
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[30] 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[30] guardaremos un objeto que tendrá todos esos píxeles. Pueden ser uno o un millón de píxeles. Nos da igual. Realmente, lo que estamos guardando es la posición (offset) desde el principio de la imagen a 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.
* @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, nosotros accederemos al mapa de comarcas como si lo hubiéramos leído a memoria, pero en realidad, lo que estará haciendo el sistema operativo es acceder al disco duro, al fichero del mapa de comarcas, de forma directa. 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 (porque a veces es imposible por el tamaño del propio fichero, 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).
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. Da igual uno u otro. Nosotros lo manejaremos en nuestro programa 'casi' de la misma manera que trabajamos con cualquier otra variable escalar de Perl.
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 por cada una de ellas.
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. Sí... lo hacemos TODO en una sola línea.
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 dicha unos párrafos más arriba, y que resuelven el problema, y cada renglón corresponde con el código Perl que tiene 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 y que corresponden a la misma zona del planeta, 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 otro piddle que contiene los píxeles que están despejados. Sólo falta que en la línea 57, nelem nos devuelva el número de elementos de ese piddle.
Ya está... el proceso se resume en que, de los 24 millones de puntos de la imagen, sólo miramos aquellos que corresponden a la comarca que estamos analizando. Y de esos pocos píxeles, 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 La enorme ventaja es que miramos los píxeles de las comarcas diréctamente, pasando por alto todos los píxeles 0 de la imagen.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 ha 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.
Using perl Syntax Highlighting
- #!/usr/bin/perl
- #
- # Cálculo de la nubosidad.
- # Programa para el cálculo de los píxeles nubosos de las comarcas.
- # Joaquín Ferrero. 2004.
- #
- use PDL;
- use PDL::IO::FastRaw;
- use strict;
- my $ancho = 6144; # Constantes
- my $alto = 4096;
- my $dir = '/almacen1/nubosidad';
- my $fichero_comarcas = 'COMARCAS.img';
- chdir $dir or die "No puedo entrar en el directorio de trabajo: $!";
- $|++; # Salida a pantalla sin buffer
- my @comarcas; # Números de las comarcas
- my @indices_comarcas; # Índices de píxeles a cada comarca
- my @pixeles_comarcas; # Número de píxeles por cada comarca
- print "Leyendo imagen de comarcas\n";
- my $imagen_comarca = mapfraw( $fichero_comarcas, { Dims => [ $ancho * $alto ], Datatype => ushort, ReadOnly => 1 } );
- @comarcas = list $imagen_comarca->uniq; # Valores únicos de las comarcas
- shift @comarcas; # Quitamos el '0', que siempre está
- foreach ( 0 .. $#comarcas ) {
- # Dónde está cada comarca
- $indices_comarcas[$_] = which( $imagen_comarca == $comarcas[$_] );
- # Cuantos píxeles hay por comarca
- $pixeles_comarcas[$_] = $indices_comarcas[$_]->nelem;
- printf "Comarca %-3d : %-3d : %d\n", $_, $comarcas[$_], $pixeles_comarcas[$_];
- }
- foreach my $fichero (<ng_*img>) # Para todos los ficheros
- {
- print "Procesando fichero $fichero...\n";
- # Leemos la imagen filtrada por nubes
- my $imagen_nubosa = mapfraw( $fichero, { Dims => [ $ancho * $alto ], Datatype => ushort, ReadOnly => 1 } );
- # Abrimos fichero de texto destino
- ( my $texto = $fichero ) =~ s/img$/txt/;
- print "Salida en fichero $texto\n";
- open SALIDA, ">$texto" or die "No puedo escribir a $texto: $!";
- foreach my $num_comarca ( 0 .. $#comarcas ) # Para cada comarca
- {
- my $despejado = # Cuánto está despejado
- nelem( # Cuántos píxeles,
- which( # de
- $imagen_nubosa->index( # la imagen nubosa
- $indices_comarcas[$num_comarca] # que pertenecen a la comarca,
- ) > 0 # son mayores que 0 (están despejados)
- )
- );
- my $ratio = $despejado / $pixeles_comarcas[$num_comarca];
- printf
- "%s comarca %3d: %5d/%5d = %05.1f%%\n",
- $fichero, $comarcas[$num_comarca],
- $despejado, $pixeles_comarcas[$num_comarca], 100 * $ratio;
- printf SALIDA "%03d %05.1f\n", $comarcas[$num_comarca], 100 * $ratio;
- }
- close SALIDA;
- }
- # Fin
Coloreado en 0.006 segundos, usando GeSHi 1.0.8.4
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 aprendes 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.
Quizás lo mejor sea que, al cabo de un tiempo de pensar en objetos en 3 dimensiones, tiendes a resolver el problema pensando de la misma forma. Mientras que en otros lenguajes como C tendemos a simplificar el problema hasta el mínimo detalle para poder resolverlo, con PDL la mejor solución es pensar que en una variable escalar está todo nuestro objeto n-dimensional.
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 y técnicas de caching, 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 no necesito usar PDL, pero algunas cosas no se olvidan fácilmente. PDL quizás se pueda olvidar rápidamente -por lo complicado que es-, pero nó sus efectos, que te dejan unos buenos recuerdos.
con las cifras que expones en tu articulo, no tengo un manejo fluido de módulos en perl y por el momento no creo poder darle una utilidad al PDL en la cual se destaque su desempeño, pero se nota que vale mucho la pena probarlo he intentar implementarlo, muchas gracias por la explicación tan sencilla a un problema tan complejo... solo me gustaria hacerte una pregunta personal, si no es molestia, ese es trabajo de algun grupo de investigación de la universidad de valladolid? 
