may 05 2008

Algunos paquetes científicos de Python

Tag: Herramientas y Programaciónadmin @ 6:09

Aquí va una lista de algunos paquetes que poco a poco van haciendo de Python un lenguaje más apto para uso en un contexto astrofísico, o científico en general. Más adelante añadiremos alguna entrada más detallada comentando algunos de estos paquetes.

  • Matplotlib: Potente librería gráfica para la realización de gráficos en 2-D de gran calidad, que utiliza una sintaxis similar a la de MATLAB. Permite hacer histogramas, gráficos de puntos, barras, curvas de nivel, degradados, etc. (ver aquí algunos ejemplos). Maneja diferentes formatos gráficos, por ejemplo png y eps.
  • ParselTongue: Se trata de una interfaz mediante el uso de Python para el software AIPS, que sirve para la reducción y análisis de datos en radioastronomía. Añade la funcionalidad y comodidad de los scripts en Python a la hora de reducir datos con AIPS, porque aunque éste sea un programa potente y ampliamente difundido ente la comunidad radioastronómica, posee una interfaz un tanto rígida.
  • PyEphem: Módulo para realizar cálculos astronómicos tales como posiciones de astros, planetas, horas de orto y ocaso de los astros, etc. También se puede calcular la posición de cometas, asteroides y satélites una vez proporcionados los elementos orbitales de los mismos.
  • PyFITS: Pequeño pero útil módulo que permite usar desde Python los ficheros con formato FITS, tanto tablas como imágenes. Se puede acceder a la información de las cabeceras y manipularla, trabajar con las imagenes convirtiendolas a una matriz para operar con ellas, etc.
  • PyGTK + Glade: Una combinación de herramientas, no necesariamente para uso científico, que permite diseñar interfaces gráficas para aplicaciones. PyGTK utiliza las librerías GTK+ desde python, y el programa Glade crea el esqueleto de la aplicación, con ventanas, botones, menú, etc.
  • PyRAF: Interfaz desde Python del programa IRAF que sirve para análisis y reducción de datos astronómicos.
  • Scipy: Paquete científico de uso general que incluye rutinas para realizar cálculo numérico, integración, optimización, estadística y un largo etcétera. Está basado en el paquete Numpy que es el sucesor de Numeric y numarray como paquete encargado del uso eficiente de matrices en Python.

dic 11 2007

Espectropolarimetría con Python

Tag: Herramientas y Programaciónfranwerst @ 7:21

Aparece hoy en arXiv.org un documento que se publicará en el libro de proceedings del congreso "Solar Polarization Workshop #5" (editores S. Berdyugina, K.N. Nagendra, R. Ramelli). Se trata de un artículo de F. Paletou, R. Rezaei y L. Léger en el que se muestra un ejemplo de reducción y análisis de datos espectropolarimétricos en el campo de la física solar. Los datos que sirven para el ejemplo fueron tomados con el telescopio franco-italiano THéMIS, en el Observatorio del Teide, Tenerife. La novedad es que se utiliza para ello software completamente basado en Python y algunas de sus librerías científicas, como PyFITS, NumPy/Scipy y matplotlib.

Reproduzco aquí el resúmen del documento:

 

Most of the solar physicists use very expensive software for data reduction and visualization. We present hereafter a reliable freeware solution based on the Python language. This is made possible by the association of the latter with a small set of additional libraries developed in the scientific community. It provides then a very powerful and economical alternative to other interactive data languages. Although it can also be used for any kind of post-processing of data, we demonstrate the capabities of such a set of freeware tools using THeMIS observations of the second solar spectrum.

Algunos enlaces relacionados:


dic 03 2007

Medir sobre figuras con python

Tag: Herramientas y Programaciónahc @ 11:42

Actualización: Hay una nueva versión de Srip disponible en la página web del proyecto.

Cuando en una publicación -del tipo que sea- queremos mostrar un conjunto de datos numéricos que guardan cierta relación entre sí, la forma habitual de hacerlo es mediante una tabla o un gráfico. A menudo, el gráfico es preferible porque la distribución de los valores salta a la vista, y permite una comprensión cualitativa de la información representada mucho más rápida y eficiente que la lectura de valores en una tabla. Por otro lado, cualquier análisis cuantitativo requerirá disponer de los valores numéricos tabulados.

 

Representar los datos gráficamente, si se dispone de una tabla, es una tarea trivial y rápida si se utiliza un ordenador. El proceso inverso -el de recuperar los valores numéricos representados en una figura- suele ser tedioso e impreciso, porque en general el ordenador no puede hacerlo y es necesario que el usuario lea los valores manualmente, uno por uno, y a menudo con la dificultad añadida de escalas logarítmicas en los ejes o pocas marcas de referencia.
captura de sripPor desgracia, es muy habitual que necesitemos utilizar información numérica que sólo está disponible en forma de gráficos. En astrofísica, casos comunes son las curvas de transmisión de filtros, de eficiencia de detectores, espectros, cartas de localización de astros, diagramas color-magnitud, etc etc.

 

Para facilitar la tarea del investigador que necesita medir valores sobre un gráfico, he escrito un programa en python que, aunque aún está en fase de desarrollo, ya ha facilitado la vida a más de uno.

 

El mecanismo de funcionamiento es muy simple: se calibra la imagen que contiene el gráfico para establecer una correspondencia entre las coordenadas en píxeles de la imagen, y las unidades físicas representadas en los ejes del gráfico. A continuación, el usuario simplemente marca los puntos de la imagen que quiere medir y obtiene sus valores en las unidades físicas.

 

La calibración consiste en marcar una serie de puntos (bastan dos para la coordenada X y otros dos para la Y) para los que conocemos su valor en unidades físicas (por ejemplo, las marcas en los ejes del gráfico), e indicarle al programa cuál es este valor.

 

Si los píxeles p0 = (px0,py0) y p1 = (px1,py1) tienen valores físicos v0 = (vx0,vy0) y v1 = (vx1,vy1), entonces para cualquier otro pixel p la correspondencia será:

 

v = v0 + (p - p0)(v1-v0)/(p1-p0)

si la escala de los ejes es lineal, y:

 

v = v0 exp[(p - p0) ln(v1/v0)/(p1-p0)]

si es logarítmica.
captura de srip
Con la imagen calibrada, las mediciones se realizan simplemente haciendo click sobre la imagen, y las coordenadas físicas de los puntos marcados van apareciendo en una tabla. La incertidumbre en la medida es del equivalente a un píxel en las unidades físicas -si se tiene buen pulso al marcar- y puede reducirse haciendo zoom sobre la imagen.

 

El código, disponible aquí, está escrito en python y usa la librería Numeric para los cálculos y pyGTK para la interfaz gráfica. Requiere ImageMagik (convert) para rotar/escalar la imagen original, y puede utilizarse con prácticamente cualquier formato de imagen.


dic 02 2007

Ajustes 1-D con Python + Scipy

Tag: Herramientas y Programaciónfranwerst @ 12:30

En ocasiones es necesario ajustar una serie de puntos a una función dada, por ejemplo cuando queremos determinar parámetros físicos de una determinada línea espectral o cuando queremos modelar la transmisión de un determinado instrumento mediante un polinomio, etc...

Los ajustes unidimensionales en Python se pueden realizar de modo sencillo utilizando el módulo optimize del paquete matemático scipy. La función leastsq implementa una modificación del algoritmo Levenberg-Marquardt para minimizar la suma de los cuadrados de M ecuaciones con N incógnitas, a partir de una estimación inicial de la solución. Para obtener más información sobre el funcionamiento de leastsq con los parámetros de entrada o salida que permite se puede simplemente recurrir a la ayuda de Python:

>>> from scipy import optimize
>>> help(optimize.leastsq)

Aquí vamos a mostrar un ejemplo de su uso aplicándolo a un problema concreto como es el ajuste de una línea espectral a dos funciones: una gaussiana y a una lorentziana. El código correspondiente a este ejemplo se puede descargar desde aquí. Primero definimos las dos funciones de Python, una por cada perfil:


def gauss(params,x):
     """ Funcion gaussiana"""
     A=params[0]
     x0=params[1]
     sigma=params[2]
     g=A*exp(-0.5*power((x-x0)/sigma,2))
     return g

def lorentz(params,x):
     """ Funcion lorentziana"""
     A=params[0]
     x0=params[1]
     gamma=params[2]
     lorentz=1.0/(power(x-x0,2)+power(gamma,2))
     lorentz_norm=lorentz/max(lorentz)
     return A*lorentz_norm

Ambas fuciones dependen de tres parámetros que se introducen mediante la lista "params", y que representan la altura máxima de la curva (A), la abscisa correspondiente al máximo (x0) y una medida de la anchura (sigma en la gaussiana o gamma en la lorentziana). El vector x representa las abscisas del dominio de la función.

También necesitaremos otras dos funciones que dado el vector "y" con los puntos experimentales y un conjunto de parámetros "params" (que representan a una única curva) nos devuelva las diferencias entre los puntos y la función, es decir, calcula los residuos:

def res_gauss(params,y,x):
     """Residuos de la gaussiana"""
     return y-gauss(params,x)

def res_lorentz(params,y,x):
     """Residuos de la lorentziana"""
     return y-lorentz(params,x)

Ahora veamos cómo leastsq utiliza estas funciones para calcular el mejor ajuste. Pero primero vamos a crear una línea sintética que tenga una forma gaussiana con un cierto ruido. Definimos un vector de abscisas con 200 puntos entre x=-20 y x=20:

x=(arange(200)-100)/5.0

Creamos también un vector de ruido mediante el uso de la función random(). La intensidad del ruido será al máximo de 0.3, y para que el vector se mantenga en torno a cero le restamos 0.15 a cada punto.

random_vector=array([0.3*random()-0.15 for i in x])

Creamos la línea sintética problema con los parámetros A=-1, x0=2 y sigma=3 y le añadimos el ruido que hemos creado anteriomente.

linea=gauss([-1.0,2.0,3.0],x)+random_vector

Comencemos con los ajustes. Necesitamos unos parámetros iniciales con los que comenzar a iterar en la búsqueda de la solución óptima. Probemos con todos los parámetros iguales a 1. El algoritmo puede no converger adecuadamente si escogemos unos parámetros iniciales poco afortunados.

guess=[1.0,1.0,1.0]

Por fin realizamos los ajustes con leastsq:

ajuste_g=optimize.leastsq(res_gauss,guess,args=(linea,x))
ajuste_l=optimize.leastsq(res_lorentz,guess,args=(linea,x))

Finalmente mostramos los coeficientes resultantes de los ajustes:

print ajuste_g[0]
print ajuste_l[0]

Vemos que leastsq requiere como parámetros de entrada la función de python que calcula los residuos (esto es lo que vamos a minimizar), un conjunto inicial de parametros y después como argumentos las abscisas y las ordenadas correspondientes a la línea a ajustar. Como salida leastsq devuelve una lista cuya primera componente son los coeficientes de la solución. En nuestro ejemplo los coeficientes obtenidos son:

[-0.9927015   1.99923233 -3.01312331]   ->  para la gaussiana
[-1.08437082  2.02586623 -2.76829756]   ->  para la lorentziana

En la gráfica se pueden ver los resultados de los ajustes. En negro vemos la línea sintética que hemos creado, en rojo la gaussiana ajustada y en azul la lorentziana. Como es lógico el ajuste gaussiano es mejor ya que la línea es en realidad una gaussiana. Para verlo podemos, por ejemplo, echar un vistazo a la suma de los residuos en ambos casos:

print sum(res_gauss(ajuste_g[0],linea,x))
print sum(res_lorentz(ajuste_l[0],linea,x))
0.189616782797  -> para la gaussiana
5.68377589754   -> para la lorentziana

Algunos enlaces relacionados: