#######################
# Features of IPython #
#######################

ipython
In [10]: x = [1, 2] # create a short list
In [11]: x?
Type: list
Base Class: <type 'list'>
String Form: [1, 2]
Namespace?: Interactive
Length: 2
Docstring:
list() -> new list
list(sequence) -> new list initialized from sequence's items
In [12]: x.<TAB>
x.__add__ x.__iadd__ x.__setattr__
x.__class__ x.__imul__ x.__setitem__
x.__contains__ x.__init__ x.__setslice__
x.__delattr__ x.__iter__ x.__str__
x.__delitem__ x.__le__ x.append
x.__delslice__ x.__len__ x.count
x.__doc__ x.__lt__ x.extend
x.__eq__ x.__mul__ x.index
x.__ge__ x.__ne__ x.insert
x.__getattribute__ x.__new__ x.pop
x.__getitem__ x.__reduce__ x.remove
x.__getslice__ x.__reduce_ex__ x.reverse
x.__gt__ x.__repr__ x.sort
x.__hash__ x.__rmul__
In [13]: x.append?
Type: builtin_function_or_method
Base Class: <type 'builtin_function_or_method'>
String Form: <built-in method append of list object at 0x15bcff0>
Namespace: Interactive
Docstring:
L.append(object) -- append object to end
In [14]: !ls
newfile.fits seq.fits pix.fits
In [15]: import pyfits
In [16]: pyfits.info 'pix.fits' # instead of pyfits.info('pix.fits')
------->pyfits.info('pix.fits')
Filename: pix.fits
In [17]:,pyfits.info pix.fits
------->pyfits.info('pix.fits')
Filename: pix.fits
No. Name Type Cards Dimensions Format
0 PRIMARY PrimaryHDU 71 (512, 512) Int16

In [18]:,print hello
-------->print('hello')
hello
In[19]: print 1
1

In[20]: print 2
2

In[21]: print 3
3

In[22]: print 4
4

In [23]: exec In[19:21]+In[22]
1
2
4
In[25]: 6+9
15
In[26]: print _
15

####################################
# A little on Python Introspection #
####################################

>>> x = [1,2] # a simple python list
>>> print x.__doc__
list() -> new list
list(sequence) -> new list initialized from sequence's items
This shows the docstring for the list creator function.
>>> print pyfits.__doc__ # print the pyfits module docstring
A module for reading ....
....
>>> print pyfits.getdata.__doc__
...

###################################
# Python looping and conditionals #
###################################

>>> pets = ['dog','cat','gerbil','canary','goldfish','rock']
>>> for pet in pets: print 'I have a pet', pet
I have a pet dog
I have a pet cat
I have a pet gerbil
I have a pet canary
I have a pet goldfish
I have a pet rock
>>> for i in range(3): print i
0
1
2
>>> x = 0
>>> if x==0:
... print "x equals 0"
...
x equals 0
>>> if x==1:
... print "x equals 1"
... elif x==0:
... print "x equals 0"
... else:
... print "x equals something other than 1 or 0"
...
x equals 0
>>> if '': print 'True' # nothing will be printed
...
>>> if '0': print 'True'
...
True
>>> x = arange(9)
>>> if x: print 'True'
[...]
RuntimeError: An array doesn't make sense as a truth value. Use
sometrue(a) or all true(a)

##################
# More on PyFITS #
##################

>>> hdr = pyfits.getheader('pix.fits')
>>> hdr.add_history('This is a test history card')
>>> hdr.add_comment('A test comment card',after='date')
>>> hdr.add_blank('',before='comment')
# To delete a keyword just use the del statement
>>> del hdr['date']
>>> hdus = pyfits.open('fuse.fits') # returns an HDUList object
>>> tabhdu = hdus[1] # the table extension
>>> hdr = tabhdu.header # getting the header from the HDU
>>> tab = tabhdu.data # and the table

>>> hdus = pyfits.open('acs.fits')
>>> sci = hdus['sci',1] # EXTNAME=SCI, EXTVER=1

>>> hdus.writeto('temp.fits') # write current form to new file
>>> hdus.flush() # update opened file with current form

>>> ff = pyfits.open('nicmos.fits')
>>> sum = 0. # plan is to sum all extension images
>>> for hdu in ff[1:]:
... if hdu.header['extname'].startswith('SCI'):
... sum = sum + hdu.data + 2.**15 # now data is read in
... hdu.data = None # done with it, delete it from memory
>>> ff.close()

>>> ff = pyfits.open('pix.fits') # returns an HDUList object
>>> hdu = ff[0] # Get first (and only) HDU
>>> sum = 0L # for keeping a running sum of all image values
>>> # now read data 64 lines at a time (although small, it illustrates
the principle)
>>> for i in range(8):
... subim = hdu.section[i*64:(i+1)*64,:]
... sum = sum + subim.sum()
...
>>> print sum
>>> ff.close()

>>> hdu = pyfits.PrimaryHDU() # create an HDU from scratch
>>> hdu.header.update('P.I.','Hubble') # Try creating an illegal keyword
ValueError: Illegal keyword name 'P.I.'
# force into header an illegal keyword
>>> card = pyfits.Card().fromstring("P.I. = 'Hubble'")
>>> hdu.header.ascardlist().append(card)
>>> hdu.writeto('pi.fits') # try writing to file
Output verification result:
HDU 0:
Card 4:
Unfixable error: Illegal keyword name 'P.I'
[...]
raise VerifyError
VerifyError
# now force it to file
>>> h.writeto('pi.fits', output_verify='ignore')
>>> hdus = pyfits.open('pi.fits')
>>> print hdus[0].header.ascardlist()
[...]
P.I. = 'Hubble'
>>> hdus[0].header['p.i.']
'Hubble'
>>> hdus.verify()
[repeat of previous error message]

##############################################################
# Intermission for 1 powerpoint slide, time for some popcorn #
##############################################################


#################################
# Some numarray module examples #
#################################

>>> import numarray.random_array as ra
>>> ra.random((3,3)) # trivial example

>>> y, x = indices((11,11))
>>> y -= 5 # example of "augmented operator, subtract 5 from y and
place result in y
>>> x -= 5
>>> im = 1000*exp(-(x**2+y**2)/2**2) # noiseless gaussian
>>> nim = ra.poisson(im)
>>> print nim

>>> from numarray.fft import fft
>>> x = sin(arange(1024)/20.) + ra.normal(0,1,shape=(1024,))
>>> # sine wave with gaussian noise
>>> from pylab import *
>>> plot(abs(fft(x)))

##################
# Avoiding loops #
##################
OR
#######################
# Stupid Array Tricks #
#######################

# Using mask and index arrays
##############################

>>> im = pyfits.getdata('pix.fits')
>>> # Using the Fortran or C approach
>>> sum = 0L
>>> for j in range(512):
... for i in range(512):
... sum = sum + im[j, i] # remember that index order!!
print sum
28394234L
>>> Using the array approach
>>> print im.sum()
28394234L

