19 de abril de 2012

perl + GNU Scientific Library

Hola,
leyendo el Linux Journal me he encontrado con este artículo que explica como usar la GNU Scientific Library en C, y me he acordado de que tenemos por el laboratorio, sin estrenar, la versión 1.12 del manual. Para qué vale esta librería? La respuesta corta es: para hacer cálculos de todo tipo con algoritmos eficientes que minimizan los errores de redondeo de los procesadores digitales.
A veces una librería es la mejor forma de emplear algoritmos extensivamente optimizados y que tal vez no sepamos cómo programar. A día de hoy la SGL incluye funciones en estas áreas:
Complex Numbers   Roots of Polynomials
Special Functions   Vectors and Matrices
Permutations   Sorting
BLAS Support   Linear Algebra
Eigensystems   Fast Fourier Transforms
Quadrature   Random Numbers
Quasi-Random Sequences   Random Distributions
Statistics   Histograms
N-Tuples   Monte Carlo Integration
Simulated Annealing   Differential Equations
Interpolation   Numerical Differentiation
Chebyshev Approximation   Series Acceleration
Discrete Hankel Transforms   Root-Finding
Minimization   Least-Squares Fitting
Physical Constants   IEEE Floating-Point
Discrete Wavelet Transforms   Basis splines
En algunos aspectos puede ser más limitada que las Numerical Recipes, pero al liberarse bajo GPL tenemos mayor libertad de uso, y sin pagar un duro. Además, gracias al trabajo de Jonathan Leto, disponemos en CPAN de una interfaz para Perl, llamada Math::GSL. Algunas de las aplicaciones de la librería se pueden conseguir con otros módulos de CPAN, pero a priori la GSL tiene a su favor: i) menor tiempo de ejecución y ii) mayor precisión numérica, una limitación de Perl ya discutida.

Pasemos a la práctica: cómo se instala esto? En mi Ubuntu 10.4 fue tan sencillo como decir $ sudo cpan -i Math::SGL y luego ir aceptando la instalación de sus dependencias.

Podemos ahora probar un ejemplo donde calcularemos una derivada numérica, mostrando algunas de las capacidades de la librería:


 #!/usr/bin/perl -w  
 # Ejemplo de derivada de ln(x) con GSL  
 # Adaptado de ejemplos de Jonathan Leto en:  
 # https://github.com/leto/math--gsl/blob/master/examples/deriv/basic  
   
 use strict;  
 use Math::GSL::Deriv qw/:all/;  
 use Math::GSL::Errno qw/:all/;  
   
 # incremento de derivada  
 my $h = 0.01;   
   
 # queremos X = pi, ejemplo constante matematica en GSL   
 my $x = $Math::GSL::Deriv::M_PI;   
    
 # derivada numerica, recuerda en perl log = ln  
 my ($status,$val,$err) = gsl_deriv_central ( sub { log($_[0]) }, $x, $h);  
   
 # derivada analitica, demostracion en  
 # http://www.math.com/tables/derivatives/more/es-ln.htm  
 sub dlndx{ return 1/$_[0] }  
   
 if($status == $GSL_SUCCESS)   
 {  
   printf("deriv(ln((%g)) = %.18g, error maximo esperado = %.18g\n", $x, $val, $err);  
   printf("dlndx(%g)   = %.18g\n" , $x, dlndx($x));  
   printf("error observado   = %.18g\n",abs($val-dlndx($x)));  
 }   
 else   
 {  
   my $gsl_error = gsl_strerror($status);  
   print "ERROR: $gsl_error (es derivable log(x) en ese punto?)\n";  
 }  

Se obtiene algo como:

deriv(ln((3.14159)) = 0.31830988618650752, error maximo esperado = 5.7425388097458009e-11
dlndx(3.14159)      = 0.31830988618379069
error observado     = 2.7168267635602206e-12

El autor del módulo tiene en la web una colección de ejemplos que sin duda amplian la descripción de las capacidad de GSL. Por ejemplo, en este ejemplo se explica como calcular productos escalares de dos maneras, encontrando que la implementación del módulo Math::GSL::BLAS es al menos el doble de rápida.

Ya que este ejemplo no representa las tareas típicas de la programación en Bioinformática os muestro un segundo mucho más terrenal, donde ordenamos vectores de números reales, comparando la implementación de mergesort de Perl con la heapsort inestable de GGL. Hay un capítulo del manual dedicado sólo a esto, pero probablemente la documentación del módulo Math::GSL::Sort sea suficiente:

 use Math::GSL::Sort qw/:all/;  
   
 my $numbers = [ map { rand(100) } (1..100000) ];  
 my $sorted = gsl_sort( $numbers,1,scalar(@$numbers) );     
 my @sortedp = sort {$a<=>$b}( @$numbers );  

Si hacéis pruebas veréis que la implementación GSL resulta ser más rápida, en mi caso aproximadamente un 20%.

Otras aplicaciones posibles incluyen el cálculo de permutaciones y combinaciones, algo ya discutido en este blog,
hasta otra,
Bruno

PD Si alguien ha probado a instalar este módulo en sistemas Windows le agradeceremos sus comentarios

1 comentario:

  1. Viendo la matriz de pruebas de instalación de Math::GSL ( http://matrix.cpantesters.org/?dist=Math-GSL+0.26 ) vemos que muy pocos se han atrevido a instalar este módulo en Windows.

    Algo que no sorprende, desde luego ;)

    En cuanto a Math::GSL, decir que también se puede usar desde la biblioteca matemática PDL: http://cpan.uwinnipeg.ca/htdocs/PDL/Modules.html#GNU_SCIENTIFIC_LIBRARY
    con la consiguiente ventaja en la velocidad de tratamiento de objetos matemáticos, de forma directa, y con conversión a objetos Perl solo al final del proceso.

    ResponderEliminar