Welcome to my lessons

Bo Zhang (NAOC, mailto:bozhang@nao.cas.cn) will have a few lessons on python.

These are very useful knowledge, skills and code styles when you use python to process astronomical data.
All materials can be found on my github page.
jupyter notebook (formerly named ipython notebook) is recommeded to use

These lectures are organized as below:

install python
basic syntax
numerical computing
scientific computing
plotting
astronomical data processing
high performance computing
version control

+ - * /



In [1]:

    
%pylab inline
a = np.arange(10)
print a









    



Populating the interactive namespace from numpy and matplotlib
[0 1 2 3 4 5 6 7 8 9]



In [2]:

    
print type(a)
print len(a)
print size(a)
print shape(a)
print a.shape
print a.reshape((-1, 5)).shape









    



<type 'numpy.ndarray'>
10
10
(10,)
(10,)
(2, 5)



In [3]:

    
b = a + np.random.randn(10)



In [4]:

    
plt.plot(b)









    Out[4]:





[<matplotlib.lines.Line2D at 0x7fd0c6507510>]



In [5]:

    
print 'a+b'
print a+b
print 'a-b'
print a-b
print 'a*b'
print a*b
print 'a/b'
print a/b
print 'a^b'
print a**b









    



a+b
[  0.23937492   1.88796844   4.83892141   6.3690742    8.66241255
  11.1787515   11.79377214  15.48062711  14.83953616  20.37811288]
a-b
[-0.23937492  0.11203156 -0.83892141 -0.3690742  -0.66241255 -1.1787515
  0.20622786 -1.48062711  1.16046384 -2.37811288]
a*b
[   0.            0.88796844    5.67784281   10.1072226    18.6496502
   30.89375752   34.76263284   59.36438979   54.71628929  102.40301594]
a/b
[ 0.          1.12616615  0.70449291  0.89045234  0.85792494  0.80922497
  1.03559475  0.82541066  1.16966996  0.79099233]
a^b
[  0.00000000e+00   1.00000000e+00   7.15484942e+00   4.05001759e+01
   6.41286443e+02   2.08334870e+04   3.22426727e+04   1.46879556e+07
   1.50216149e+06   7.20242906e+10]

numpy tools



In [6]:

    
%pylab inline
np.__version__









    



Populating the interactive namespace from numpy and matplotlib






    Out[6]:





'1.11.1'



In [7]:

    
print 1/2
print 1./2



In [8]:

    
np.argsort(np.random.randn(100,))









    Out[8]:





array([13, 82, 57, 29, 95, 22, 18, 66, 89, 34, 85, 20,  5, 46, 42, 73,  4,
       38,  9, 19, 76, 70, 77, 47,  7, 61, 26,  8, 49, 52, 48, 37, 50, 98,
        0, 12, 54, 14, 91, 83, 28, 88, 78, 45, 44, 24, 39, 72, 41, 75, 32,
       43, 96, 53, 56, 68,  3, 58, 23, 33, 59, 35, 81, 16, 55, 87, 64, 69,
       97, 86,  1, 40, 84, 93, 30,  6, 17, 21, 27,  2, 25, 71, 92, 99, 63,
       94, 15, 74, 11, 51, 36, 79, 60, 10, 80, 31, 65, 62, 90, 67])

Vectorization()



In [9]:

    
a = np.arange(100)
b = np.zeros(a.shape)
a_list = [_ for _ in a]
print a_list









    



[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99]



In [10]:

    
print a_list*3









    



[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99]



In [11]:

    
%%time
for i in xrange(len(a)):
    b[i] = a[i]*10









    



CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 51 µs



In [12]:

    
%%time
c = a*10









    



CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 49.1 µs

matrix manipulation

python use element-wise multiplication



In [13]:

    
a = np.random.randn(5, 10)
b = np.random.randn(10, 5)



In [14]:

    
%%time
print np.dot(a, b)









    



[[ 1.25388074 -0.33695197  0.97878567  1.43996735 -1.04775042]
 [ 3.23627877  1.95764194  1.7910861  -1.09184855  0.02772932]
 [-3.00779501 -1.47261213  0.45334396  0.39091005 -1.95127891]
 [-0.9040374   1.12676008 -2.25067679 -1.58093382 -3.78433705]
 [-3.10202347  0.3307634   1.18192323  1.97328469 -0.0785904 ]]
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 753 µs



In [15]:

    
%%time
print a*b









    



---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-15-f3ad7429ce2f> in <module>()
----> 1 get_ipython().run_cell_magic(u'time', u'', u'print a*b')

/usr/local/lib/python2.7/dist-packages/IPython/core/interactiveshell.pyc in run_cell_magic(self, magic_name, line, cell)
   2101             magic_arg_s = self.var_expand(line, stack_depth)
   2102             with self.builtin_trap:
-> 2103                 result = fn(magic_arg_s, cell)
   2104             return result
   2105 

<decorator-gen-59> in time(self, line, cell, local_ns)

/usr/local/lib/python2.7/dist-packages/IPython/core/magic.pyc in <lambda>(f, *a, **k)
    186     # but it's overkill for just that one bit of state.
    187     def magic_deco(arg):
--> 188         call = lambda f, *a, **k: f(*a, **k)
    189 
    190         if callable(arg):

/usr/local/lib/python2.7/dist-packages/IPython/core/magics/execution.pyc in time(self, line, cell, local_ns)
   1174         else:
   1175             st = clock2()
-> 1176             exec(code, glob, local_ns)
   1177             end = clock2()
   1178             out = None

<timed exec> in <module>()

ValueError: operands could not be broadcast together with shapes (5,10) (10,5)



In [ ]:

HOMEWORK

try to use np.sin, np.cos, np.dot to build a function which can translate a sightline from Equatorial coordinates to Galactic coordinates. Usage:
```
l, b = func(ra=194., dec=30.5)
print l, b
```
Given a mask image (True/False), generate a list of pixels met requirements:
1. The pixel should be on the image
2. The distance from the generated pixels to at least 1 of the True pixel is within a given upper limit, say d0
3. Usage:

>>> def random_dots(mask_image, d_upper=20., n_dots=100):
>>>    pass

>>> mask_ = np.loadtxt('./data/mask_image/mask_image.dat')>0
>>> d_upper = 20.
>>> n_dots = 100
>>> x_rand, y_rand = random_dots(mask_image, d_upper=d_upper, n_dots=n_dots)
>>> plt.plot(x_rand, y_rand, 'rx')



In [16]:

    
# import numpy as np
len1 = 700
len2 = 900
mask_image_data = np.random.randn(len1, len2)
mesh_x, mesh_y  = meshgrid(np.arange(len2), np.arange(len1))

mask = mask_image_data>4
print '# of center pixels: ', np.sum(mask)
plt.imshow(mask, cmap='jet')
plt.plot(mesh_x[mask], mesh_y[mask], 'o', ms=30, mfc='w', markeredgecolor='w')
plt.plot(mesh_x[mask], mesh_y[mask], 'kx')
plt.xlim(0, len2)
plt.ylim(0, len1)
plt.title('scheme of valid random dots area');









    



# of center pixels:  26



In [17]:

    
np.savetxt('./data/mask_image/mask_image.dat', mask, fmt='%d')



In [ ]: