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 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: