my_int = 5
print my_int, type(my_int)

my_float = 5.0
print my_float, type(my_float)

my_boolean = False
print my_boolean, type(my_boolean)

my_string = 'hello'
print my_string, type(my_string)

my_list = [my_int, my_float, my_boolean, my_string]
print type(my_list)

for element in my_list:
    print type(element)

print my_list[1]

my_list[1] = 3.0
my_sublist = my_list[1:3]

print my_sublist
print type(my_sublist)

# Function with a required and an optional argument
def regress(x, c=0, b=1):
    return (x*b)+c

print regress(5)         # Only required argument
print regress(5, 10, 3)  # Use argument order
print regress(5, b=3)    # Specify the name to skip an optional argument

# Function without return argument
def divisible(a,b):
    if a%b:
        print str(a) + " is not divisible by " + str(b)
    else:
        print str(a) + " is divisible by " + str(b)

divisible(9,3)
res = divisible(9,2)
print res

# Function with multiple return arguments
def add_diff(a,b):
    return a+b, a-b

# Assigned as a tuple
res = add_diff(5,3)
print res

# Directly unpacked to two variables
a,d = add_diff(5,3)
print a
print d

my_list = [1, False, 'boo']
my_list.append('extra element')
my_list.remove(False)
print my_list

return_arg = my_list.append('another one')
print return_arg
print my_list

my_string = 'kumbaya, milord'
return_arg = my_string.replace('lord', 'lard')
print return_arg
print my_string

subj_length = [180.0,165.0,190.0,172.0,156.0]
subj_weight = [75.0,60.0,83.0,85.0,62.0]
subj_bmi = []

# EXERCISE 1: Try to compute the BMI of each subject, as well as the average BMI across subjects
# BMI = weight/(length/100)**2

subj_bmi = subj_weight/(subj_length/100)**2 
mean_bmi = mean(subj_bmi)

import numpy as np

# Create a numpy array from a list
subj_length = np.array([180.0,165.0,190.0,172.0,156.0])
subj_weight = np.array([75.0,60.0,83.0,85.0,62.0])

print type(subj_length), type(subj_weight)

# EXERCISE 2: Try to complete the program now!
# Hint: np.mean() computes the mean of a numpy array
# Note that unlike MATLAB, Python does not need the '.' before elementwise operators

# Multi-dimensional lists are just nested lists
# This is clumsy to work with
my_nested_list = [[1,2,3],[4,5,6]]

print my_nested_list

print len(my_nested_list)
print my_nested_list[0]
print len(my_nested_list[0])

# Numpy arrays handle multidimensionality better
arr = np.array(my_nested_list)

print arr         # nicer printing
print arr.shape   # direct access to all dimension sizes
print arr.size    # direct access to the total number of elements
print arr.ndim    # direct access to the number of dimensions

arr3d = np.array([ [[1,2,3],[4,5,6]] , [[7,8,9],[10,11,12]] ])

print arr3d

print arr3d.shape
print arr3d.size
print arr3d.ndim

# The type of a numpy array is always... numpy.ndarray
arr = np.array([[1,2,3],[4,5,6]])
print type(arr)

# So, let's do a computation
print arr/2

# Apparently we're doing our computations on integer elements!
# How do we find out?
print arr.dtype

# And how do we fix this?
arr = arr.astype('float')    # Note: this is not an in-place function!
print arr.dtype
print arr/2

# Alternatively, we could have defined our dtype better from the start
arr = np.array([[1,2,3],[4,5,6]], dtype='float')
print arr.dtype

arr = np.array([[1.,2.,3.],[4.,5.,6.]])
print arr.dtype

arr = np.array([[1,2,3],[4,5,6],[7,8,9]], dtype='float')

# Indexing and slicing
print arr[0,0]      # or: arr[0][0]
print arr[:-1,0]

# Elementwise computations on slices
# Remember, the LAST dimension is the INNER dimension
print arr[:,0] * arr[:,1]
print arr[0,:] * arr[1,:]

# Note that you could never slice across rows like this in a nested list!

# This doesn't work
# print arr[1:,0] * arr[:,1]

# And here's why:
print arr[1:,0].shape, arr[:,1].shape

# This however does work. You can always use scalars as the other operand.
print arr[:,0] * arr[2,2]

# Or, similarly:
print arr[:,0] * 9.

# EXERCISE 3: Create a 2x3 array containing the column-wise and the row-wise means of the original matrix
# Do not use a for-loop, and also do not use the np.mean() function for now.
arr = np.array([[1,2,3],[4,5,6],[7,8,9]], dtype='float')

# 1-D array, filled with zeros
arr = np.zeros(3)
print arr

# Multidimensional array of a given shape, filled with ones
# This automatically allows you to fill arrays with /any/ value
arr = np.ones((3,2))*5
print arr

# Sequence from 1 to AND NOT including 16, in steps of 3
# Note that using a float input makes the dtype a float as well
# This is equivalent to np.array(range(1.,16.,3))
arr = np.arange(1.,16.,3)
print arr

# Sequence from 1 to AND including 16, in 3 steps
# This always returns an array with dtype float
arr = np.linspace(1,16,3)
print arr

# Array of random numbers between 0 and 1, of a given shape
# Note that the inputs here are separate integers, not a tuple
arr = np.random.rand(5,2)
print arr

# Array of random integers from 0 to AND NOT including 10, of a given shape
# Here the shape is defined as a tuple again
arr = np.random.randint(0,10,(5,2))
print arr

arr0 = np.array([[1,2],[3,4]])
print arr0

# 'repeat' replicates elements along a given axis
# Each element is replicated directly after itself
arr = np.repeat(arr0, 3, axis=-1)
print arr

# We may even specify the number of times each element should be repeated
# The length of the tuple should correspond to the dimension length
arr = np.repeat(arr0, (2,4), axis=0)
print arr

print arr0

# 'tile' replicates the array as a whole
# Use a tuple to specify the number of tilings along each dimensions
arr = np.tile(arr0, (2,4))
print arr

# 'meshgrid' is commonly used to create X and Y coordinate arrays from two vectors
# where each array contains the X or Y coordinates corresponding to a given pixel in an image
x = np.arange(10)
y = np.arange(5)
print x,y

arrx, arry = np.meshgrid(x,y)
print arrx
print arry

arr0 = np.array([[1,2],[3,4]])
arr1 = np.array([[5,6],[7,8]])

# 'concatenate' requires an axis to perform its operation on
# The original arrays should be put in a tuple
arr = np.concatenate((arr0,arr1), axis=0)
print arr # as new rows

arr = np.concatenate((arr0,arr1), axis=1)
print arr # as new columns

# Suppose we want to create a 3-D matrix from them,
# we have to create them as being three-dimensional
# (what happens if you don't?)
arr0 = np.array([[[1],[2]],[[3],[4]]])
arr1 = np.array([[[5],[6]],[[7],[8]]])
print arr0.shape, arr1.shape
arr = np.concatenate((arr0,arr1),axis=2)
print arr

# hstack, vstack, and dstack are short-hand functions
# which will automatically create these 'missing' dimensions
arr0 = np.array([[1,2],[3,4]])
arr1 = np.array([[5,6],[7,8]])

# vstack() concatenates rows
arr = np.vstack((arr0,arr1))
print arr

# hstack() concatenates columns
arr = np.hstack((arr0,arr1))
print arr

# dstack() concatenates 2D arrays into 3D arrays
arr = np.dstack((arr0,arr1))
print arr

# Their counterparts are the hsplit, vsplit, dsplit functions
# They take a second argument: how do you want to split
arr = np.random.rand(4,4)
print arr
print '--'

# Splitting int equal parts 
arr0,arr1 = hsplit(arr,2)
print arr0
print arr1
print '--'

# Or, specify exact split points
arr0,arr1,arr2 = hsplit(arr,(1,2))
print arr0
print arr1
print arr2

arr0 = np.arange(10)
print arr0
print '--'

# 'reshape' does exactly what you would expect
# Make sure though that the total number of elements remains the same
arr = np.reshape(arr0,(5,2))
print arr

# You can also leave one dimension blank by using -1 as a value
# Numpy will then compute for you how long this dimension should be
arr = np.reshape(arr0,(-1,5))
print arr
print '--'

# 'transpose' allows you to switch around dimensions
# A tuple specifies the new order of dimensions
arr = np.transpose(arr,(1,0))
print arr

# For simply transposing rows and columns, there is the short-hand form .T
arr = arr.T
print arr
print '--'

# 'flatten' creates a 1D array out of everything
arr = arr.flatten()
print arr

# EXERCISE 4: Create your own meshgrid3d function
# Like np.meshgrid(), it should take two vectors and replicate them; one into columns, the other into rows
# Unlike np.meshgrid(), it should return them as a single 3D array rather than 2D arrays
# ...do not use the np.meshgrid() function

def meshgrid3d(xvec, yvec):
    # fill in!

xvec = np.arange(10)
yvec = np.arange(5)
xy = meshgrid3d(xvec, yvec)
print xy
print xy[:,:,0] # = first output of np.meshgrid()
print xy[:,:,1] # = second output of np.meshgrid()

arr = np.random.rand(5)
print arr

# Sorting and shuffling
res = arr.sort() 
print arr # in-place!!!
print res

res = np.random.shuffle(arr) 
print arr # in-place!!!
print res

# Min, max, mean, standard deviation
arr = np.random.rand(5)
print arr

mn = np.min(arr)
mx = np.max(arr)
print mn, mx

mu = np.mean(arr)
sigma = np.std(arr)
print mu, sigma

# Some functions allow you to specify an axis to work along, in case of multidimensional arrays
arr2d = np.random.rand(3,5)
print arr2d

print np.mean(arr2d, axis=0)
print np.mean(arr2d, axis=1)

# Trigonometric functions
# Note: Numpy works with radians units, not degrees
arr = np.random.rand(5)
print arr

sn = np.sin(arr*2*np.pi)
cs = np.cos(arr*2*np.pi)
print sn
print cs

# Exponents and logarithms
arr = np.random.rand(5)
print arr

xp = np.exp(arr)
print xp
print np.log(xp)

# Rounding
arr = np.random.rand(5)
print arr

print arr*5
print np.round(arr*5)
print np.floor(arr*5)
print np.ceil(arr*5)

# EXERCISE 5: Make a better version of Exercise 3 with what you've just learned
arr = np.array([[1,2,3],[4,5,6],[7,8,9]], dtype='float')

# What we had:
print np.array([(arr[:,0]+arr[:,1]+arr[:,2])/3,(arr[0,:]+arr[1,:]+arr[2,:])/3])

# Now the new version:

# EXERCISE 6: Create a Gabor patch of 100 by 100 pixels
import numpy as np
import matplotlib.pyplot as plt

# Step 1: Define the 1D coordinate values
# Tip: use 100 equally spaced values between -np.pi and np.pi

# Step 2: Create the 2D x and y coordinate arrays
# Tip: use np.meshgrid()

# Step 3: Create the grating
# Tip: Use a frequency of 10

# Step 4: Create the Gaussian
# Tip: use np.exp() to compute a power of e

# Step 5: Create the Gabor

# Visualize your result
# (we will discuss how this works later)
plt.figure(figsize=(15,5))
plt.subplot(131)
plt.imshow(grating, cmap='gray')
plt.subplot(132)
plt.imshow(gaussian, cmap='gray')
plt.subplot(133)
plt.imshow(gabor, cmap='gray')
plt.show()

# Check whether each element of a 2x2 array is greater than 0.5
arr = np.random.rand(2,2)
print arr
res = arr>0.5
print res
print '--'
 
# Analogously, check it against each element of a second 2x2 array
arr2 = np.random.rand(2,2)
print arr2
res = arr>arr2
print res

# We can use these boolean arrays as indices into other arrays!

# Add 0.5 to any element smaller than 0.5
arr = np.random.rand(2,2)
print arr
res = arr<0.5
print res
arr[res] = arr[res]+0.5
print arr

# Or, shorter:
arr[arr<0.5] = arr[arr<0.5] + 0.5

# Or, even shorter:
arr[arr<0.5] += 0.5

arr = np.array([[1,2,3],[4,5,6]])

# The short-hand forms for elementwise boolean operators are: & | ~ ^
# Use parentheses around such expressions
res = (arr<4) & (arr>1)
print res
print '--'

res = (arr<2) | (arr==5)
print res
print '--'

res = (arr>3) & ~(arr==6)
print res
print '--'

res = (arr>3) ^ (arr<5)
print res

# To convert boolean indices to normal integer indices, use the 'nonzero' function
print res
print np.nonzero(res)
print '--'

# Separate row and column indices
print np.nonzero(res)[0]
print np.nonzero(res)[1]
print '--'

# Or stack and transpose them to get index pairs
pairs = np.vstack(np.nonzero(res)).T
print pairs

import numpy as np

# We will keep track of the sum of first occurence positions,
# as well as the number of positions entered into this sum.
# This way we can compute the mean.
sum111 = 0.
n111 = 0.
sum123 = 0.
n123 = 0.

for sim in range(5000):
    
    # Keep track of how far along we are in finding a given pattern
    d111 = 0
    d123 = 0
    
    for throw in range(2000):
        
        # Throw a die
        die = np.random.randint(1,7)
        
        # 111 case
        if d111==3:
            pass
        elif die==1 and d111==0:
            d111 = 1
        elif die==1 and d111==1:
            d111 = 2
        elif die==1 and d111==2:
            d111 = 3
            sum111 = sum111 + throw
            n111 = n111 + 1            
        else:
            d111 = 0
      
        # 123 case
        if d123==3:
            pass
        elif die==1:
            d123 = 1
        elif die==2 and d123==1:
            d123 = 2
        elif die==3 and d123==2:
            d123 = 3
            sum123 = sum123 + throw
            n123 = n123 + 1
        else:
            d123 = 0
        
        # Don't continue if both have been found
        if d111==3 and d123==3:
            break

# Compute the averages
avg111 = sum111/n111
avg123 = sum123/n123
print avg111, avg123   

# ...can you spot the crucial difference between both patterns?

# EXERCISE 7: Vectorize the above program

# You get these lines for free...
import numpy as np
throws = np.random.randint(1,7,(5000,2000))
one = (throws==1)
two = (throws==2)
three = (throws==3)

# Find out where all the 111 and 123 sequences occur
find111 = 
find123 = 

# Then at what index they /first/ occur in each sequence
first111 = 
first123 = 

# Compute the average first occurence location for both situations
avg111 = 
avg123 = 

# Print the result
print avg111, avg123

from PIL import Image

# Opening an image is simple enough:
# Construct an Image object with the filename as an argument
im = Image.open('python.jpg')

# It is now represented as an object of the 'JpegImageFile' type
print im

# There are some useful member variables we can inspect here
print im.format  # format in which the file was saved
print im.size    # pixel dimensions
print im.mode    # luminance/color model used

# We can even display it
# NOTE this is not perfect; meant for debugging
im.show()

# Alternative quick-show method
from Tkinter import Tk, Button
from PIL import ImageTk

def alt_show(im):
    win = Tk()
    tkimg = ImageTk.PhotoImage(im)
    Button(image=tkimg).pack()
    win.mainloop()

alt_show(im)

# We can convert PIL images to an ndarray!
arr = np.array(im)
print arr.dtype   # uint8 = unsigned 8-bit integer (values 0-255 only)
print arr.shape   # Why do we have three layers?

# Let's make it a float-type for doing computations
arr = arr.astype('float')
print arr.dtype

# This opens up unlimited possibilities for image processing!
# For instance, let's make this a grayscale image, and add white noise
max_noise = 50
arr = np.mean(arr,-1)
noise = (np.random.rand(arr.shape[0],arr.shape[1])-0.5)*2
arr = arr + noise*max_noise

# Make sure we don't exceed the 0-255 limits of a uint8
arr[arr<0] = 0
arr[arr>255] = 255

# When going back to PIL, it's a good idea to explicitly 
# specify the right dtype and the mode.
# Because automatic conversions might mess things up
arr = arr.astype('uint8')
imn = Image.fromarray(arr, mode='L')
print imn.format
print imn.size
print imn.mode  # L = greyscale
imn.show() # or use alt_show() from above if show() doesn't work well for you

# Note that /any/ 2D or 2Dx3 numpy array filled with values between 0 and 255
# can be converted to an image object in this way

im = Image.open('python.jpg')

# Make the image smaller
ims = im.resize((800,600))
ims.show()

# Or you could even make it larger
# The resample argument allows you to specify the method used
iml = im.resize((1280,1024), resample=Image.BILINEAR)
iml.show()

# Rotation is similar (unit=degrees)
imr = im.rotate(10, resample=Image.BILINEAR, expand=False)
imr.show()

# If we want to lose the black corners, we can crop (unit=pixels)
imr = imr.crop((100,100,924,668))
imr.show()

# 'convert' allows conversion between different color models
# The most important here is between 'L' (luminance) and 'RGB' (color)
imbw = im.convert('L')
imbw.show()
print imbw.mode

imrgb = imbw.convert('RGB')
imrgb.show()
print imrgb.mode

# Note that the grayscale conversion of PIL is more sophisticated
# than simply averaging the three layers in Numpy (it is a weighted average)

# Also note that the color information is effectively lost after converting to L

from PIL import Image, ImageFilter

im = Image.open('python.jpg')
imbw = im.convert('L')

# Contour detection filter
imf = imbw.filter(ImageFilter.CONTOUR)
imf.show()

# Blurring filter
imf = imbw.filter(ImageFilter.GaussianBlur(radius=3))
imf.show()

from PIL import Image, ImageDraw

im = Image.open('python.jpg')

# You need to attach a drawing object to the image first
imd = ImageDraw.Draw(im)

# Then you work on this object
imd.rectangle([10,10,100,100], fill=(255,0,0))
imd.line([(200,200),(200,600)], width=10, fill=(0,0,255))
imd.text([500,500], 'Python', fill=(0,255,0))

# The results are automatically applied to the Image object
im.show()

# PIL will figure out the file type by the extension
im.save('python.bmp')

# There are also further options, like compression quality (0-100)
im.save('python_bad.jpg', quality=5)

# EXERCISE 8: Visualize the difference between the PIL conversion to grayscale, and a simple averaging of RGB
# Display pixels where the average is LESS luminant in red, and where it is MORE luminant in shades green
# The luminance of these colors should correspond to the size of the difference
#
# Extra 1: Maximize the overall contrast in your image
#
# Extra 2: Save as three PNG files, of different sizes (large, medium, small)

import numpy as np
from PIL import Image
import matplotlib.pyplot as plt

# As data for our plots, we will use the pixel values of the image
# Open image, convert to an array
im = Image.open('python.jpg')
im = im.resize((400,300))
arr = np.array(im, dtype='float')

# Split the RGB layers and flatten them
R,G,B = np.dsplit(arr,3)
R = R.flatten()
G = G.flatten()
B = B.flatten()

# QUICKPLOT 1: Correlation of luminances in the image

# This works if you want to be very quick:
# (xb means blue crosses, .g are green dots) 
plt.plot(R, B, 'xb')
plt.plot(R, G, '.g')

# However we will take a slightly more disciplined approach here
# Note that Matplotlib wants colors expressed as 0-1 values instead of 0-255

# Create a square figure
plt.figure(figsize=(5,5))

# Plot both scatter clouds
# marker: self-explanatory
# linestyle: 'None' because we want no line
# color: RGB triplet with values 0-1
plt.plot(R, B, marker='x', linestyle='None', color=(0,0,0.6))
plt.plot(R, G, marker='.', linestyle='None', color=(0,0.35,0))

# Make the axis scales equal, and name them
plt.axis([0,255,0,255])
plt.xlabel('Red value')
plt.ylabel('Green/Blue value')

# Show the result
plt.show()

# QUICKPLOT 2: Histogram of 'red' values in the image
plt.hist(R)

# ...and now a nicer version

# Make a non-square figure
plt.figure(figsize=(7,5))

# Make a histogram with 25 red bins
# Here we simply use the abbreviation 'r' for red
plt.hist(R, bins=25, color='r')

# Set the X axis limits and label
plt.xlim([0,255])
plt.xlabel('Red value', size=16)

# Remove the Y ticks and labels by setting them to an empty list
plt.yticks([])

# Remove the top ticks by specifying the 'top' argument
plt.tick_params(top=False)

# Add two vertical lines for the mean and the median
plt.axvline(np.mean(R), color='g', linewidth=3, label='mean')
plt.axvline(np.median(R), color='b', linewidth=1, linestyle=':', label='median')

# Generate a legend based on the label= arguments
plt.legend(loc=2)

# Show the plot
plt.show()

# QUICKPLOT 3: Bar chart of mean+std of RGB values
plt.bar([0,1,2],[np.mean(R), np.mean(G), np.mean(B)], yerr=[np.std(R), np.std(G), np.std(B)])

# ...and now a nicer version

# Make a non-square-figure
plt.figure(figsize=(7,5))

# Plot the bars with various options
# x location where bars start, y height of bars
# yerr: data for error bars
# width: width of the bars
# color: surface color of bars
# ecolor: color of error bars ('k' means black)
plt.bar([0,1,2], [np.mean(R), np.mean(G), np.mean(B)], 
         yerr=[np.std(R), np.std(G), np.std(B)],
         width=0.75,
         color=['r','g','b'], 
         ecolor='k')

# Set the X-axis limits and tick labels
plt.xlim((-0.25,3.))
plt.xticks(np.array([0,1,2])+0.75/2, ['Red','Green','Blue'], size=16)

# Remove all X-axis ticks by setting their length to 0
plt.tick_params(length=0)

# Set a figure title
plt.title('RGB Color Channels', size=16)

# Show the figure
plt.show()

# So, copy-paste this line into the box above, before the plt.show() command
plt.savefig('bar.png') 

# There are some further formatting options possible, e.g.
plt.savefig('bar.svg', dpi=300, bbox_inches=('tight'), pad_inches=(1,1), facecolor=(0.8,0.8,0.8))

# A simple grayscale luminance map
# cmap: colormap used to display the values
plt.figure(figsize=(5,5))
plt.imshow(np.mean(arr,2), cmap='gray')
plt.show()

# Importantly and contrary to PIL, imshow luminances are by default relative
# That is, the values are always rescaled to 0-255 first (maximum contrast)
# Moreover, colormaps other than grayscale can be used
plt.figure(figsize=(5,5))
plt.imshow(np.mean(arr,2)+100, cmap='jet') # or hot, hsv, cool,...
plt.show() # as you can see, adding 100 didn't make a difference here

# 'Figure' objects are returned by the plt.figure() command
fig = plt.figure(figsize=(7,5))
print type(fig)

# Axes objects are the /actual/ plots within the figure
# Create them using the add_axes() method of the figure object
# The input coordinates are relative (left, bottom, width, height)
ax0 = fig.add_axes([0.1,0.1,0.4,0.7], xlabel='The X Axis')
ax1 = fig.add_axes([0.2,0.2,0.5,0.2], axisbg='gray')
ax2 = fig.add_axes([0.4,0.5,0.4,0.4], projection='polar')
print type(ax0), type(ax1), type(ax2)

# This allows you to execute functions like savefig() directly on the figure object
# This resolves Matplotlib's confusion of what the current figure is, when using plt.savefig()
fig.savefig('fig.png')

# It also allows you to add text to the figure as a whole, across the different axes objects
fig.text(0.5, 0.5, 'splatter', color='r')

# The overall figure title can be set separate from the individual plot titles
fig.suptitle('What a mess', size=18)

# show() is actually a figure method as well
# It just gets 'forwarded' to what is thought to be the current figure if you use plt.show()
fig.show()

# Create a new figure
fig = plt.figure(figsize=(15,10))

# As we saw, many of the axes properties can already be set at their creation
ax0 = fig.add_axes([0.,0.,0.25,0.25], xticks=(0.1,0.5,0.9), xticklabels=('one','thro','twee'))
ax1 = fig.add_axes([0.3,0.,0.25,0.25], xscale='log', ylim=(0,0.5))
ax2 = fig.add_axes([0.6,0.,0.25,0.25])

# Once you have the axes object though, there are further methods available
# This includes many of the top-level pyplot functions
# If you use for instance plt.plot(), Matplotlib is actually 'forwarding' this
# to an Axes.plot() call on the current Axes object
R.sort()
G.sort()
B.sort()
ax2.plot(R, color='r', linestyle='-', marker='None')  # plot directly to an Axes object of choice
plt.plot(G, color='g', linestyle='-', marker='None')  # plt.plot() just plots to the last created Axes object
ax2.plot(B, color='b', linestyle='-', marker='None')

# Other top-level pyplot functions are simply renamed to 'set_' functions here
ax1.set_xticks([])
plt.yticks([])

# Show the figure
fig.show()

# Create a new figure
fig = plt.figure(figsize=(15,5))

# Specify the LAYOUT of the subplots (rows,columns)
# as well as the CURRENT Axes you want to work on
ax0 = fig.add_subplot(231)

# Equivalent top-level call on the current figure
# It is also possible to create several subplots at once using plt.subplots()
ax1 = plt.subplot(232)      

# Optional arguments are similar to those of add_axes()
ax2 = fig.add_subplot(233, title='three') 

# We can use these Axes object as before
ax3 = fig.add_subplot(234)
ax3.plot(R, 'r-') 
ax3.set_xticks([])
ax3.set_yticks([])

# We skipped the fifth subplot, and create only the 6th
ax5 = fig.add_subplot(236, projection='polar')

# We can adjust the spacings afterwards
fig.subplots_adjust(hspace=0.4)

# And even make room in the figure for a plot that doesn't fit the grid
fig.subplots_adjust(right=0.5)
ax6 = fig.add_axes([0.55,0.1,0.3,0.8])

# Show the figure
fig.show()

# EXERCISE 9: Plot y=sin(x) and y=sin(x^2) in two separate subplots, one above the other
# Let x range from 0 to 2*pi

# This uses the result of the exercise above
# You have to copy-paste it into the same code-box, before the fig.show()

# Add horizontal lines
ax0.axhline(0, color='g')
ax0.axhline(0.5, color='gray', linestyle=':')
ax0.axhline(-0.5, color='gray', linestyle=':')
ax1.axhline(0, color='g')
ax1.axhline(0.5, color='gray', linestyle=':')
ax1.axhline(-0.5, color='gray', linestyle=':')

# Add text to the plots
ax0.text(0.1,-0.9,'$y = sin(x)$', size=16)   # math mode for proper formula formatting!
ax1.text(0.1,-0.9,'$y = sin(x^2)$', size=16)

# Annotate certain points with a value
for x_an in np.linspace(0,2*np.pi,9):
    ax0.annotate(str(round(sin(x_an),2)),(x_an,sin(x_an)))
    
# Add an arrow (x,y,xlength,ylength)
ax0.arrow(np.pi-0.5,-0.5,0.5,0.5, head_width=0.1, length_includes_head=True)

# This uses the result of the exercise above
# You have to copy-paste it into the same code-box, before the fig.show()

# For instance, fetch the X axis
# XAxis objects have their own methods
xax = ax1.get_xaxis()
print type(xax)

# These methods allow you to fetch the even smaller building blocks
# For instance, tick-lines are Line2D objects attached to the XAxis
xaxt = xax.get_majorticklines()
print len(xaxt)

# Of which you can fetch AND change the properties
# Here we change just one tickline into a cross
print xaxt[6].get_color()
xaxt[6].set_color('g')
xaxt[6].set_marker('x')
xaxt[6].set_markersize(10)

# This uses the result of the exercise above
# You have to copy-paste it into the same code-box, before the fig.show()

# Another example: fetch the lines in the plot
# Change the color, change the marker, and mark only every 100 points for one specific line
ln = ax0.get_lines()
print ln
ln[0].set_color('g')
ln[0].set_marker('o')
ln[0].set_markerfacecolor('b')
ln[0].set_markevery(100)

# Finally, let's create a graphic element from scratch, that is not available as a top-level pyplot function
# And then attach it to existing Axes
# NOTE: we need to import something before we can create the ellipse like this. What should we import?
ell = matplotlib.patches.Ellipse((np.pi,0), 1., 1., color='r')
ax0.add_artist(ell)
ell.set_hatch('//')
ell.set_edgecolor('black')
ell.set_facecolor((0.9,0.9,0.9))

# EXERCISE 10: Add regression lines
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
import matplotlib.lines as lines

# Open image, convert to an array
im = Image.open('python.jpg')
im = im.resize((400,300))
arr = np.array(im, dtype='float')

# Split the RGB layers and flatten them
R,G,B = np.dsplit(arr,3)
R = R.flatten()
G = G.flatten()
B = B.flatten()

# Do the plotting
plt.figure(figsize=(5,5))
plt.plot(R, B, marker='x', linestyle='None', color=(0,0,0.6))
plt.plot(R, G, marker='.', linestyle='None', color=(0,0.35,0))

# Tweak the plot
plt.axis([0,255,0,255])
plt.xlabel('Red value')
plt.ylabel('Green/Blue value')

# Fill in your code...

# Show the result
plt.show()

import numpy as np
import scipy.stats as stats

# Generate random numbers between 0 and 1
data = np.random.rand(30)

# Do a t-test with a H0 for the mean of 0.4
t,p = stats.ttest_1samp(data,0.4)
print p

# Generate another sample of random numbers, with mean 0.4
data2 = np.random.rand(30)-0.1

# Do a t-test that these have the same mean
t,p = stats.ttest_ind(data, data2)
print p

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

# Simulate the size of the F statistic when comparing three conditions
# Given a constant n, and an increasing true effect size. 
true_effect = np.linspace(0,0.5,500)
n = 100
Fres = []

# Draw random normally distributed samples for each condition, and do a one-way ANOVA
for eff in true_effect:
    c1 = stats.norm.rvs(0,1,size=n)
    c2 = stats.norm.rvs(eff,1,size=n)
    c3 = stats.norm.rvs(2*eff,1,size=n)
    F,p = stats.f_oneway(c1,c2,c3)
    Fres.append(F)

# Create the plot
plt.figure()
plt.plot(true_effect,Fres,'r*-')
plt.xlabel('True Effect')
plt.ylabel('F')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

# Compute the pdf and cdf of normal distributions, with increasing sd's
# Then plot them in different colors
# (of course, many other distributions are also available)
x = np.linspace(-5,5,1000)
sds = np.linspace(0.25,2.5,10)
cols = np.linspace(0.15,0.85,10)

# Create the figure
fig = plt.figure(figsize=(10,5))
ax0 = fig.add_subplot(121)
ax1 = fig.add_subplot(122)

# Compute the densities, and plot them
for i,sd in enumerate(sds):
    y1 = stats.norm.pdf(x,0,sd)
    y2 = stats.norm.cdf(x,0,sd)
    ax0.plot(x,y1,color=cols[i]*np.array([1,0,0]))
    ax1.plot(x,y2,color=cols[i]*np.array([0,1,0]))

# Show the figure
plt.show()

import numpy as np
import scipy.fftpack as fft

# The original data: a step function
data = np.zeros(200, dtype='float')
data[25:100] = 1

# Decompose into sinusoidal components
# The result is a series of complex numbers as long as the data itself
res = fft.fft(data)

# FREQUENCY is implied by the ordering, but can be retrieved as well
# It increases from 0 to the Nyquist frequency (0.5), followed by its reversed negative counterpart
# Note: in case of real input data, the FFT results will be symmetrical
freq = fft.fftfreq(data.size)

# AMPLITUDE is given by np.abs() of the complex numbers
amp = np.abs(res)

# PHASE is given by np.angle() of the complex numbers
phase = np.angle(res)

# We can plot each component separately
plt.figure(figsize=(15,5))
plt.plot(data, 'k-', lw=3)
xs = np.linspace(0,data.size-1,data.size)*2*np.pi
for i in range(len(res)):   
    ys = np.cos(xs*freq[i]+phase[i]) * (amp[i]/data.size)
    plt.plot(ys.real, 'r:', lw=1)
plt.show()

# Can you then plot what the SUM of all these components looks like?

# ifft = inverse fft
reconstructed = fft.ifft(res)

plt.figure(figsize=(15,5))
plt.plot(data,'k-', lw=3)
plt.plot(reconstructed, 'r:', lw=3)
plt.show()

# Note that /some/ information has been lost, but very little
print 'Total deviation:', np.sum(np.abs(data-reconstructed))

# Set components with low frequencies equal to 0
resf = res.copy()
mask = np.abs(fft.fftfreq(data.size)) < 0.25
resf[mask] = 0

# ifft the modified components
filtered = fft.ifft(resf)

# And the result is high-pass filtered
plt.figure(figsize=(15,5))
plt.plot(data,'k-', lw=3)
plt.plot(filtered, 'r:', lw=3)
plt.show()

import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
import scipy.fftpack as fft

im = Image.open('python.jpg').convert('L')
arr = np.array(im, dtype='float')
res = fft.fft2(arr)

# Just set all amplitudes to one
new_asp = np.ones(res.shape, dtype='float')

# Then recombine the complex numbers
real_part = new_asp * np.cos(np.angle(res))
imag_part = new_asp * np.sin(np.angle(res))
eq = real_part+(imag_part*1j)

# And do the ifft
arr_eq = fft.ifft2(eq).real

# Show the result
# Clear the high frequencies dominate now
plt.figure()
plt.imshow(arr_eq, cmap='gray')
plt.show()