Sherpa modeling and fitting package on Ubuntu

Sherpa is a set of fitting and modelling routines used by the people from the Chandra X-Ray Observatory. It is a very robust software with its own interface, although recently they developed a python implementation which is very easy to use. I guess x-ray astronomers must know what they are doing, don't you?

You can find the python package here.
When I tried to install it on Ubuntu 10.04 or 10.10 I faced a really nasty compilation problem with the file Simplex.cc. To fix it you must edit the following file before compiling:

sherpa-4.2.2/sherpa/optmethods/Simplex.cc

Just add this line to the beginning

#include <stdio.h> 

If you have all dependencies installed the compilation should run just fine. In the following weeks I will post some examples with Sherpa.

Python: how to start learning?

Many (most?) people that want to start learning Python are confused about where to start. So many options! Motivated by this, I list in this post the references that I used to learn Python (and object-oriented programming as well), which can serve as a starting point for other people. I had scientists in mind when I wrote this post.

Beginner material

Official Python documentation

Learned the basic syntax and capabilities of the language with the official Python tutorial. You can download all of this as PDF files. I suggest this for people with previous programming experience. For absolute beginners, have a look at the Think Python book below.

Computação científica com Python, Fabrício Ferrari (portuguese only)

Introductory lecture about Python, its syntax and science applications. It shows what Python is capable of for data analysis and plotting. Inspiring. The audio is also available for download as a MP3 file.

Using Python for interactive data analysis, Greenfield & Jedrzejewski (STSCI)

Tutorial on using Python for data analysis! How to on how to replace IDL/Matlab with Python, essentially. Includes: plotting, FITS files, signal processing.

Introduction to Programming using Python, Programming Course for Biologists at the Pasteur Institute, Schuerer et al.

I learned object-oriented programming using this material. Very clear and "application-oriented" approach. You don't need to be a biologist to understand this.

Think Python: How to think like a computer scientist

Longer introduction for people with no previous extensive programming experience.

Quick reference

Python data analysis reference card

Reference card: Python commands equivalent to IDL/Matlab ones

Migrating from IDL/Matlab to Python.

IPython quick reference card

If you are going to do serious stuff with Python, I suggest using the enhanced interactive Python terminal IPython.

Longer introductory books

Learning Python, Mark Lutz

A primer on scientific programming with Python, Hans Petter Langtangen

Longer reference books

Python essential reference, David Beazley

Collection of Python-related links

New STSCI Python package (2.10)

Those who updated Ubuntu from version 9.10 (Karmic Koala) to 10.04 (Lucid Lynx) might have noticed that Python 2.5 is not available in the new default repository. This is a big issue for the users such as the users of the Space Telescope Science Institute Python packages (PyRAF, MultiDrizzle, etc...) or some older google engine apps.

Fortunately a new version of the stsci_python is out today! With full support to Python 2.6 (3.0 will take longer).

http://www.stsci.edu/resources/software_hardware/pyraf/stsci_python/current/download

Update (21/06)

Some may have the following error when importing the iraf modules using

from pyraf import iraf
  File "", line 38
    iraf.set(as = 'host$as/')
              ^
SyntaxError: invalid syntax
/usr/local/lib/python2.6/dist-packages/pyraf/irafimport.py:54: RuntimeWarning: Parent module 'pyraf' not found while handling absolute import
  return _originalImport(name, globals, locals, fromlist, level)

To fix it the trailing slash on the /iraf/iraf/unix/hlib/zzsetenv file on line 38 should be removed. Everything runs smoothly.

Select fits files based on header keyword

The task of selecting a list of fits files based on some header keyword value is much easier to do with the IRAF task hselect. Although this example is very useful for understanding how to use pyfits and the sys module for passing parameters via command line (like bash our Perl's @ARGV).

#!/usr/bin/python

import sys
import pyfits
import re

for file in sys.argv[1:]:  #exclude the first element of sys.argv
     hdulist = pyfits.open(file)  #open fits
     key = hdulist[1].header['OBJECT'] #reads a keyword from the fits header
     if (re.search('standard',key)):  #search the keyword value for some string
             print file

To execute just run on the command line:

# python script.py file1 file2 ...

or

# chmod +x script.py
# ./script.py file1 file2 ...

All bash's wildcards (*, ?, etc) are accepted.

Using f2py

The purpose of the F2PY --Fortran to Python interface generator-- is to provide connection between Python and Fortran languages. Here I will show a simple example of how to use this tool and compare the its efficiency against a pure python code.

First, lets write our fortran routine:

        Subroutine test(a,b,n,c)

        Implicit none
        Integer::i
        Integer, intent(in)::n
        Real, dimension(n), intent(in)::a,b
        Real, dimension(n), intent(out)::c
!f2py   depend(n)::a,b,c       
!f2py   intent(in)::a,b
!f2py   intent(out)::c

        Do i=1,n
          c(i) = sqrt(a(i)**2.0 + b(i)**2.0)
        EndDo

        End subroutine test

This code simply calculate the distance to the origin of a bunch of points. The input is two vectors with identical dimension n. Note that you must specify what goes in and what goes out to the f2py parser, this is done by using the a comment line followed by f2py:

!f2py   depend(n)::a,b,c       
!f2py   intent(in)::a,b
!f2py   intent(out)::c

Now lets compile and link it to python using:

f2py -m func -h func.pyf func.f --overwrite-signature
f2py --fcompiler=gfortran -c func.pyf func.f

Unfortunately f2py does not easily support the Intel Fortran Compiler (ifort), although gfortran does a good job on most of the cases.

It ridiculously easy to call this routine inside python. The script bellow call the fortran routine and then compares its execution time with an identical pure python routine plus the same done using Numpy.

#!/usr/bin/python
import numpy as np
import time
import func

n = 10000
#generates a normal distribution with n points
x = np.random.randn(n)
y = np.random.randn(n)

# Using the fortran routine
t1 = time.time()
f = func.test(x,y,n)
t2 = time.time()

t = (t2-t1)*1000

print 'Fortran runtime (ms): ',t

# Using pure python
t1 = time.time()
r = np.zeros(n)
for i in range(n):
 r[i] = np.sqrt(x[i]**2 + y[i]**2)
t2 = time.time()

t = (t2-t1)*1000

print 'Python runtime (ms): ',t

# Using pure python
t1 = time.time()
r = np.sqrt(x**2 + y**2)
t2 = time.time()
t = (t2-t1)*1000

print 'Python + numpy runtime (ms): ',t

If everything goes OK the result should be something like this:

Fortran runtime (ms):  0.150918960571
Python runtime (ms):  76.4260292053
Python + numpy runtime (ms):  0.485181808472

The advantages of using f2py are obvious, the runtime of the pure python code is almost 500 times greater! When using numpy to handle the arrays we get a much better execution time, but still slower when compared to fortran.

I hope this post was useful to the people thinking of migrating to python in the near future.

Mais Python?

Só pra mostrar que eu faço alguma coisa da vida, ai vai mais um post com um tutorial em python.

Agora vou mostrar como interagir com gráficos. Isso na verdade é bastante complicado, requer um entendimento forte de orientação a objetos(OO). Nesse exemplo eu tentei evitar OO, justamente porque não estou tão familiarizado com isso.

A proposta é: dado um conjunto de pontos como selecionar somente aqueles dentro de um polígono definido pelo usuário interagindo com o gráfico? A solução é o seguinte a seguir. Como de costume, precisa o NumPy, MatPlotLib e PyLab.

O programa funciona da seguinte maneira. Enquanto você vai com o mouse no ponto onde quer definir um dos vértices do polígono e aperta alguma tecla do teclado que não seja a tecla 'q', um ponto será marcado. Ao apertar 'q' você sai do programa e os pontos dentro do polígono são marcados pela cor vermelha.

#!/usr/bin/python
from pylab import *
import numpy as np
import matplotlib.nxutils as nx

poly=np.zeros([100,2]) #limits the number of vertices to 100 (how to ovecome this?) 
i = int(0)             #not very python-ish
x = np.random.rand(1000)
y = np.random.rand(1000)
f = np.array([x,y])
f = np.transpose(f)

def click(event):
    global poly,f,i,g
    tb = get_current_fig_manager().toolbar
    if event.key != 'q' and event.inaxes and tb.mode == '':
      poly[i,0],poly[i,1] = event.xdata,event.ydata
      plot(poly[0:i+1,0],poly[0:i+1,1],'rs-')
      draw()
      i += 1
    else:
      poly[i,0],poly[i,1] = poly[0,0],poly[0,1]
      poly = np.resize(poly,(i+1,2))
      inside = nx.points_inside_poly(f[:,0:2], poly)
      g = f[nonzero(inside),]   #this adds an extra dimension to the array (why?!)
      plot(g[:,:,0],g[:,:,1],'r.',markersize=6)
      return g

fig = figure()
ax = fig.add_subplot(111,autoscale_on=True)
ax.plot(f[:,0],f[:,1],'k.',markersize=4)

#connects python with the display
cid = connect('key_press_event', click)
show()
disconnect(cid)

Felizmente, nesse exemplo não precisa baixar nada. Basta copiar, colar e rodar. Ah sim, o highlight da sintaxe python foi feito nesse site.

Ia escrever outro post sobre um resolvedor de simulação 3 corpos em python, mas já ta no meu site e eu cansei do blogspot, queria mesmo é o wordpress. Nesse exemplo tem até instruções pra criar uma animação.

PS.: caiu fora do post as linhas muito longas, não to com saco pra arrumar agora

Python \o/

Então, nas ultimas semanas estive envolvido com Python bastante. O motivo disso é que o Dark Energy Survey (DES) adota o python como sua linguagem de frontend. Como eu faço parte do grupo de populações estelares do DES fui obrigado a lidar com Python em diversos níveis. Então vou explicar uma das coisas que tive que desenvolver nesse projeto.

O objetivo desse tuturial é gerar uma carta do céu e identificar certas estrelas nesta carta. O propósito é o planejamento de uma observação astronômica, nesse caso do aglomerado globular NGC 2298. Com essa carta o observador vai saber se o telescópio está apontando na direção que deveria estar apontando. O campo é uma simulação do campo de visão do instrumento MOSAIC2 do telescópio Blanco de 4m.

Bom, vamos ao Python... Para rodar esse programa você vai precisar dos módulos NumPy e PyLab, ambos disponíveis através do repositório Debian/Ubuntu e do módulo APLpy que deve ser instalado manualmente. Além disso para rodar você precisa baixar e descompactar os arquivos específicos para esse problema aqui.

O produto final, se tudo der certo, deve ser algo do tipo:

#!/usr/bin/python
from pylab import *
import aplpy
import numpy as np

Alphabet = map(chr, range(65, 91))   # Capitalized alphabet 

fig = figure(figsize=(10,10))
gc = aplpy.FITSFigure('field.fits',figure=fig, subplot=[0.15,0.1,0.8,0.8])
gc.show_grayscale()
gc.show_grid()
gc.set_grid_alpha(0.2)

# read and plot the GSC stars
f = np.loadtxt('gsc.dat')
ra,dec = f[:,0],f[:,1]
gc.show_markers(ra,dec,edgecolor='red',marker='o',s=50)

#this convert from world coordinates to pixel (canvas) coordinates
newra,newdec = gc.world2pixel(ra+0.025,dec)

#place a tag near each marked point in the image
for i in range(f.shape[0]):
  text(newra[i],newdec[i],Alphabet[i],fontsize=16,color='red')

#save as png
gc.save('findingchart.png')