Python course

I'm not sure I'll follow with this blog...
I'm now involved in teaching Python to astronomers colleagues... and I set up another one, where I share the course I'm giving. If you are interested, follow it THERE.

calling shell commands from within python

That's quite easy to interact with the shell from within python, it only needs to import the subprocess module.
If the matter is just executing a command, use the call method:

>>>>subprocess.call('ls')
Applications    SharedLin       Volumes         etc             mach_kernel     sbin            var
Developer       SharedXP        bin             home            net             sw
Library         System          cores           iraf            opt             tmp
Network         Users           dev             lost+found      private         usr

If some arguments need to be used, you cannot put them in the same string as the command, need to use:

>>>>subprocess.call(['ls','-l'])
total 40733
drwxrwxr-x+ 134 root                admin      4556 28 jui 00:21 Applications
drwxrwxr-x@  18 root                admin       612 18 mar 20:03 Developer
[...]

or you'll need to execute the command in a shell:
>>>>subprocess.call('ls -l', shell=True)

If one want to deal with the output, use the Popen and communicate methods :

>>>>ls = subprocess.Popen(["ls"], stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
>>>>stdout, stderr = ls.communicate()
You can access the distinct values of the stdout by splitting them:
stdout.split()

Now, one can need to execute the command with a givenfile as input and redirect the output to another file.
That's possible using the Popen method:

input_file = open('model.in','r')
output_file = open('model.out','wb')
p = subprocess.Popen('cloudy.exe', stdout=output_file, stdin = input_file)
input_file.close()
output_file.close()
That's it!

References:
http://docs.python.org/library/subprocess.html

http://jimmyg.org/blog/2009/working-with-python-subprocess.html

.r and %run and Debugging

From IDL, it's easy to run a program at main level and stay after the execution with all the variables as they are defined by the program.
In (i)python, if you interactively import myprog it will execute it, but let you without any access to what has been set within the program.
There is a way to do the equivalent of .r in IDL, but you needs to run ipython.

ipython:
%run myprog
..... executing the myprog.py and give you the prompt back so you can explore the variables and objects to see if they are as they are supposed to be...

The %run (the % is needed, it's not a prompt) function is recompiling the code each time it is used, so changes are taken into account (it's not the case with import).

Nice tip for developpers...

BTW, there is also a debugger mode : from ipython, just type %pdb and then import or %run what you want. If something wrong happens, you will be let where it happened, with access to all the variables... like in IDL ;-) You have this activated by default using ipython -pdb.

More and more important informations  on ipython here (I think I must read these pages every week for the next 3 months!):
http://ipython.org/documentation.html

Broadcasting: adding arrays of different shapes

Something very different between python-numpy and IDL is the way they are both dealing with linear algebra in case of operations on different size or shape arrays.

In IDL, trying to operate on 2 vectors with different size reduces the oprtation to the lowest size:
IDL> a=[1,2,3]
IDL> b=[10,20,30,40]
IDL> print,a+b
       11       22       33
As we can see, the latest element is purely omitted.
IDL> c=[[1,2,3,4]]
IDL> d=[[10],[20],[30],[40]]
IDL> print,c
       1       2       3       4
IDL> print,d
      10
      20
      30
      40
The shape of the first operand is conserved:
IDL> print,c+d
      11      22      33      44
IDL> print,d+c
      11
      22
      33
      44
And we can loose part of the vector:
IDL> print,a+c
       2       4       6
IDL> print,a+d
      11      22      33


On the contrary, Python-numpy is adding some information, this is the broadcasting.

import numpy as np
a=np.array([1,2,3.])
b=np.array([10,20,30.,40])
a*b
# Traceback (most recent call last):
#  File "<stdin>", line 1, in <module>
#ValueError: shape mismatch: objects cannot be broadcast to a single shape


BUT:

a = a.reshape(3,1)
a
#array([[ 1.],
#      [ 2.],
#       [ 3.]])
b
#array([ 10.,  20.,  30.,  40.])

a*b
#array([[  10.,   20.,   30.,   40.],
#       [  20.,   40.,   60.,   80.],
#       [  30.,   60.,   90.,  120.]])

Another example where Python is guessing in which dimension the array must be extended:
c= a*b
c.shape
# (3, 4)
c*a
#array([[  10.,   20.,   30.,   40.],
#       [  40.,   80.,  120.,  160.],
#       [  90.,  180.,  270.,  360.]])
c*b
#array([[  100.,   400.,   900.,  1600.],
#       [  200.,   800.,  1800.,  3200.],
#       [  300.,  1200.,  2700.,  4800.]])
b2==np.array([10,20,30.])
#array([False, False, False], dtype=bool)
d=a*b2
d.shape
#(3, 3)
d
#array([[ 1.,  2.,  3.],
#       [ 2.,  4.,  6.],
3       [ 3.,  6.,  9.]])
d*a
#array([[  1.,   2.,   3.],
#       [  4.,   8.,  12.],
#       [  9.,  18.,  27.]])
d*b2
#array([[  1.,   4.,   9.],
#       [  2.,   8.,  18.],
#       [  3.,  12.,  27.]])

More details and nice figures in:

http://www.scipy.org/EricsBroadcastingDoc

Structure-like object (2): record arrays, access and views

Some complement to the previous message on structured arrays and record arrays:

import numpy as np
a = np.zeros((10,),dtype=[('name', str), ('ra', float), ('dec', float)])
a['ra'] = np.random.random_sample(10)*360
a['dec'] = np.random.random_sample(10)*180-90
tt = ((a['ra'] > 5.) & (abs(a['dec']) < 10.))
b = a[tt]
a.size
#10
tt.size
#10
tt.sum()
#1
b.size
#1

If no names age given to the different tags, it is set by default to f0,f1,...fN.
One can access the data without knowing the names of the tags:
a[a.dtype.names[2]]
is more complicated that IDL>a.(2), but at least it's possible ;-) And it as some advantage: you can change the names:
a.dtype.names = ('obs','ra','dec')

Be careful with subset, they are views:
b=a[1]
a['dec'][1]
# 0.0
b['dec']
# 0.0
b['dec'] = 2
a['dec'][1]
# 2.0
But it's not so easy to see it:
b['dec'] is a['dec'][1]
# False

Now we can turn the structured array a into a recarray:

a2 = a.view(np.recarray)

This add a new access mode for the data:
>>> a2.dec
array([ 32.61106119,  82.72958898, -18.46190884,  44.79729473,
       -54.65838972, -23.78818937,   3.56472044, -79.63061338,
        15.81108779,  73.37221597])
Be careful, this is a view, it means that the data are NOT duplicated, they are the same:
>>> a2.dec[1] = 2
>>>
>>> a2.dec[1]
2.0
>>>
>>> a[1]['dec']
2.0
>>>
>>> a['dec'][1]
2.0

This can slow down the access to a2 AND to a !!! So not so useful, or for small tables.

Refs:
http://docs.scipy.org/doc/numpy/user/basics.rec.html
http://www.scipy.org/Cookbook/Recarray

Structure-like object

I really like the structure variable in IDL. Especially the arrays of structure. It allows to search for elements using the where function and to extract a sub-structure matching a give condition.
I mean, if I have an dataset like this:
IDL> a = replicate({name:'',ra:0.0,dec:0.0},1000)
I can search for all the elements matching for exemple:
IDL> tt = where(a.ra gt 10. and abs(a.dec) gt 5.)
IDL> b = a[tt]


The same (more or less) can be made using numpy (imported as np):

a = np.zeros((1000,),dtype=[('name', str), ('ra', float), ('dec', float)])
tt = ((a['ra'] > 5.) & (abs(a['dec']) < 10.))
b = a[tt]

I don't really know if this is the best way. And also don't know how to access for example the second tag without naming it, like in IDL a.(1)...

I can also create an object:
class obs(object):
    def __init__(self,name='',ra=0.,dec=0.):
        self.name=name
        self.ra=ra
        self.dec=decAnd even define an array of objects:
colec = np.empty( (3,3), dtype=object)
And then put the objects in the colec:
colec[:,:] = obs()
BUT this will create a collection of 9 times the same object!!!
colec[0,0].ra = 5.5
colec[1,1].ra
>>>5.5
Some loop needed here. But the worst is that one will loose all the power of linear algebra from numpy.
So stay with the first approach for now.

Row- or Column major arrays and loops order.

After a too long period of silence, I'm coming back to learn Python. I just want to point out a little difference between IDL and Python in the order arrays are stored in the computer memory.
There is two ways arrays can be stored: row- or column major. It has a direct impact on the way one has to loop on the arrays. IDL is like Fortran (column major) and Python is like C (row major). It means that in Python, as you move linearly through the memory of an array, the second dimension (rigthmost) changes the fastest, while in IDL the first (leftmost) dimension changes the fastest.
Consequence on the loop order:

for i in range(0,5):
 for j in range(0,4):

              ... a[i,j] ...

list comprehension

There is a way to write loops and test on lists in Python which is very condensate: the so-called list comprehension.
An example found in the French forum: http://www.developpez.net/forums/d962510/autres-langages/python-zope/general-python/operations-listes-dictionnaires/.
A guy wanted to transform a dictionary like this:
d= {'Ei': (1,3,4,4,6) , 'id' : ('r','r','t','t','t')}
into the list of occurence of each id:

[('r', [1, 3]), ('t', [4, 4, 6])]

One of the answers was:

L=[(d['id'][i],d['Ei'][i]) for i in xrange(0,len(d['Ei']))] 

R=[(x,[y[1] for y in L if y[0]==x]) for x in set(d['id']) ] 

Quite efficient!

The first step is to build a dictionary with the correspondances between the id and the Ei, the L:

[('r', 1), ('r', 3), ('t', 4), ('t', 4), ('t', 6)]
This is done with something like this:

In [8]: L2 = []

In [9]: for i in xrange(0,len(d['Ei'])):
   ...:     L2.append((d['id'][i],d['Ei'][i]))
Which can effectively be compacted as:

L=[(d['id'][i],d['Ei'][i]) for i in xrange(0,len(d['Ei']))] 

The second step (R) is the list comprehension form of the following, where we are counting the occurences of each id. First find the set of the different and uniq values of id: set() gives the solution.

Then looping on these values and find the occurences for each values. Finally put all this into a dictionary. The expand form would be:

R2=[]

uniqid = set(d['id'])
for x in uniqid:
    for y in L:
        if y[0]==x:
            R2.append([x,[y[1]]])

What can also be obtained like this:

R=[(x,[y[1] for y in L if y[0]==x]) for x in set(d['id']) ] 

This is more compact, in some sense more elegant, and I think more efficient (quicker), but not really sure of that...

Save and Restore (2)

I guess I got it, not the same syntax as in IDL, but quite as fast in reading/writing:
def save(file,**kwargs):
    """
    Save the value of some data in a file.
    Usage: save('misdatos.pypic',a=a,b=b,test=test)
    """
    import cPickle
    f=open(file,"wb")
    cPickle.dump(kwargs,f,protocol=2)
    f.close
def restore(file):
    """
    Read data saved with save function.
    Usage: datos = restore('misdatos.pypic')
    """
    import cPickle
    f=open(file,"rb")
    result = cPickle.load(f)
    f.close
    return result
For example (notice I  use to import my stuff as CM):
CM.save('data3.pypic',data3=data3,a=a,b=b)
and
dd=CM.restore('data3.pypic')
take a small second to write and read.
Something interesting:
dd is a dictionary containing the 3 variables data3, a and b.
Now if I want to use data3, I can "extract" it:

data3 = dd['data3']

Nice, and very quick, as it's not a copy, but rather an object poiting to the same memory place:
In [15]: data3 is dd['data3']
Out[15]: True
In [16]: id(data3)
Out[16]: 4393756504
In [17]: id(dd['data3'])
Out[17]: 4393756504
Nice way to do the things ;-)

ADD:
Some people can find easier to save data using

save(file,"data")
others will prefere:
save(file,data=data)
You can use both with the following save function:

def save(file,*args,**kwargs):
  """
  Save the value of some data in a file.
  Usage: save('misdatos.pypic','a',b=b)
  """
  import cPickle
  f=open(file,"wb")
  dico = kwargs
  for name in args:
    dico[name] = eval(name)
  cPickle.dump(dico,f,protocol=2)
  f.close

Save and Restore

I use a lot the save/restore facilities in IDL. BTW, this is the lack of a complete save/restore set of tools that avoid me to use GDL in some projects. I was looking for an equivalent in Python, and realized that there are some few tools, the main problem now is choosing the right one ;-)
I found this page: http://kbyanc.blogspot.com/2007/07/python-serializer-benchmarks.html
and made a few tests with cPickle, which is part of Python (no need to install extra module).
It is actually very efficient when used with the protocol=2 option.
Marshal seems to be even faster, but cannot handle the rec.array objects (at least it doesn't work for me...).
Example of the use:
In [2]: import CMorisset as CM
In [3]: data3 = CM.ReadFortran('test3.dat','a10,1x,f6.2,1x,f6.2,1x,i2',['name', 'ra', 'dec','mag'])
In [4]: import cPickle
In [5]: cPickle.dump(data3, open("data3.pickle", "wb"),protocol=2)
In [6]: data3=cPickle.load(open("data3.pickle","rb"))

Reading the 27Mo of the test3.dat file take me 30 seconds with the Fortran format, and only 1 with the cPickle function! The main problems are that 1) one need to know the name of the store variable, and 2)  only one object can be saved at a time.
I think it can be bypassed using a dictionary containing the variables and the names.
My main issue now is to build the dicionnary from the arguments passed to a function, so that I could just have:

save(file='data',dat1,dat2,dat3)

Will se this latter...

Parameters in functions: take care!

I found a nice page with something really strange for me, I'll have to keep this in my mind because it's totally different from IDL.
The page: http://hetland.org/writing/instant-python.html
The "problem", as I rewrote it:
Let's define 2 functions:
def test(x):
    x=2
def test2(x):
    x[0]=2
Now, we'll call test and test2 with different variables:
In [3]: y=1
In [4]: test(y)
In [5]: y
Out[5]: 1
y is not changed. BUT:
In [7]: y=[1,2,3]
In [8]: test2(y)
In [9]: y
Out[9]: [2, 2, 3]
Now y[0] is changed!!! And even worst:
In [10]: y=[[4,5,6],[14,15,16],[24,25,26]]
In [11]: test2(y)
In [12]: y
Out[12]: [2, [14, 15, 16], [24, 25, 26]]
And the best for the end:

In [13]: y=[[4,5,6],[14,15,16],[24,25,26]]

In [15]: test2(y[0])
In [16]: y
Out[16]: [[2, 5, 6], [14, 15, 16], [24, 25, 26]]
So part of the table can be changed, but a single variable not: exactly the opposite of IDL...

My first module! Still on reading ascii formatted file.

I did my first program in Python! Here it is:

def ReadFortran(file,format,names,comment="#"):
    """
    Read a file using a Fortran-style format.
    Return a NumPy rec.array with each column named following the given names.
    Example: data = ReadFortran('test2.dat','a10,1x,f6.2,1x,f6.2,1x,i2',['name', 'ra', 'dec','mag'],comment="#")
    Morisset, IA-UNAM, Oct. 2010
    """
    import Scientific.IO.FortranFormat as FF
    import numpy.core.records as nprec
    FFformat = FF.FortranFormat(format)
    f=open(file,'r')
    rows=[]
    for line in f:
        if line[0] != comment:
            row = FF.FortranLine(line,FFformat)
            rows.append(row.data)
    f.close()
    return nprec.fromrecords(rows, names=names)

OK, it's not a very big one, but it took me a lot of time trying to avoid the list.append command. And I didn't found. But it seems that most of the execution time is on the FortranLine command.
It is very slower than the same in IDL: it reads a 1000000 lines file in some 45 seconds, while IDL take 5... The csv2rec takes 20 secs.
Perhaps one of these days I'll try to call a fortran routine to read the file...

ADD:
It seems that a more compact and pythonesk way of writing the loop is to change:
rows=[]
    for line in f:
        if line[0] != comment:
            row = FF.FortranLine(line,FFformat)
            rows.append(row.data)
into:
rows = [FF.FortranLine(line,FFformat).data for line in f if line[0] != comment]

The map function could also be used:

rows = map(lambda line:FF.FortranLine(line,FFformat).data,f)

BUT in this latest case we can't manage the comment parameter.

Got the tips from http://jaynes.colorado.edu/PythonIdioms.html

Reading formated ascii file a la fortran

I first learned to program in Fortran (after some introductions to Basic, ADA, Turbo pascal, in the early 80's). Then I meet IDL in 1994 (thanks to mi friend Philippe) and the life changed! Interactive + Data + Language was exactly what I needed. But as already said at the beginning of this blog, I now want to change for a free access open language.
But I feel very difficult this change, I'm like a baby learning to walk and talk... For example, I was looking since 2 weeks a way to read a simple formated ascii file, like I used to do in IDL.
The file is just:

     alpha 193.63  18.40 19
      beta 280.12   0.52 16
     gamma 206.59   0.06 17
     delta  23.74  17.92 19
       eta  18.10  10.07 19
and the IDL process is:

data = replicate({name:'',ra:0.0,dec:0.0,mag:0},n_lines)

openr,lun,/get_lun,file

readf,lun,data,format='(a10,1x,f6.2,1x,f6.2,1x,i2)'

et voilà!

The format string is using the Fortran convention, that is quite powerful in describing quite any fixed format. It seems that it's was not possible to do this in Python, 'till I found a module including this facility! Developped by french people from CNRS in Orleans, it is avilable here:
http://dirac.cnrs-orleans.fr/plone/software/scientificpython/

The part of the module I want is this one:

Module FortranFormat

Fortran-style formatted input/output

This module provides two classes that aid in reading and writing Fortran-formatted text files.

Examples:

Input::

   >>>s = '   59999'
   >>>format = FortranFormat('2I4')
   >>>line = FortranLine(s, format)
   >>>print line[0]
   >>>print line[1]

 prints::

   >>>5
   >>>9999


 Output::

   >>>format = FortranFormat('2D15.5')
   >>>line = FortranLine([3.1415926, 2.71828], format)
   >>>print str(line)

 prints::

   '3.14159D+00    2.71828D+00'

I used it to read the same file as previously in IDL
data=numpy.rec.array(['          ',0.,0.,0], names=['name', 'ra', 'dec','mag'])
from Scientific.IO import FortranFormat as FF
format = FF.FortranFormat('a10,1x,f6.2,1x,f6.2,1x,i2')
f=open('test1.dat','r')
for line in f:
#    data['name'],data['ra'],data['dec'],data['mag'] = FF.FortranLine(line,format)
    data.name,data.ra,data.dec,data.mag = FF.FortranLine(line,format)
    print data
The main problem is that I don't know how to have the whole array in the data variable. Anyway, the problem of reading fixed formatted ascii file is solved ;-)

Arrays and types: the rec.array type

Here are some example of unexpected behaviors, when coming from IDL:
In [4]: a
Out[4]: [1.0, 2.0, 3] #a is a list
In [5]: b
Out[5]: array([ 0.84147098,  0.90929743,  0.14112001])
In [6]: 2*b
Out[6]: array([ 1.68294197,  1.81859485,  0.28224002])
In [7]: 2*a
Out[7]: [1.0, 2.0, 3, 1.0, 2.0, 3] #this double the list, not the elements!!!
In [10]: c=array([1,2.,3])
In [11]: 2*c
Out[11]: array([ 2.,  4.,  6.]) #now it's OK
In [12]: d=[1,2.,'3'] #mixing int, float and string
In [13]: e=array(d)
In [14]: e
Out[14]: array(['1', '2.0', '3'], dtype='|S8') #converted to a single type: string

As you can see, the [] are not always defining arrays, sometime it needs explicitly the function array().
A list can contain different types, not an array.

Let's see now how the csv-reader can deal with different types in the same line (as when using structures in IDL):
The file to read is:
int  flt  str  int
1 1 "test1" 1
2 2.3 "tralala" 2
5 3.14 "" 3
6 1e3 "double " 4
7 1e79 "big one" 5

The reading process:

rec=csv2rec('test1.csv',delimiter=" ")

And the type is a rec.array, with mixing types:

In [43]: rec
Out[43]:
rec.array([(1, 1.0, 'test1', 1), (2, 2.2999999999999998, 'tralala', 2),
       (5, 3.1400000000000001, '', 3), (6, 1000.0, 'double ', 4),
       (7, 9.9999999999999997e+78, 'big one', 5)],
      dtype=[('int', '<i8'), ('flt', '<f8'), ('str', '|S7'), ('int_1', '<i8')])

Note that no problems with the value of 1 for the 2nd elements in the first data row: it's an integer, but as the element of the same column in the 2nd data row is 2.3, it is not wrongly considered that the type of this row is integer (contrary to IDL read_ascii(), which only consider the first data row to determine the data types).
Also note the headers, defined using the first row: int is used twice, the second time the name is converted to int_1. So cute!

the basic syntax (1)

Following the IDL cookbook, I first explore some syntax of Python.

I use the ipython, so the keyboards inputs are transcript with "In [n]" and the corresponding output with "Out [n]". To print something, just type it and press Enter:

In [4]: 'Hello mundo;-)'
Out[4]: 'Hello mundo;-)'

In [5]: a=1
In [6]: a
Out[6]: 1
In [8]: a?
Type:        int
Base Class:    <type 'int'>
String Form:    1
Namespace:    Interactive
Docstring:
    int(x[, base]) -> integer

.... 

In [35]: 3*exp(-12/(27.*log(1.3)))
Out[35]: 0.55135006225210048

In [36]: 3*exp(-12/(27.*log10([12.,14.,15])))
Out[36]: array([ 2.5086749 ,  2.53502116,  2.54592125])
 

Special characters:

;  is use to have 2 comands on the same line (and is NOT starting a comment!) 

a=5 ; b= 3

which actually can be done with:

a,b = 5,3

No matter spaces IF NOT AT THE BEGINNING OF THE LINE:

a , b = 5 , 3  

BTW this can be very powerful:

a, b, c = b, a+b, c+1  

Do you remember this one? Fibonnacci!

# is starting a comment 

a = b # a and b are pointing to different memory place

Parameters are separated by "," on the same line:

can = Canvas(f,width =250, height =250, bg ='ivory')

The blocks are defined by indentation. No begin, no end. In the following, the ....: are printed by ipython. The indentation is automatic in ipython, but not in python: you have to type at least one space to define the IF block. To finish a block, just empty lines (one in python, 2 in ipython...)
In [70]: a,b=5,6

In [71]: if a<b:
   ....:     print a+b
   ....:    
   ....:    
11

In [65]: def f(x):
   ....:     return x*2
   ....: 

In [67]: f(5)
Out[67]: 10

The same module in different libraries

One thing that I found nice in Python compared to IDL, is the amount of very complete and complex libraries that have been developed and are accessible so easily... And one thing I quickly found horrible is the amount os libraries!...

Just an example: I installed the Mayavi package (it means that I have to also install the vtk package, but hopefully all is already in the macport repositories...). and began to try playing with it.
And I found that sometimes some tutorial are talking about functions I don't have... For example, the points3d() function is not available. Well, I realized that it is in the mlab module and that actually I have... 6 mlab.py files in the directory where Python is looking for!!!
Just asking from a command line shell:

find /opt/local/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages mlab.py

give me:
./enthought/mayavi/mlab.py
./enthought/mayavi/tools/mlab.py
./enthought/tvtk/tools/mlab.py
./matplotlib/mlab.py
./numpy/numarray/mlab.py
./numpy/oldnumeric/mlab.py
So I have to be very careful with which module I download and in which order, to be sure that I have the mlab module I want.
Or the best is to call it explicitely and to check what is in each:

from enthought.mayavi import mlab as mlab_mayavi 

from enthought.tvtk.tools import mlab as mlab_tvtk
To obtain the list of the functions available from one module, just type for example mlab_mayavi.  and use the TAB to obtain the possible completions.
We can very quickly realize that the needed points3d() function is absent in all mlab modules except the mayavi one.
So be very careful with the modules, mlab is not always mlab...

Arrays bounds using NumPy

NumPy is load by default when calling ipython -pylab. If not, it's always possible to import all the module by from numpy import *  or import numpy as N (then you need to start any numpy function with N.).
To easy create arrays, just use:

a=r_[1:10]

It will create an array from 1 to... 9 !!! It is actually a kind of shortcut for:

a=arange(1,10)

which also provides an array of 9 elements, from 1 to 9.
The step can be specified:
a=r_[1:10:.1] or a=arange(1,10,.1).
In this case, the array latest element is 9.9, not 10.
It is also possible to give the number of elements, for this a "j" is added to the 3rd parameter:
b=r_[1:10:100j] or b=linspace(1,10,100).
BUT in these latest two equivalent cases, the latest element is actually 10!...
Another way to say this:
a=r_[1:10:.1] and b=r_[1:9.9:90j] are equivalent. Not obvious for me...
BTW, the same philosophy apply when selecting part of an array:
In [71]: a=r_[1:10]
In [72]: a
Out[72]: array([1, 2, 3, 4, 5, 6, 7, 8, 9])
In [73]: a[0:2]
Out[73]: array([1, 2])
I certainly would expect the answer to be a 3 elements sub-array. I'm sure this will be a big source of bug in my future Python codes ;-) 
I found something quite powefull: accessing the latest elements of an array:
In [77]: a[-1]
Out[77]: 9
In [78]: a[-5:-1]
Out[78]: array([5, 6, 7, 8])
In [79]: a[-5::]
Out[79]: array([5, 6, 7, 8, 9])

I also saw that the indexing of arrays can be quite sophisticated: http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html

Quick function definition: lambda and "closure"

There is 2 ways to define a function in Python. The classical:

def fitfunc(p,x):

         return p[0]*exp(-(x-p[1])**2/(2.0*p[2]**2))
and with the use of lambda:
fitfunc = lambda p, x: p[0]*exp(-(x-p[1])**2/(2.0*p[2]**2))

Both can be used in any place:
plot(x,fitfunc([1.,0.7,.3],x),'.')

From what I saw on the net, the compact lambda definition is more used when the function is small and used locally...

BTW, while Googling around this concept, I found this:

In [114]: def summer(nombre):
   .....:     def sumit(value):
   .....:         return value+nombre
   .....:     return sumit
   .....: 

In [115]: sum10 = summer(10)

In [116]: sum10(20)
Out[116]: 30

The summer function return a function, which depend on the value of the nombre variable. And the sum10 function keep in mind the value 10, even after leaving the summer function. This is what the computing science people name a "closure"...

filtering data: the where() function

I'm now trying to use my favorite IDL function in Python: the where() function. It is used to select the subscripts of an array where a condition is fulfilled.
As I'm using the pylab version of ipython, it comes with a where function that can match more or less the IDL one. But it seems that another method can also be used...
Let's do some examples: I want to select from a table of 1000000 lines and 10 columns the elements that have a value for the 3rd and 7th columns bigger than 0.5.


I'm first creating a 2D-table on which I will apply my filter.

a=random((1e6,10))

I first had to realize that the order of the subscripts are in the inverse order than in IDL: first rows, then columns.
So the elements of the 3rd columns are a[:,2], and the ones for the 7th columns are a[:,6].
I first try:

tt = where(a[:,2] > 0.5 and a[:,6] > 0.5)
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

OK, I'm not yet at the level of understanding what Python told me, but clearly it's not correct. I finally found that with parenthesis and & instead og and, things are going well:
tt = where((a[:,2] > 0.5) & (a[:,6] > 0.5))

tt is a so-called "tuple", which seems to be an array, it has a size:
In [90]: size(tt)
Out[90]: 250248

It's a credible value given the condition and the size of the input table.

I can now use this variable to extract the values from the initial table:

x=a[tt,2]

y=a[tt,6]

The main problem here is that the table I have in x and y are not correctly shaped:
In [112]: x.shape
Out[112]: (1, 250248)

If I try to plot this (using plot(x,y,'.')), it doesn't work (well, after some minutes waiting for a result I killed the plot!).
I can reshape the x and y tables using for example the transpose function:
plot(transpose(x),transpose(y),'.')
but the best is to transpose the filter before using it:
tt2 = transpose(where((a[:,2] > 0.5) & (a[:,6] > 0.5)))
plot(a[tt2,2],a[tt2,6],'.')
is working fine. BTW, tt2 is not anymore a tuple, it's not an array of integers, so for example:

In [146]: tt2.size
Out[146]: 250248

Another way of filtering the data is to generate a table of booleans:
tt3 = (a[:,2] > 0.5) & (a[:,6] > 0.5)
This is an array:
In [147]: tt3.size
Out[147]: 1000000

In [148]: tt3.dtype
Out[148]: dtype('bool')
Contrary to IDL, it can directly be used in a table:
In [149]: a[tt3,2].size
Out[149]: 250248
In [150]: a[tt3,2].shape
Out[150]: (250248,)

Fine, the plot works also with this:
In [154]: plot(a[tt3,2],a[tt3,6],'.')

I tried both methods (where and boolean) on big table, but didn't saw any real difference on time execution. If some readers could tell me which one is really the more "python-way"...

Reading an ascii file

I used to use the read_ascii() function in IDL (BTW I use my own version, which return an array of structures and not a structure containing arrays...), and it seems there is not an equivalent tool... or I didn't found it yet!
Nevertheless, I tried to make the job simpler by reading a Comma Separated Variable file (csv), as it can be output by any spreadsheet or downloaded from the net (e.g. Vizier tables).
I first found a module able to read such file: csv (!).
This is the first solution I used:

import csv
from scipy import *
rr = csv.DictReader(open('DIGEDA1.1.csv'))
x=[]
y=[]
for row in rr:
    x.append(row['5007'])
    y.append(row['3727'])
plot(x,y)

The problem is that I didn't succeed to apply simple filter to the data: I want to plot the log of the values, so I have to filter out the negative values... It seems that the format used for the x and y variables (namely "list") is not really the best one to be used with where.
 
Finally, I found another function to read the file, which outputs something very similar to a structure.
It is part of the pylab distribution (matplolib), so no need to import it as I use the pylab argument when calling ipython. So here my second (and successful) try to read, filter and plot my data:

rec=csv2rec('DIGEDA1.1.csv')
tt = where(logical_and(rec['3727']>0,rec['5007']>0))
rec2=rec[tt]
x=rec2['5007']
y=rec2['3727']
plot(log(y),log(x),'o')

The rec variable contains the whole file, in such a way very similar to IDL structure. A very nice thing is that the names of the different columns are directly taken from the first line of my file, which contains this information.

It is finally quite easy (well, it took me 3 hours to find all this in the Google jungle... it's why I'm doing this blog!).