################################################################################
# Import statements. Python's library system is module based.  The
# ones we are interested in are numpy, scipy, and matplotlib's pylab
# plotting module.

import numpy as np
import scipy as sp
import matplotlib.pylab as plt

################################################################################
# As a calculator

2*3
2**8  # Exponentiation
3.6 / 1.2
13 % 4  # modulus

3/2  # Note -- be careful with integer types, behavior is as in C
3.0 / 2

np.sin(3.14159)
np.sqrt(16)
16 ** 0.5

################################################################################
# Basic data types.

####################
# Numbers
x = 3
x += 5
x 

####################
# Lists:

L1 = []           # empty list
L2 = [1,2,3]      # Contains numbers
L3 = [1,2,L2]     # Contains any other types
print L3

L2[0] = 5         # Note that all indexing is done based off of zeros
print L3          # Note that all objects are passed by reference.

L4 = range(10)    # Range constructs a list of the numbers 0,1,...,9
print L4

# Indexing
L4[:4]    
L4[1:4]
L4[::2]
L4[1:5:2]
L4[-1]
L4[:-1]
L4[5:1:-1]

# extending
L4 += [1,2]
print L4

L4.append(100)
print L4

# list based processing
L5 = [x**2 for x in L4]
L6 = [x**2 for x in L4 if x <= 5]

####################
# Tuples:  Like lists, but can't change the elements; these are immutable

t1 = (1,2,3)
t2 = (t1, "blah", 4)

####################
# Strings

s1 = "My name"
s2 = s1 + " is Inigo Montoya."
s2

s2.upper()  # Many string manipulation tools

s2[::-1]

answer = 42

# Formatting can be done using c-style syntax
print "The answer is %d." % answer  
print "The answer is",  answer, "." 

####################
# Dictionaries: Hash tables mapping from any immutable type to anything

d = {}
d["one"] = 1
d[2] = "two"
d[(1, "one")] = L3

print d

d2 = {2 : 4, 5 : 10}
d2

# Create dicts from length 2 tuples also
d3 = dict( (x, x**2) for x in range(10) )

################################################################################
# Program flow

x = 1 + 1

if x == 2:  # Also have <, >, <=, >=, != 
    print "1 + 1 = 2"


def f(x):
    if x <= 1:
        print "%d is too small." % x
    elif all([x % y != 0 for y in range(2, x)]):
        print "%d is prime." % x 
    else:
        print "%d is not prime." % x
    
for i in range(15):
    f(i)

def collatz(n):
    while n != 1:
        print n,   # Extra comma means a newline isn't printed
        n = n / 2 if n % 2 == 0 else 3*n + 1
    print '1 ...end!'

collatz(5)
collatz(20)
collatz(17)

################################################################################
# Numpy

# The basic class for numpy is the n-dimensional array.  You can
# create the array in many ways

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

X.shape # tuple containg the shape 
X.ndim  # number of dimensions
X.size  # number of elements

X.sum(axis=0)
X.sum(axis=1)
X.sum()

X.mean(axis=0)
X.mean(axis=1)

# Indexing
X[:2, :]
X[:, 1:3]

# Array manipulation
X2 = X.copy()
X2[1, 1:3] += 2
X2[:, ::2] = 0

# Combine things
X2 * (X2 + 1)

# Masks
X3 = X.copy()
X3 > 5          # true / false mask
X3[X3 > 5] = 0

X4 = X.copy()
X4[((X4 > 5) & (X4 < 10)) | (X4 == 3)] = 0

####################
# Other functions

np.cos(X * (3.14159 / 12) )

####################
# Other ways of making arrays
np.arange(10)   # Just like range, but creates arrays & can use fractional imputs
np.arange(5,10)
np.arange(5,10,0.5)

np.arange(12).reshape(3,4)

np.linspace(0, 5, 5)

np.linspace(0, 5, 6)

####################
# Random number generation 

np.random.seed(0)

x = np.random.uniform()
A = np.random.uniform( size = (2, 5) )

################################################################################
# Plotting -- matplotlib.pylab

# Interactive plotting
plt.figure()

plt.plot( [1,2,3], [1,3,2] )
plt.plot( [3,2,3], [1,3,2], 'xr')

plt.title("A random plot.")
plt.xlabel("Input.")
plt.ylabel("Value.")

plt.clf()

x = np.linspace(0,10, 500)

plt.plot(x, np.sin(x), label = "sin(x)")
plt.plot(x, np.cos(x), label = "cos(x)")

plt.legend()
plt.xlabel("x")
plt.ylabel("f(x)")

plt.show()

#######################################################################
# Also an object oriented version

fig = plt.figure()
ax = fig.add_subplot(2,1,1)

x1 = np.random.normal(0,1,50)
x2 = np.random.normal(1,1,50)
x3 = np.random.normal(2,1,50)

ax.set_title("The normal density -- 3 means.")
ax.boxplot([x1,x2,x3])

ax = fig.add_subplot(2,1,2)

x = np.random.normal(0,1,1000)
num_bins = 50
ax.hist(x,num_bins,color='green')

plt.show()