Wednesday, August 20, 2014

Introduction to Numpy

2014 Python lecture. Part II


The introduction to Numpy can be seen here:
https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/intro_numpy.ipynb

The topics that are presented are:

  • The Array class
    • create an array
    • 1D, 2D 3D arrays
    • creating array from scratch
    • arrays share memory (views)
  • random generator
  • timing a command
  • slicing arrays
  • assignments
  • using masks
  • the where function
  • some operations with arrays
  • broadcasting
  • calling scripts
  • structured arrays and record arrays
  • NaN other ANSI values.
Any comments are welcome.
Chris.Morisset a t Gmail.com

Wednesday, August 13, 2014

Python: Basics

2014 Python lecture. Part I


The introduction to Python I'm giving at IA-UNAM is accessible here:

https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/intro_Python.ipynb

I will modify this notebook during the lecture (August 2014), so reload it to have the latest version.

The topics of this first lecture are:
  • Using python as a calculator
  • assignments
  • comments
  • types
  • complex numbers
  • booleans
  • printing strings
  • strings
  • Tuples, lists and dictionaries
  • Blocks
  • List and dictionary comprehension
  • Functions, procedures
  • Scripting
  • Importing libraries
If you want to have an interactive session with this lecture using the ipython notebook facilities, follow the link above and download the ipynb file (download button at the right top of the web page). Save the file in a directory from where you execute the following (you must have a recent version of ipython installed):
ipython notebook
It should open a new tab in your web browser, with the list of ipynb files in the directory. Click on the one you want, will open a new tab similar to the first one, but this one is executed on YOUR computer, it means you are able to interact with the commands. You can change the commands, and execute a cell by SHIFT-ENTER. You can add comments in new cells, and save the result.

Any comments are welcome.

Thursday, August 7, 2014

Brief introduction to Python

2014 Python lecture. Part 0


Back to the Python lecture, I want to share here the very quick introduction I gave before starting to play with python: https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Intro_1.pdf
You may want to install python from Ureka from this site: http://ssb.stsci.edu/ureka/

Wednesday, March 19, 2014

Using ipython Notebook to teach scientific python

A very good and efficient way to teach python and python related tools, is to use ipython Notebook: http://ipython.org/notebook.html

As example of this use, the following link is a collection of lectures on python, numpy, scipy, matplotlib, use of Fortran from python, etc:
https://github.com/jrjohansson/scientific-python-lectures
Enjoy them.

Friday, June 28, 2013

Installing python and nice lectures on scientific Python.

Since more than one year without any message!... And the new one is almost nothing from me, just links to good pages.

Two easy ways to install python+ipython+numpy+matplotlib+scipy:

Anaconda from continuum (but only 64bit version for OSX, which can be a problem for MySQLdb):
Ask for the academic licence if you can.
http://continuum.io/
Once installed, you will need to add the anaconda/bin directory to your PATH and the anaconda directory to your PYTHONPATH.

Canopy from Entought:
https://www.enthought.com/products/canopy/
Once installed, you must setup the virtual environment by adding to your .tcshrc:
setenv VIRTUAL_ENV /Users/YOURNAME/Library/Enthought/Canopy_32bit/User
setenv PATH  $VIRTUAL_ENV/bin:$PATH

!!! Warning !!! When using this virtual environment, don't install using pip with the --user option. Install directly, as for example:
pip install pyfits

UPDATE: Another (better) package comes from STSCI, it's UREKA: http://ssb.stsci.edu/ureka/

Here follows good pages to learn python. Those pages are made using Notebook, which is very efficient to show/share python programs.

Here are the links:

http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-0-Scientific-Computing-with-Python.ipynb
http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-1-Introduction-to-Python-Programming.ipynb
http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-2-Numpy.ipynb
http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-3-Scipy.ipynb
http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-4-Matplotlib.ipynb
http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-6B-HPC.ipynb

Saturday, May 5, 2012

Playing with arrays: slicing, sorting, filtering, where function, etc.

You can also read the more recent post on Numpy here: http://python-astro.blogspot.mx/2014/08/introduction-to-numpy.html


Most of the data we have to manage are in form of arrays. That's why it's quite always necessary to import the numpy library to easily deal with tables, arrays, images, cubes, etc...
In this post I will discuss some of the methods and functions one may have to use in this context.

In Python, one can use lists, tuples and dictionaries to put different elements together. They even can contain elements of different types. But nothing better than numpy arrays to really "play" with the data, extract subsets, combine them using arithmetic operations.

We already discuss a little some array commands in a previous post: read-ascii-file-cont. Here we are more complete, but not exhaustive, as numpy is a whole world... Have a look at the reference guide: http://docs.scipy.org/doc/numpy/reference/ and other references at the end of this post.

Create arrays

a = np.array([1, 2, 3, 7, 5])
a
array([1, 2, 3, 7, 5])

numpy's arrays cannot contain values of different types, the values are transformed if necessary, from integer into real, and from real into string:
b = np.array([1, 2, 3, 4.])
b
array([ 1.,  2.,  3.,  4.])
c = np.array([1, 2, 3., '4'])
c
array(['1', '2', '3', '4'],
      dtype='|S1')

