C++ STL en paralelo

Hola,
como contábamos en una entrada anterior, últimamente hemos estado liados trabajando con archivos FASTQ de varios GBs, con decenas de millones de lecturas o reads. Para algunas de las tareas que hemos tenido que realizar, como ordenar o renombrar, nos hemos dado cuenta que lenguajes interpretados como Perl o Python son varias veces más lentos que caballos de carreras compilados como C/C++, como también comenta el autor de readfq, por ejemplo.
Por eso he estado refrescando la Standard Template Library (STL) de C++, que para alguien acostumbrado a programar en lenguajes de alto nivel es esencial, como recurso para crear estructuras de datos flexibles y eficientes en programas escritos en C/C++. Como dice Scott Meyers en su libro Effective STL, probablemente los contenedores más populares de la STL son vectores y strings. Sólo por estas dos estructuras dinámicas, a las que yo añado también las tablas asociativas (maps), vale la pena aprender a usar la STL. Sin embargo, hay mucho más que esto, y por eso escribo esta entrada, porque la STL implementada en el compilador GNU gcc (libstdc++) incluye una biblioteca de algoritmos en paralelo de propósito general, que podemos usar fácilmente en nuestros programas para sacarle el jugo a los procesadores multicore y resolver de manera más eficiente múltiples problemas. Entre otros, la versión actual de libstdc++ paraleliza los siguientes algoritmos de la STL, todos ellos útiles para tareas habituales en programación:

• std::count
• std::find
• std::search
• std::replace
• std::max_element
• std::merge
• std::min_element
• std::nth_element
• std::set_union
• std::set_intersection
• std::sort

Esquema de sort parelelo y merge posterior,
tomado de http://javaero.org/tag/parallel-merge-sort

El siguiente ejemplo de sort en paralelo, que depende de la librería OpenMP, se compila con $ g++ -O3 -fopenmp -o testP testP.cc para generar un ejecutable que empleará 4 cores para ordenar un millón de palabras de 25 caracteres de longitud:

 
#include <stdlib.h>  
#include <omp.h>  
#include <vector>  
#include <parallel/algorithm>  
using namespace std;  

// g++ -O3 -fopenmp -o testP testP.cc  
// http://gcc.gnu.org/onlinedocs/libstdc++/manual/parallel_mode.html  
#define MAXTHREADS 4;  
 
string randomStrGen(int length) 
{
   //http://stackoverflow.com/questions/2146792/how-do-you-generate-random-strings-in-c    
   static string charset = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ1234567890;:|.@";  
   string result;  
   result.resize(length);  
   for (int i = 0; i < length; i++) result[i] = charset[rand() % charset.length()];  
   return result;  
}  

int main()  
{  
   unsigned threads = MAXTHREADS;  
   omp_set_dynamic(false);  
   omp_set_num_threads(threads);  
   srand(time(NULL));  
   std::vector<string> data(1000000);  
   for(int i=0;i<data.size();i++) data[i] = randomStrGen(25);  
   __gnu_parallel::sort(data.begin(), data.end()); //std::less<string>() , std::smaller<string>() );  
   printf("%s\n%s\n%s\n%s\n",data[0].c_str(),data[1].c_str(),data[2].c_str(),data[data.size()-1].c_str());  
   return 0;  
}  

Referencias más avanzadas:
[1] http://algo2.iti.kit.edu/singler/mcstl/parallelmode_se.pdf

[2] http://ls11-www.cs.uni-dortmund.de/people/gutweng/AD08/VO11_parallel_mode_overview.pdf

Hasta luego,
Bruno

Potenciando el entorno Python I: webapps con CherryPy

Buenas!

Todos necesitamos utilizar distintos lenguajes de programación y distintas sintaxis y convenciones para muy diversas tareas. A saber: bash scripting, sed, awk, C, Perl, Python, PHP, Java, XML, HTML, JavaScript, JSON, AJAX, etcetc