>>> x = arange(10)
>>> ind = where(x > 5)
>>> x[ind] # this works
array([6, 7, 8, 9])
>>> ind.shape # this doesn't
[...]
AttributeError: 'tuple' object has no attribute 'shape'
>>> ind
(array([6, 7, 8, 9],)
>>> ind[0].shape # this does
(4,)
>>> len(ind) # this is effectively the number of dimensions
1
>>> len(ind[0]) # this is the number of points and is likely what was
desired
4

>>> cim = im.copy() # copy for future examples
>>> bigvalues = im[im>100] # use of mask array to index
>>> im[im>100] = 100
>>> im[where(im>100)] = 100 # same effect, but slower

>>> indices = where(cim > 300) # where function returns tuple of index
arrays
>>> indices
(array([ 31, 31, 32, ..., 505, 506, 506]),
array([319, 320, 318, ..., 208, 206, 207]))
>>> cim[indices[0], indices[1]] = 100 # the explicit approach of using
# the component index arrays
>>> cim[indices] = 100 # the "magic" approach, exactly equivalent to
the explicit form

# examples
##########

>>> # linear interpolation
>>> x = arange(10.)* 2
>>> y = x**2
>>> import numarray.random_array as ra
>>> xx = ra.uniform(0., 18., shape=(100,)) # 100 random numbers between
0 and 18
>>> xind = searchsorted(x, xx)-1 # indices where x is next value below
that in xx
>>> xfract = (xx - x[xind])/(x[xind+1]-x[xind]) # fractional distance
between x points
>>> yy = y[xind] + xfract*(y[xind+1]-y[xind]) # the interpolated values
for xx


>>> # coin tossing simulation, odds of HH before TT given prob of heads
= p
>>> p = .4
>>> n = 1000000 # initial ensemble size
>>> n2heads = 0; n2tails = 0 # totals for each case
>>> import numarray.random_array as ra
>>> prev = ra.random((n,)) < p # true if heads, false if not
>>> while n>0:
... next = ra.random((n,)) < p
... notdone = (next != prev)
... n2heads += sum(1.*(prev & next)) # careful about sums!
... n2tails += sum(1.*logical_not(prev | next))
... prev = next[notdone]
... n = len(prev)
... print n
...
480311
239661
[...]
1
1
0
>>> n2heads, n2tails
(336967.0, 663033.0)

>>> # calculate average radial profile
>>> y, x = indices((512,512)) # first determine radii of all pixels
>>> r = sqrt((x-257.)**2+(y-258)**2
>>> ind = argsort(r.flat) # get sorted indices (could use sort
# if not needed to arange image values too
>>> sr = r.flat[ind] # sorted radii
>>> sim = im.flat[ind] # image values sorted by radii
>>> ri = sr.astype(Int16) # integer part of radii (bin size = 1)

>>> deltar = ri[1:] - ri[-1] # assume all radii represented
# (more work if not)
>>> rind = where(deltar)[0] # location of changed radius
>>> nr = rind[1:] - rind[:-1] # number in radius bin
>>> sim = sim*1. # turn into double to avoid integer overflow
# (and single precision rounding)
>>> csim = cumsum(sim) # cumulative sum to figure out sums for each
radii bin
>>> tbin = csim[rind[1:]] - csim[rind[:-1]] # sum for image values in
radius bins
>>> radialprofile = tbin/nr # the answer
>>> semilogy(radialprofile)

>>> # solve Laplace's equation for 500x500 region
>>> # boundary conditions: outer edge = 0
>>> # central 100x100 = 100
>>> grid = zeros((500,500),Float32)
>>> prevgrid = grid.copy()
>>> grid[200:300,200:300] = 100.
>>> done = False
>>> i = 0
>>> while not done:
... i +=1
... prevgrid[:,:] = grid # just reuse arrays
... # new value is average of 4 neighboring values
... grid[1:-1,1:-1] = 0.25*(
... grid[:-2,1:-1]
... + grid[2:, 1:-1]
... + grid[1:-1,:-2]
... + grid[1:-1,2:])
... grid[200:300,200:300] = 100.
... diffmax = abs(grid - prevgrid).max()
... print i, diffmax
... if diffmax < 0.1: done = True
1 25.0
2 12.5
[...]
243 0.0996131896973

>>> # find nearest neighbors
>>> x = array([1.1, 1.8, 7.3, 3.4])
>>> y = array([2.3, 9.3, 1.5, 5.7])
>>> deltax = x - reshape(x,(len(x),1) # computes all possible
combinations
>>> deltay = y - reshape(y,(len(y),1) # could also use subtract.outer()
>>> dist = sqrt(deltax**2+deltay**2)
>>> dist = dist + identity(len(x))*dist.max() # eliminate self matching
>>> # dist is the matrix of distances from one coordinate to any other
>>> print argmin(dist) # the closest points corresponding to each
coordinate
[3 3 3 1]

>>> # ways of avoiding new array creation
>>> im = 2*im # replaces the array im used to refer to
# with a new array
>>> im[:,:] = 2*im # writes back into the array im
>>> im *= 2 # also reuses im array
# (but be careful of type downcasts!)
>>> multiply(im, 1, im) # use of optional output array