shapes and sizes of the arrays are obtained using method and functions:
print c.shape
(4,)
print len(c)
4

Different ways to create arrays:
a = np.arange(1, 10, 0.5)
print a
[ 1.   1.5  2.   2.5  3.   3.5  4.   4.5  5.   5.5  6.   6.5  7.   7.5  8.
  8.5  9.   9.5]

TAKE CARE, THE LATEST ELEMENT IS NOT WHAT YOU MAY EXPECT...
It's because the first element in Python is indexed by 0. The simple use of np.arange is :
a = np.arange(10)
print a
[0 1 2 3 4 5 6 7 8 9]
a starts at 0, and has 10 elements, that's ok.

Easy to create linearly spaced arrays:
c = np.linspace(0, 1, 4)
print c
[ 0.          0.33333333  0.66666667  1.        ]
or log spaced:

c = np.logspace(0, 1, 4)
c
array([  1.        ,   2.15443469,   4.64158883,  10.        ])
which is actually the same as :
c = 10**np.logspace(0, 1, 4)

To create arrays of 0. or 1.:
a = np.zeros(10)
b = np.ones((5, 4))
Here b is a 2D array.
b.shape
(5, 4)

You can use a list (or a tuple, or an array) and replicate it:
a = np.array([1, 2, 3])
b = np.tile(a, 5)
print b
[1 2 3 1 2 3 1 2 3 1 2 3 1 2 3]
c = np.tile(a, (5, 1))
print c
[[1 2 3]
 [1 2 3]
 [1 2 3]
 [1 2 3]
 [1 2 3]]

You can also create 2D arrays from 1D vectors:
x = np.linspace(-1, 1, 100)
y = x
X, Y = np.meshgrid(x, y) #this create 2D arrays containing x and y for each pixel
print X.shape
(100, 100)
Outer products are obtain using... np.outer:
a = np.array([1, 2, 3])
b = np.array([1, 10, 100])
np.outer(a, b)
array([[  1,  10, 100],
       [  2,  20, 200],
       [  3,  30, 300]])
np.outer(b, a)
array([[  1,   2,   3],
       [ 10,  20,  30],
       [100, 200, 300]])

Broadcasting

Python numpy is able to add dimensions to arrays to perform operations:
a = np.array([1, 2, 3, 4])
b = np.ones((6, 4))
a * b
array([[ 1.,  2.,  3.,  4.],
       [ 1.,  2.,  3.,  4.],
       [ 1.,  2.,  3.,  4.],
       [ 1.,  2.,  3.,  4.],
       [ 1.,  2.,  3.,  4.],
       [ 1.,  2.,  3.,  4.]])
Nice page on this:
http://www.scipy.org/EricsBroadcastingDoc

Indexing and Slicing

The access to elements of arrays is done using []:
a = np.arange(10)
a[5]
0
a[-1]  # last elements
9

In case of N-dims arrays, one can extract slides using :
a = np.array([1, 2, 3, 4, 5])
b = np.array([1, 10, 100, 1000])
c = np.outer(a, b)
c.shape
(5, 4)
print c 
[[   1   10  100 1000]
 [   2   20  200 2000]
 [   3   30  300 3000]
 [   4   40  400 4000]
 [   5   50  500 5000]]
print c[:,2]
[100 200 300 400 500]
print c[-1,:]
[   5   50  500 5000]
print c[-1, -1]
5000

There is a lot of methods in the array object, the best page to learn more is there: http://scipy-lectures.github.com/intro/numpy/numpy.html#the-numpy-array-object

Filtering

Sometimes we need to extract elements from arrays following some criteria. There is basically two ways to do this: defining a set of indices where the condition is completed (a la WHERE in IDL), or defining a boolean mask.

Let's first define a 2D array made of 10 times 1000 random values:
a = np.random.random((10, 1000))

We want to extract the values where the 2nd and the 4th 1000-elements vectors are greater than 0.5.

Using the numpy.where function:
w1 = np.where((a[1,:] > 0.5) & (a[3,:] > 0.5))
as the result is a tuple of indices, and in this case the dimension of the result is 1, better extract on the fly the array of indices from the tuple:
w2 = np.where((a[1,:] > 0.5) & (a[3,:] > 0.5))[0]
len(w2)
248  # your result may differ form this value, but must be close to 250.

We can now reduce the initial array to the desired values:
b = a[:, w2]
b.shape
(10, 248)

One can also build an array of boolean:
mask = (a[1,:] > 0.5) & (a[3,:] > 0.5)
which is actually similar to
mask2 =  np.where((a[1,:] > 0.5) & (a[3,:] > 0.5), True, False)

It can be used the same way as before:
b = a[:, mask]
b.shape
(10, 248)

These latter cases result in a 1000-element array, filled with True and False. To know the number of True values, just sum the array:
mask.sum()
248  # your result may differ form this value, but must be close to 250.

One big advantage of the mask technique is that you can combine different masks:
mask1 = a[1,:] > 0.5
mask2 = a[3,:] > 0.5
mask_total = mask1 & mask2
b = a[:, mask_total]