Yo soy partidario de tratar de evitar esto lo máximo posible. Así que últimamente estoy tratando de potenciar mi entorno de trabajo Python y moverme a otros lenguajes sólo cuando sea realmente necesario. He hecho algunas nuevas adquisiciones a mi repertorio últimamente, aunque todavía sin mucho dominio de ellas. A pesar de esto, los resultados me han parecido muy satisfactorios. Aquí voy a hablar de un Framework para aplicaciones web en Python, CherryPy (a minimalist python web framework).

Cuales son las claves para que me haya decidido por utilizar CherryPy:
- Python no incorpora las utilidades de aplicación web que vienen incluídas en, por ejemplo, PHP. En mi caso, especialmente, control de sesión HTTP.
- Necesitaba trabajar tras un servidor Apache y tanto el módulo CGI de Python, como el mod_python de Apache son desaconsejados.
- Una aplicación CherryPy es un servidor HTTP en sí mismo, por lo que puede trabajarse muy cómodamente sobre la aplicación sin necesidad de involucrar otros sistemas en un principio.
- Permite trabajar única y exclusivamente en Python. Incluso la configuración se puede realizar con diccionarios de Python y decorators, si se prefiere esto a tener ficheros de configuración externos. La petición HTTP (la request) y sus parámetros se manejan directamente desde métodos Python.
- Es una Framework ligera, que no se inmiscuye en creación de vistas, control de usuarios, etc. Si bien permite la creación y utilización de herramientas y plugins complementarios. Siempre me han atraido las pequeñas cosas modulares, por encima de las grandes aplicaciones todo en uno. Las primeras me sugieren versatilidad y claridad; las segundas dependencia.

Quizás la parte más complicada para tirarse a trabajar con CherryPy es la configuración. Aunque más que por su complejidad, es debido a la errática documentación existente y los ligeros cambios entre versiones de la framework que para un novicio pueden ser difíciles de atribuir a una u otra versión*. En mi caso, lo mejor ha sido obviar a menudo la documentación de la página principal y trabajar directamente con el libro CherryPy Essentials y apoyándome en la lista de usuarios de CherryPy.

Sin embargo, ya digo que no es muy complicado. Una vez uno se tira a escribir sobre la aplicación y va viendo que las cosas funcionan, en seguida se pasa a trabajar con la aplicación en Python y se olvida de lo que hay detrás. Eso hasta que ejecutamos la aplicación y voilá! Ya podemos hacer peticiones desde un navegador directamente a nuestros métodos Python. Pero ¿cómo se hace una aplicación con CherryPy?

Lo primero es descargar e instalar CherryPy, para lo cual hay información detallada y fácil de llevar a cabo en la documentación oficial.
Lo segundo es crear nuestro fichero principal de Python que hará las veces de servidor. En mi caso "server.py".
A éste fichero le añadimos el import de cherrypy:

    import cherrypy

Luego hacemos el diseño de las clases y métodos que serán dianas de las distintas URLs. Por ejemplo, para:

    http://mysite/myapp/exec

con un parámetro "myparam" que viene por ejemplo de un submit en HTML o en un GET:


class NameOfMyClass(Base):
    @cherrypy.expose
    def exec(self, myparam):
        # tratar los datos de la request, por ejemplo añadirlos a sesión
        cherrypy.session['param_01'] = myparam
        # realizar operaciones
        result = mymodel.mymethod(myparam)
        # devolver una respuesta, por ejemplo en formato HTML
        return html_generator(result)

Especialmente fijarse en que el decorator (@cherrypy.expose) hace al método exec visible para el mapeador de URLs, es decir, para el cliente HTTP en definitiva.
El nombre de la clase no tiene que coincidir con el de la aplicación en éste caso, precisamente por ser la clase de la aplicación, cuyo mapeo desde URL definiremos a continuación:

    cherrypy.config.update(conf)

    cherrypy.tree.mount(NameOfMyClass(), script_name='/myapp', app_conf)

    cherrypy.engine.start()

    cherrypy.engine.block()

Faltaría añadir aquí el contenido de las variables "conf" y "app_conf". En ambos casos pueden ser rutas a ficheros de configuración o diccionarios de Python. Sobre la configuración mejor consultar aquí o el libro que ya hemos comentado. Pero un ejemplo con diccionarios podría ser:


conf = {
    'global': {
        'server.socket_host': '127.0.0.1',
        'server.socket_port': 9091,
    },
}

app_conf = {
    '/style.css': {
        'tools.staticfile.on': True,
        'tools.staticfile.filename': os.path.join(_curdir,
        '/mydeploydir/css/style.css'),
    }
}

Y con esto sería suficiente. Ejecutamos:

    python server.py

Y ya deberíamos poder hacer peticiones desde un navegador a:

    http://127.0.0.1:9091/myapp/exec?myparam=eead

Por último, tocaba desplegar la aplicación tras el servidor Apache, de forma que este siguiera sirviendo otras aplicaciones (escritas en PHP, Perl, etc.) y CherryPy sirviera mi aplicación Python. Hay muchos ejemplos sobre como desplegar CherryPy tras Apache en el libro, y también en algunas páginas en caché de Google ;) pero para un profano de Apache no quedaba claro cuál de ellas era mejor para que ambos sirvieran páginas, ni cómo modificar los ejemplos para ello. Finalmente, para el esquema que muestro a continuación, añadiendo 2 líneas a la configuración de Apache fue suficiente.

    ProxyPass /myapp http://127.0.0.1:9091/myapp/
    ProxyPassReverse /myapp http://127.0.0.1:9091/myapp/

Esto es, un proxy reverso para la aplicación creada en CherryPy con script_name en el método cherrypy.tree.mount a "/myapp".

Y esto es todo lo que quería contar sobre CherryPy y servir aplicaciones web escritas en Python. Pythonic uh? :)

oferta de empleo en Bioinformática

Puesto de Postdoc en laboratorio de Genética de Plantas del CRAG. Se trata del grupo de genómica de melón del Dr. Jordi García-Mas:

Funciones: análisis bioinformático de datos de secuenciación de RNA de melón, en el marco de un proyecto PLANT-KBBE (SAFQIM: Sugar and Fruit Quality in Melon). Además, participará en la mejora de una base de datos de genómica de melón.

·         Requisitos: PhD en Biología o similar. Gran experiencia en bioinformática y Linux.

·         Se ofrece: contrato de 6 meses, renovable por 6 meses más. 37’5h semanales.

·         Inscripción: CRAG, Ref 18/2013 (antes del 30 de junio de 2013)

FASTQ sort + parallel

Buenas,
recientemente en el laboratorio hemos estado manipulando archivos de pares de secuencias (pair-end reads) en formato FASTQ. Una de las tareas habituales ha sido limpiar los archivos de secuencias de baja calidad, ya sea recortando o eliminando directamente, para luego volver a definir parejas entre las secuencias que superaron el corte. Una estrategia posible para esta tarea es simplemente linearizar las secuencias, de manera que cada una ocupe ahora una sola línea, separando con tabuladores la cabecera, la secuencia, el separador y las calidades. Por ejemplo, la siguiente secuencia:

@SEQ_ID
GATTTGGGGTTCAAAGCAGTA...
+
!''*((((***+))%%%++)(...

quedaría así:

@SEQ_ID  GATTTGGGGTTCAAAGCAGTA...   +   !''*((((***+))%%%++)(...

Tras esta transformación ya sí es posible ordenar un archivo FASTQ con GNU sort, que viene instalado en cualquier sistema linux (y se puede instalar en Windows). GNU sort es ideal para ordenar conjuntos de datos que no caben en memoria (M), como a menudo ocurre con los archivos FASTQ, porque de manera implícita divide el problema inicial, de tamaño N, en N/M trozos que luego mezcla (merge) de manera externa.

En nuestra experiencia GNU sort es significativemente más eficiente que nuestros scripts para este tipo de problemas, puesto que ya trae de fábrica toda la lógica para partir el problema en trozos y luego mezclar las soluciones parciales. Sólo hay que tener cuidado de asignar la variable de ambiente LC_ALL, por ejemplo con:

$ export LC_ALL=POSIX

y echar a andar. Muy bien. Pero te quedas con la duda de si estás sacando el máximo partido a tu CPU multicore, podremos optimizar sort en paralelo? Y si invocamos a GNU parallel (mira el vídeo)?

Nos ponemos manos a la obra y hacemos pruebas con un archivo FASTQ linearizado real, de 576Mb:

$ ls -lh /tmp/unsortedXXpJ6CaB
-rw-------. 1 576M Jun  5 10:15 /tmp/unsortedXXpJ6CaB

Ahora lo ordenamos con GNU sort, dándole un máximo de 500Mb de área en RAM para trabajar (la M de antes):

$ time sort -k 1,1 -u -S 500M /tmp/unsortedXXpJ6CaB > /tmp/unsortedXXpJ6CaB.S
real 0m7.628s
user 0m5.143s
sys 0m2.373s

Finalmente probamos ahora con parallel, agrupando las secuencias en grupos de 100.000 elementos (ojo con esto, puedes llegar a obtener resultados parcialmente desordenados porque la segunda llamada a parallel puede recibir más argumentos de los que el shell soporta). Cambiando este valor a 10000 o a 10E6 los resultados son similares:

$ time cat /tmp/unsortedXXpJ6CaB | parallel -N 100000 --pipe --files sort -k 1,1 | 

   parallel -Xj1 sort -k 1,1 -u -m {} ';' rm {} > /tmp/unsortedXXpJ6CaB.P
real 0m15.451s
user 0m9.919s
sys 0m8.371s

Comprobamos que los resultados son idénticos:

$ diff /tmp/unsortedXXpJ6CaB.S /tmp/unsortedXXpJ6CaB.P

Conclusión de estas pruebas: no vale la pena complicarse con parallel para ordenar grandes archivos FASTQ, ya que probablemente el cuello de botella sea el merge final, y eso parece resolverlo mejor directamente GNU sort. De todos modos es posible que haya otras maneras de invocar a parallel más ventajosas
Si en vuestras pruebas obtenéis resultados distintos por favor escribid,

un saludo,

Bruno

Multicore con threads, threads::shared y Thread::Queue para controlar la RAM utilizada

Multicore con threads, threads::shared y Thread::Queue para controlar la RAM utilizada

Hace unos dias tanto Alvaro como Bruno publicaron unos blogs incluyendo scripts, tanto en Perl como en Python, en los cuales discutían sobre el uso y los problemas operativos al utilizar los modulos threads y threads::shared. Basicamente el problema principal es que, a veces es imposible utilizarlos ya que los problemas de duplicación de memoria hace que el programador no se pueda beneficiar de utilizar todos los procesadores de un ordenador, ya que se ve limitado por la cantidad de RAM que tiene disponible. Esa es una de las premisas que comentaba Bruno en su post. Yo me he encontrado multiples veces con este problema, por ejemplo procesando ficheros de genomas de Uniprot, que pueden tener ~ 30GB en sus casos mas extremos. En estos casos, al utilizar los modulos threads y threads::shared, normalmente lo que ocurre es que la memoria se va llenando mientras avanza el trabajo de las threads, y como resultado final el sistema operativo termina matando al proceso porque ocupa toda la RAM. He encontrado, mezclando código de varias personas, una forma de resolver este asunto utilizando estos modulos, y sin depender de la RAM disponible, combinando el uso de threads y threads::shared con Thread::Queue. De hecho, lo que se hace es generar una cola (en este caso una FIFO) en la cual un número predeterminado de 'worker threads' vayan sacando trabajo mientras vayan terminando. La idea es simultanear los procesos de llenado de la cola y trabajo, lo cual tambien permite controlar el tamaño de la cola, y por lo tanto la memoria utilizada. La idea sería, que esta aproximación podría combinarse con cualquier tarea que se necesite hacer, de hecho el código que compartiré a continuación, se podría adaptar fácilmente al compartido por Bruno para multi Blast.
Básicamente el programa a continuación hace lo que comentaba arriba, leer un fichero en formato Swissprot, y escanear las secuencias buscando un bias de glutamina y asparagina. El tamaño de la cola y por tanto la memoria utilizada se fija con la constante (Q_limit), en este caso definida para guardar 10000 lineas, aunque esto se puede modificar a conveniencia. Luego se generan 8 'worker threads', definidos por la constante (CORES), que van sacando trabajo mientras se llena la FIFO. Luego las 'worker threads' devuelven predicciones de secuencias con este tipo de bias y el programa principal las imprime en un fichero texto. Basado en este ejemplo, es muy facil adaptar esto a cualquier tarea que implique escanear secuencias en la búsqueda de cualquier propiedad, en ficheros enormes. Por ejemplo, un programa similar a este procesa en aproximadamente 2 horas un fichero de 35 GB con 64 CPUs y una utilización máxima de RAM de 5 GB.
Bueno aqui va el código, y espero sea útil la idea:

#!/usr/bin/perl -w

#Loading all the libraries to run in a multithreading environment
use warnings;
use threads;
use threads::shared;
use Thread::Queue;

my %aa_scores = (
   'A' => -0.149,
   'C' => -1.000,
   'D' => -0.395,
   'E' => -0.726,
   'F' => -0.125,
   'G' => 0.010,
   'H' => -0.034,
   'I' => -0.397,
   'K' => -0.496,
   'L' => -0.408,
   'M' => 0.044,
   'N' => 0.659,
   'P' => 0.059,
   'Q' => 0.536,
   'R' => -0.314,
   'S' => 0.192,
   'T' => -0.070,
   'V' => -0.450,
   'W' => -0.908,
   'Y' => 0.206,
   'O' => 0,
   'U' => 0,
   'B' => 0,
   'Z' => 0,
   'X' => 0,
);

use constant WINDOW => 60;
use constant CUTOFF => 12;
use constant CORES => 8; #Set the number of cores to use
use constant Q_limit => 10000; #Set the number of elements in the Queue

my $Q = new Thread::Queue (); #Initializing the Queue obect
my @threads;

#Loading the Swissprot file from command line argument
my $dat_file = shift;
$/ = "//\n"; #Reading the multiple lines of a Swissprot entry in a single line. The (//) line corresponds to the entry end makr 
#As this huge files are ussually distributed compressed here I pipe with zcat to avoid decompresing prior to processing
open (DAT_FILE, "zcat $dat_file |") or die "Sequence file couldn't be opened for reading\n";

#First generating the Queue
my $builder = threads->create (\&Queue_generator);

#And then setting up the worker threads
push @threads, threads->create (\&parse_swissprot) for (1 .. CORES);

#Wait until all the threads finish their work
$builder->join;

#Now that the work is done, collect the returned info extracted by the worker threads
#It is neccesary to fill a hash with all the partial results to be able to print the
#results with the scanning of the complete file in an organized fashion
my %domains;
foreach my $thread (@threads)
{
   my %partial_results = %{$thread->join};
 
   while (my ($organism, $predictions_ref) = each (%partial_results))
   {
      my @predictions;
      while (my ($seq_id, $info_ref) = each (%{$predictions_ref}))
      {
         my $score = $info_ref->{'Score'};
         my $domain = $info_ref->{'Seq'};
         my $window_pos = $info_ref->{'Window'};
   
         my $protein_id = $1 if ($seq_id =~ m/^>(\w+);/);
         my $protein_ac = $1 if ($seq_id =~ m/^>$protein_id;\s(\w+);\s/);
   
         push @predictions, "\t$protein_id|$protein_ac\tWindow Position=$window_pos; Score=$score | Predicted Domain: $domain\n";
      }
      push @{$domains{$organism}}, @predictions;#Loading all the predictions returned by each thread for a given organism
   }
}

  
#And finally printing the predictions of domains with a specific amino acid bias in the complete file
#organized by organisms
open (OUTFILE, ">predictions.txt");
while (my ($organism, $predictions_ref) = each (%domains))
{
   my @predictions = @{$predictions_ref};
 
   my $total_predictions = scalar (@predictions);
   print OUTFILE ">$organism: Total=$total_predictions\n@predictions\n";
}
close (OUTFILE);
close (DAT_FILE);

#Here I construct the Queue. The idea is that while the generator is feeding the queue other working threads
#are retrieving (emptying) from the queue and the memory used can be controlled during execution
sub Queue_generator
{
   while (my $line = <DAT_FILE>)
   {
      chomp ($line);
      while (1)
      {
         #The line is only send to processing if the number of items in the queue is less than a fixed size defined by Q_limit
         #thus a line is trapped in this infinite loop until the queue has space to receive it
         if ($Q->pending () < Q_limit) 
         {
            $Q->enqueue ($line);
            last; #Once the line has been loaded in the queue we need to break out of the infinite loop to read the next line from the file
         }
      }
   }
 
   #Once the file is exhausted it is needed to alert the worker threads that the work is over
   #When a thread catch and UNDEF variable it stops working
   $Q->enqueue (undef) for (1 .. CORES);
}

#This is the code of the worker threads
sub parse_swissprot
{
   my %predictions; #The variable to return to the main program
 
   #Here the thread DEQUEUE one line from the queue stack
   while (my $line = $Q->dequeue)
   {
      my @patterns = ();
      #Getting the Entry Name
      my $seq_id = $1 if ($line =~ m/^ID\s+(\w+)\s+/);
      
      # And now getting the Accession Number lines
      @patterns = $line =~ m/^AC\s+(.*)/mg;
      # and finally just getting the primary AC -i.e. the first of all the possible multiple ACs-
      my $seq_ac = $1 if ($patterns[0] =~ m/^(\w+);/);
      
      #And now getting the rest of the Entry info
      @patterns = $line =~ m/^DE\s+(.*)/mg;
      my $description = join ' ', @patterns;
      @patterns = $line =~ m/^SQ\s+(.*)/mg;
      my $seq_info = join ' ', @patterns;
      @patterns = $line =~ m/^OS\s+(.*)/mg;
      my $organism = join ' ', @patterns;
      @patterns = $line =~ m/^\s+(.*)/mg;
      my $seq = join '', @patterns;
      $seq =~ s/\s+//g;
      
      my $id = '>' . $seq_id . '; ' . $seq_ac . '; ' . $description . ' ' . $seq_info . "  [$organism]";
      
      next unless (length ($seq) >= WINDOW);
      my %results;
      
      #Loop to scan the complete sequence of the corresponding entry
      #After scanning the sequence only the highest scoring stretch is returned
      for my $i (0 .. (length ($seq) - WINDOW))
      {
         my $domain = substr ($seq, $i, WINDOW);
         my @domain = split (//, $domain);
   
         my $domain_score = 0;
         foreach my $aa (@domain)
         {
            $domain_score += $aa_scores{$aa};
         }
         $domain_score = sprintf ("%.3f", $domain_score);
   
         #Updating the $results variable to store only the highest scoring stretch
         if (exists ($results{$id}))
         {
            $results{$id} = {
               'Score' => $domain_score,
               'Seq' => $domain,
               'Window' => $i,
               'CompleteSeq' => $seq,
            } if ($domain_score > $results{$id}{'Score'});
         }
         else
         {
            $results{$id} = {
               'Score' => $domain_score,
               'Seq' => $domain,
               'Window' => $i,
               'CompleteSeq' => $seq,
            };
         }
      }
  
      #Loading the entry in the variable returned to the main program only if the highest scoring stretch passes the cutoff
      if ($results{$id}{'Score'} >= CUTOFF)
      {
         $predictions{$organism}{$id} = $results{$id};
      }
      $seq_id = $seq_ac = $description = $seq_info = $organism = $seq = undef;
   }
   return (\%predictions);
}