There is currently a large, and growing, community of Python programmers and scientists, and there is a plethora of free technical code available. Doing science with Python appears to be a good way to contribute to free scientific programming, rather than catering to (and putting money in the pockets of) companies like Research Systems and Mathworks. Furthermore, as Python is older than many commercial packages, these programs appear to have borrowed considerably from Python syntax, so that in many ways Python programming is similar to IDL or Matlab programming, whereas the causality is the other way around. However, there are some packages (e.g. matplotlib) which have been specifically designed to use syntax similar to familiar commercial products.
Python is a very complicated language, but one in which simple things can be done fairly quickly. It is probably the fullest, most complex (does that mean 'best'?), object oriented language I have seen, and is more similar to common Lisp than it is to C++ or Java (or even Smalltalk or Eiffel). One could spend the rest of one's life fully learning Python and doing only very esoteric, ivory tower, 'pure computing' things with it. C is a much more 'down to Earth', 'nuts and bolts' language. However, that isn't really a problem, for reasons discussed below. To begin learning Python I recommend the following:
1: ~ > python Python 2.3.5 (#1, Jan 13 2006, 20:13:11) [GCC 4.0.1 (Apple Computer, Inc. build 5250)] on darwin Type "help", "copyright", "credits" or "license" for more information. >>>
for a Python 'Easter egg', type the following and hit return:
>>> import this
I won't tell you what this does, as it should be an exercise for the reader. Cognoscenti should detect the (somewhat snide) reference to Perl, and another to Python's creator. Here's another funny from the end of the Programming Python book and The Empire Strikes Back. To quit out of Python, type control-D at the triple prompt.
Unfortunately, version 2.3.5 isn't the one you want for scientific computing, and the Mac installation comes with no additional packages, so the next thing you have to do is download and install some software. And then download and install some software. And then download, compile, and install some more software. And then... Getting a usable scientific Python system is not a trivial task. I've been working on this for about a week now, and I'm not done yet, although most of the things I want to do finally seem to work. As this is Unix software, the {download, compile, install} pocess is fairly straightforward, although you do have to do things in specific order --which only becomes apparent when something you spent an hour downloading and installing doesn't work because you didn't download and install something else first, or you downloaded and installed some things out of order. Also, it sometimes seems that to install a peice of Unix software requires you to first install two other peices of software, which requires you to first etc etc...
I seem to have a fairly stable system at present, so I hope that my past week's experiences can help you to get up and running a bit faster than I did. However, what would be really great would be if someone took all the separate packages, libraries, external programs, etc..., and packed them up as a single static stand-alone application that worked right out of the box. Until then, here's my best guess at what you will need to do to install a nominal Python scientific computing environment. I've also created CDs of all the various software and manuals I've acquired, and of the final Python 'site-packages' and /usr/local directories that I created, which should at least save you some time looking for things on the web.
(NOTE: You may want to skip a bit here and scroll down the page to look at the examples and screenshots first.)
3: ~ > python2.4 Python 2.4.4 (#1, Oct 18 2006, 10:34:39) [GCC 4.0.1 (Apple Computer, Inc. build 5341)] on darwin Type "help", "copyright", "credits" or "license" for more information. >>>
Entering just 'python' will probably still get you Python 2.3.5. To fix this, cd to /usr/bin and list all the pythons:
4: ~ > cd /usr/bin 5: /usr/bin > ls -al python* lrwxr-xr-x 1 root wheel 24 Apr 17 10:24 python -> python2.3 lrwxr-xr-x 1 root wheel 72 Feb 19 2006 python2.3 -> ../../System/Library/Frameworks/Python.framework/Versions/2.3/bin/python lrwxr-xr-x 1 root wheel 25 Apr 18 15:13 pythonw -> pythonw2.3 -rwxr-xr-x 1 root wheel 34216 Jan 13 2006 pythonw2.3
You need to replace the python and pythonw links with ones to python2.4. To do this, you will probably have to be Super User. If so, type 'su' and the root password. Then remove the python and pythonw links with 'rm python' and 'rm pythonw' and make then point to python2.4 with 'ln -s /usr/local/bin/python2.4 python' and 'ln -s /usr/local/bin/pythonw2.4 pythonw'. If all went well, your files should now look like:
lrwxr-xr-x 1 root wheel 24 Apr 17 10:24 python -> /usr/local/bin/python2.4 lrwxr-xr-x 1 root wheel 72 Feb 19 2006 python2.3 -> ../../System/Library/Frameworks/Python.framework/Versions/2.3/bin/python lrwxr-xr-x 1 root wheel 25 Apr 18 15:13 pythonw -> /usr/local/bin/pythonw2.4 -rwxr-xr-x 1 root wheel 34216 Jan 13 2006 pythonw2.3
and you can then exit root by pressing control-D. Now saying 'python' in a
shell (and in a '#!/usr/bin/python' or '#!/usr/bin/env python' line in a
program file) should get you Python 2.4.4.
libfreetype: Freetype fonts system http://www.freetype.org/ libjpeg: JPEG image support http://www.ijg.org/ libpng: PNG image support (also in /usr/lib/libz) http://www.libpng.org/pub/png/
In general, for astronomical computing, you should also have the following libraries:
libcfitsio: Read/write FITS files http://heasarc.gsfc.nasa.gov/docs/software/fitsio/ libfftw3: Fourier transforms http://www.fftw.org/
Get these libraries as .tar.gz source files. They (and most Unix software) are
all built in similar ways:
(cd to directory with configure script) ./configure make make install (as root)
Usually these libraries are installed in /usr/local/lib, with headers in
/usr/local/include, and any programs in /usr/local/bin. Some configuration
scripts may require a location prefix flag: e.g. './configure
--prefix=/usr/local'. See the README or INSTALL files in each module for
specific instructions. If the libraries and headers don't appear in /usr/local
as desired, you may have to copy them manually from wherever they were built.
Usually, you need to be root to copy to /usr/local. You may also need to run
ranlib on a moved library to update its timestamp, if the linker complains.
backend : TkAgg numerix : numpy # numpy, Numeric or numarray interactive : True # see http://matplotlib.sourceforge.net/interactive.html
This sets the GUI to use Tk, arrays as ndarrays, and interactive mode active.
This assumes that you have compiled matplotlib with the default settings in
setup.py:
BUILD_AGG = 1 ... BUILD_TKAGG = 'auto'
matplotlib can also be built to use gtk or wx widgets, if you have those
installed. Python comes with Tk.
numerix : numarray # numpy, Numeric or numarray
I then set a shell environment variable for PyFITS:
> setenv NUMERIX numarray
I then launch SAOImageDS9 under X-Windows, and run the following Python commands:
> python
Python 2.4.4 (#1, Oct 18 2006, 10:34:39)
[GCC 4.0.1 (Apple Computer, Inc. build 5341)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import pyfits
>>> from numarray import *
>>> import numdisplay
>>> im = pyfits.getdata('c1230.en.map.fits')
>>> im.__class__
<class 'numarray.numarraycore.NumArray'>
>>> numdisplay.display(im)
Image displayed with Z1: 0.0 Z2: 1000.0
The FITS file is imported as a NumArray object, and is sent to DS9 for display. I can also display the file in a matplotlib/pylab Tk window:
>>> from pylab import * >>> imshow(im)
I can also display the file in a Cocoa window (although this is currently a C program launched from the shell):
>>> import os
>>> os.system('~/Projects/FITS/XC/plotfits c1230.en.map.fits')
To use the Cocoa display as a Python function, I would have to package the C program as a Python extension. However, Cocoa windows and graphics offer significant advantages over either Tk or X-Windows, but can only be run under OS X.
To summarize, there are at least 3 ways to display images from Python:
FITS arrays can also be manipulated and plotted in other Tk windows. Some array syntax will be demonstrated below, but in general it is similar to either IDL or Matlab, although more powerful than either:
>>> im.shape (700, 700) >>> plot(im[:, 350])
>>> plot(im)
> python Python 2.4.4 (#1, Oct 18 2006, 10:34:39) [GCC 4.0.1 (Apple Computer, Inc. build 5341)] on darwin Type "help", "copyright", "credits" or "license" for more information. >>> from numpy import * >>> a = range(5); a [0, 1, 2, 3, 4] >>> a = [range(5), range(5, 10), range(10, 15)]; a [[0, 1, 2, 3, 4], [5, 6, 7, 8, 9], [10, 11, 12, 13, 14]]
Create an array out of ranges, get its shape, then create a new real array from a shaped range:
>>> a = array([range(5), range(5, 10), range(10, 15)]); a
array([[ 0, 1, 2, 3, 4],
[ 5, 6, 7, 8, 9],
[10, 11, 12, 13, 14]])
>>> a.__class__
<type 'numpy.ndarray'>
>>> a.ndim, a.shape
(2, (3, 5))
>>> a = arange(15).reshape(3, 5) * 1.0; a
array([[ 0., 1., 2., 3., 4.],
[ 5., 6., 7., 8., 9.],
[ 10., 11., 12., 13., 14.]])
Reset the value of a row and column:
>>> a[:, 2] = -a[:, 2]; a
array([[ 0., 1., -2., 3., 4.],
[ 5., 6., -7., 8., 9.],
[ 10., 11., -12., 13., 14.]])
>>> a[1, :] = -a[1, :]; a
array([[ 0., 1., -2., 3., 4.],
[ -5., -6., 7., -8., -9.],
[ 10., 11., -12., 13., 14.]])
Try to take the real sqrt(), and then the complex sqrt():
>>> sqrt(a)
Warning: invalid value encountered in sqrt
Warning: invalid value encountered in reduce
Warning: invalid value encountered in reduce
array([[ 0. , 1. , nan, 1.73205081, 2. ],
[ nan, nan, 2.64575131, nan, nan],
[ 3.16227766, 3.31662479, nan, 3.60555128, 3.74165739]])
>>> sqrt(a + 0j)
array([[ 0. +0.j , 1. +0.j ,
0. +1.41421356j, 1.73205081+0.j , 2. +0.j ],
[ 0. +2.23606798j, 0. +2.44948974j,
2.64575131+0.j , 0. +2.82842712j, 0. +3.j ],
[ 3.16227766+0.j , 3.31662479+0.j ,
0. +3.46410162j, 3.60555128+0.j , 3.74165739+0.j ]])
Get some stats:
>>> print a.min(), a.max(), a.mean(), a.sum(), a.std() -12.0 14.0 1.4 21.0 8.1059648819
Find some values satisfying conditions:
>>> where(a > 0, 1, 0)
array([[0, 1, 0, 1, 1],
[0, 0, 1, 0, 0],
[1, 1, 0, 1, 1]])
>>> a > 0
array([[False, True, False, True, True],
[False, False, True, False, False],
[ True, True, False, True, True]], dtype=bool)
>>> a[a > 0]
array([ 1., 3., 4., 7., 10., 11., 13., 14.])
Can do the same on a 1d version of arrays:
>>> a.ravel()
array([ 0., 1., -2., 3., 4., -5., -6., 7., -8., -9., 10.,
11., -12., 13., 14.])
>>> a.ravel() > 0
array([False, True, False, True, True, False, False, True, False,
False, True, True, False, True, True], dtype=bool)
>>> (a.ravel() > 0).nonzero()
(array([ 1, 3, 4, 7, 10, 11, 13, 14]),)
>>> a.ravel()[(a.ravel() > 0).nonzero()]
array([ 1., 3., 4., 7., 10., 11., 13., 14.])
>>> a.ravel()[(a.ravel() > 0).nonzero()] = 999; a
array([[ 0., 999., -2., 999., 999.],
[ -5., -6., 999., -8., -9.],
[ 999., 999., -12., 999., 999.]])
Compute the FFT of an array:
>>> from numpy.fft import *
>>> a = 1.0 * array([1, -1] * 4 * 8).reshape(8, 8); a
array([[ 1., -1., 1., -1., 1., -1., 1., -1.],
[ 1., -1., 1., -1., 1., -1., 1., -1.],
[ 1., -1., 1., -1., 1., -1., 1., -1.],
[ 1., -1., 1., -1., 1., -1., 1., -1.],
[ 1., -1., 1., -1., 1., -1., 1., -1.],
[ 1., -1., 1., -1., 1., -1., 1., -1.],
[ 1., -1., 1., -1., 1., -1., 1., -1.],
[ 1., -1., 1., -1., 1., -1., 1., -1.]])
>>> print abs(fft(a))
[[ 0. 0. 0. 0. 8. 0. 0. 0.]
[ 0. 0. 0. 0. 8. 0. 0. 0.]
[ 0. 0. 0. 0. 8. 0. 0. 0.]
[ 0. 0. 0. 0. 8. 0. 0. 0.]
[ 0. 0. 0. 0. 8. 0. 0. 0.]
[ 0. 0. 0. 0. 8. 0. 0. 0.]
[ 0. 0. 0. 0. 8. 0. 0. 0.]
[ 0. 0. 0. 0. 8. 0. 0. 0.]]
>>> print abs(fftshift(fft(a)))
[[ 8. 0. 0. 0. 0. 0. 0. 0.]
[ 8. 0. 0. 0. 0. 0. 0. 0.]
[ 8. 0. 0. 0. 0. 0. 0. 0.]
[ 8. 0. 0. 0. 0. 0. 0. 0.]
[ 8. 0. 0. 0. 0. 0. 0. 0.]
[ 8. 0. 0. 0. 0. 0. 0. 0.]
[ 8. 0. 0. 0. 0. 0. 0. 0.]
[ 8. 0. 0. 0. 0. 0. 0. 0.]]
The NumPy array object and packages contain dozens of additional functions for performing statistics, array manipulations, convolutions, polynomials, etc...
>>> from numpy.random import *
>>> a = matrix(rand(4, 4)); a
matrix([[ 0.81634207, 0.1337322 , 0.14083475, 0.21813705],
[ 0.56612496, 0.8580282 , 0.30047178, 0.3911366 ],
[ 0.88954812, 0.02283352, 0.86269536, 0.05555847],
[ 0.19919564, 0.0890195 , 0.53704153, 0.32616056]])
>>> a.I
matrix([[ 1.17141069, -0.11853493, 0.27883583, -0.68879087],
[-1.09956122, 1.43630936, 0.32853997, -1.04301787],
[-1.28867443, 0.11722304, 0.98403588, 0.55367178],
[ 1.70656493, -0.51263568, -1.88023184, 2.85966003]])
>>> print a.I * a
[[ 1.00000000e+00 2.77555756e-17 1.11022302e-16 5.55111512e-17]
[ 1.11022302e-16 1.00000000e+00 1.11022302e-16 1.11022302e-16]
[ -1.52655666e-16 0.00000000e+00 1.00000000e+00 0.00000000e+00]
[ -1.11022302e-16 -5.55111512e-17 -4.44089210e-16 1.00000000e+00]]
5: ~/Projects/Gauss-Jordan > date ; gj686 1000 0 ; date Sun Apr 22 17:13:00 MDT 2007 Sun Apr 22 17:13:45 MDT 2007
An equivalent Python function might be:
>>> def foo(n):
... os.system('date')
... a = matrix(rand(n, n))
... ai = a.I
... b = a * ai
... c = ai * a
... os.system('date')
... return c
...
>>> c = foo(1000)
Sun Apr 22 17:14:00 MDT 2007
Sun Apr 22 17:14:06 MDT 2007
>>> c
matrix([[ 1.00000000e+00, -7.49400542e-14, -1.73194792e-14, ...,
-1.36366612e-13, -9.56942858e-14, -5.10216869e-14],
[ 1.03028697e-13, 1.00000000e+00, 1.77635684e-15, ...,
5.63334102e-14, 2.99066327e-14, 8.91994811e-14],
[ 3.21964677e-15, 3.87849475e-14, 1.00000000e+00, ...,
1.02633180e-13, 6.93282237e-14, 6.65396557e-14],
...,
[ 2.07056594e-14, 1.44328993e-15, 1.35724765e-14, ...,
1.00000000e+00, 9.30366895e-14, -4.82253126e-15],
[ 6.59194921e-15, -2.49800181e-15, 2.21697660e-15, ...,
-1.39905448e-15, 1.00000000e+00, 2.49800181e-16],
[ -3.49997809e-14, -2.10942375e-14, -3.20229954e-14, ...,
-8.14452672e-14, -7.89646126e-14, 1.00000000e+00]])
>>> print c.dtype, c.shape
float64 (1000, 1000)
>>> c = foo(2000)
Sun Apr 22 17:26:00 MDT 2007
Sun Apr 22 17:26:44 MDT 2007
>>> c
matrix([[ 1.00000000e+00, -1.51732793e-13, -5.23192600e-14, ...,
-6.43009951e-14, -2.83165852e-13, -2.89716168e-14],
[ -2.05835349e-13, 1.00000000e+00, -3.13832293e-13, ...,
-2.08277839e-13, -9.95731275e-14, 2.25167107e-15],
[ 3.22519789e-14, 1.42719170e-13, 1.00000000e+00, ...,
2.22835639e-13, -4.53699578e-14, -1.32887625e-13],
...,
[ -1.85962357e-15, 1.60288449e-15, 7.93462518e-15, ...,
1.00000000e+00, 2.08166817e-16, -3.89271948e-15],
[ 1.57235336e-14, 1.14075416e-14, 1.83004653e-14, ...,
-2.17187379e-15, 1.00000000e+00, 9.55659163e-15],
[ 3.62210262e-15, -1.08940634e-15, -7.79584730e-15, ...,
1.76733628e-14, 1.25663369e-14, 1.00000000e+00]])
>>> print c.dtype, c.shape
float64 (2000, 2000)
Of course this isn't a direct comparison of Python's actual speed, since the inversion is certainly being done in compiled C code, and probably uses a method better than GJ. But, as a first pass, it suggests that Python is probably 'fast enough' for calculations. However, as additional information, consider the following:
>>> def foo(n):
... os.system('date')
... x = fft(rand(n, n))
... os.system('date')
... return x
...
>>> foo(1000)
Sun Apr 22 17:34:31 MDT 2007
Sun Apr 22 17:34:31 MDT 2007
array([[ 4.96001627e+02 +0.j , -5.96842134e+00 +3.5645091j ,
-1.26021109e+00 -9.10648446j, ..., 1.32986754e-01 +0.98455753j,
-1.26021109e+00 +9.10648446j, -5.96842134e+00 -3.5645091j ],
[ 5.12303163e+02 +0.j , 9.62648276e+00 +2.0764284j ,
-1.38441143e+01 -6.20285508j, ..., -1.51497810e+00 +7.30754387j,
-1.38441143e+01 +6.20285508j, 9.62648276e+00 -2.0764284j ],
[ 4.97417526e+02 +0.j , -9.06651415e-02 +8.20669437j,
5.74951458e+00 -4.5737919j , ..., -1.54866395e+01 -2.73620996j,
5.74951458e+00 +4.5737919j , -9.06651415e-02 -8.20669437j],
...,
[ 5.13899215e+02 +0.j , 1.53185048e+00 +1.18847739j,
-5.24489950e+00 -0.7996017j , ..., 1.29652505e+00 +1.90884849j,
-5.24489950e+00 +0.7996017j , 1.53185048e+00 -1.18847739j],
[ 5.03255483e+02 +0.j , 4.98178367e+00 +0.33379411j,
-4.45263312e+00 -0.24665999j, ..., -6.23934729e+00 -6.99924003j,
-4.45263312e+00 +0.24665999j, 4.98178367e+00 -0.33379411j],
[ 4.98274067e+02 +0.j , 3.75069495e+00+16.3833213j ,
-6.77387368e+00 -9.41446747j, ..., -8.00904375e+00 -1.64761378j,
-6.77387368e+00 +9.41446747j, 3.75069495e+00-16.3833213j ]])
>>> y = foo(1000)
Sun Apr 22 17:34:40 MDT 2007
Sun Apr 22 17:34:41 MDT 2007
>>> y = foo(5000)
Sun Apr 22 17:34:50 MDT 2007
Sun Apr 22 17:35:02 MDT 2007
>>> y = foo(10000)
Sun Apr 22 17:35:11 MDT 2007
Python(340) malloc: *** vm_allocate(size=1600000000) failed (error code=3)
Python(340) malloc: *** error: can't allocate region
Python(340) malloc: *** set a breakpoint in szone_error to debug
Traceback (most recent call last):
File "", line 1, in ?
File "", line 3, in foo
File "/Library/Frameworks/Python.framework/Versions/2.4/lib/python2.4/site-packages/numpy/fft/fftpack.py", line 89, in fft
return _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftf, _fft_cache)
File "/Library/Frameworks/Python.framework/Versions/2.4/lib/python2.4/site-packages/numpy/fft/fftpack.py", line 64, in _raw_fft
r = work_function(a, wsave)
MemoryError
Python/Tk also appears to be able to handle the complexity necessary to create something like Scientist's Workbench:
It seems clear that Python is capable of performing the calculations of either FITSRegister or FITSFlow, but it is not clear if Tk or another Python window/widget set is good enough to create the kind of GUI that can be made with Interface Builder and the Cocoa Application Kit classes, either from Objective C or from Python. It might be worthwhile to spend a little time investigating this question, as it should be fairly quickly obvious whether there is much chance of this or not. And if it is not possible with Tk, perhaps one of the other window/widget sets (e.g. wx, gtk) might be more appropriate. One thing I can say at present, however, is that I would not want to try to create a program like FITSRegister or FITSFlow (or even DS9) using only the X11 library. I also have some doubts that a good GUI can be created entirely programmatically --relying on geometry managers, packers, and the like--, rather than via an interactive (and manual) layout program such as IB, although Java and Swing seem to be able to produce reasonable results (on the Mac at least). I suppose that is something else to consider, whether we might be better off to move to Java, rather than Python, especially if interplatform operability is a goal.
I'm going to spend another week, at least, working with Python, and should have more info then. In the face of doubt, I would suggest trying to write some non-trivial program in Python as a test case. I would pick something that would be useful to the Venus project, and something that involved graphics and UI, but which would be useful even without a GUI. I would pick something that could be accomplished during the next month (i.e. May), as that would not represent too great of an investment. If that goes well, we could continue using Python. If not, then we can either fall back to C, or try something else.
İSky Coyote 2007