Sorting

a = np.array([1,2,45,2,1,46,7,-1])
b1 = np.sort(a)
ib = np.argsort(a)
b2 = a[ib]
print b1
[-1  1  1  2  2  7 45 46]
print b2 
[-1  1  1  2  2  7 45 46]

More on numpy arrays

http://scipy.org/Numpy_Example_List_With_Doc
http://scipy-lectures.github.com/intro/numpy/numpy.html#the-numpy-array-object

Sunday, April 8, 2012

Plotting


Don't miss the update: "2014 Python lecture on Matplotlib" on this blog: http://python-astro.blogspot.mx/2014/09/2014-python-lecture.html

Plotting is certainly one of the most common task one would do in Python (at least for astronomers).
There is various libraries to draw plots in Python, but the mostly used and powerful is perhaps matplotlib.
The plotting part is pyplot, so before any use, one must import it, with the (widely used) alias plt:
import matplotlib.pyplot as plt

The best to see how it works is to dive into the website where there is a lot of examples of plots with the code used to generate them: http://matplotlib.sourceforge.net/, especially the gallery: http://matplotlib.sourceforge.net/gallery.html
In the following we will present a few examples to help you to start.

X-Y plot

Let's say you have 2 vectors x and y and wanna plot y vs. x. Here are some examples of how to do this.
First generate x and y (don't forget to import numpy as np):
x = np.linspace(0., 4 * np.pi, 100)
y = np.sin(x)
plt.plot(x, y)
plt.show()

This latest command is not always necessary, if you run ipython with the --pylab option.
As you can see, python automatically define the axis ranges and draw the plot using a default blue color.
If you want to overplot another function, just call plt.plot again:
plt.plot(x, y**2)
The figure can be cleaned before another plot using
plt.clf()
You may need to draw various figure in different windows at the same time. Every call to
plt.figure()
will start a new figure where the next plot(s) will be drawn. The figures (each one in a separate window) are identified with numbers, the first one being 1.
When calling plt.figure(N), a new figure is created if figure N doesn't exist, and focus will be on figure N otherwise.

The color of the line is defined using c or color keyword:
plt.plot(x, y**3, color = 'red')
plt.plot(x, y**4, c='b')
Abbreviation Color
b blue
g green
r red
c cyan
m magenta
y yellow
k black
w white


Symbols and line style can also be easily defined:
plt.plot(x, y**3, color='r', marker='o', linestyle=':')
Symbol Description
- solid line
-- dashed line
-. dash-dot line
: dotted line
. points
, pixels
o circle symbols
^ triangle up symbols
v triangle down symbols
< triangle left symbols
> triangle right symbols
s square symbols
+ plus symbols
x cross symbols
D diamond symbols
d thin diamond symbols
1 tripod down symbols
2 tripod up symbols
3 tripod left symbols
4 tripod right symbols
h hexagon symbols
H rotated hexagon symbols
p pentagon symbols
| vertical line symbols
_ horizontal line symbols
steps use gnuplot style ‘steps’ # kwarg only

Other useful line properties:

Property Value
alpha alpha transparency on 0-1 scale
antialiased True or False - use antialised rendering
color matplotlib color arg
data_clipping whether to use numeric to clip data
label string optionally used for legend
linestyle one of - : -. -
linewidth float, the line width in points
marker one of + , o . s v x > <, etc
markeredgewidth line width around the marker symbol
markeredgecolor edge color if a marker is used
markerfacecolor face color if a marker is used
markersize size of the marker in points

If you need more control on the markers, better use scatter. For example, if one need ot change the size of the symbol according to a function:
plt.scatter(x, y**3, c=abs(y), marker='s', s=10+abs(y)*100, edgecolors='none')

The labels are set after the plot is done:
plt.xlabel('X')
plt.ylabel('Some trig function')
LaTex fans are welcome:
plt.title(r'Example of LaTex $\alpha_{\beta}$') #Notice the r before the string


One can change the axis ranges afterward:
plt.xlim((0, 10))
plt.ylim((-2, 2))

Multiple plots

To draw multiple plots of the same figure:
plt.subplot(N_y, N_x, N)

for i in np.arange(9):
    plt.subplot(3,3,i+1)
    plt.plot(x, y**i)
    plt.ylim((-2, 2))

log plots

using plt.semilogx, plt.semilogy and plt.loglog
plt.loglog(x)
plt.grid(True, which='minor')

Contours

Contour plots are done with plt.contour and plt.contourf.
x = np.linspace(-1, 1, 100)
y = x
X, Y = np.meshgrid(x, y) #this create 2D arrays containing x and y for each pixel
dist =  (X**2 + Y**2)**0.5
plt.contourf(X, Y, dist)
plt.colorbar()
CS = plt.contour(X, Y, dist, colors ='black', linewidths = 4)
plt.clabel(CS) #to print label on each contour



Saving the plot

The result of the plot can be save in PDF, EPS, JPG, BMP format, using:
plt.savefig('fig1.pdf')

More on plotting:

http://scipy-lectures.github.com/intro/matplotlib/matplotlib.html