%pylab inline
# shouldn't be needed if you started the server with `ipython notebook --pylab inline` but otherwise...

import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['image.interpolation']= 'nearest'

import numpy as np

# np.array() - takes a list or tuple and makes an array.  core creation function
A = np.array([1, 2, 3, 4, 5])
A

B = np.array([[5.,6.,7.],[8.,9.,10.]])
print B

plt.figure()
plt.plot(A, 'o', ls='')

plt.figure()
plt.imshow(B)
plt.colorbar()

# np.arange() - works like range()
np.arange(5)

np.arange(2,30,4) # 2 to 30, every 4th value

# np.zeros(), np.ones(), np.empty() - take a tuple, each entry is the size in that dimension
A = np.zeros((4,10))
B = np.ones((15,2))
C = np.empty((5))
A, B, C

# np.zeros_like(), np.ones_like(), np.empty_like() - take their shape from what you pass in
# note that this also uses the datatype of the array, as well
B_2 = np.ones_like(A)
A, B_2

A = np.random.random((10,10))
plt.imshow(A)
plt.colorbar()
A[0,0]

# and assignment obviously works
A[0,0] = 1.0
plt.imshow(A)
plt.colorbar()

# using slices to get part of an array
# the notation is start:stop:interval, once for each dimension
# if you leave off a slice, numpy assumes you mean all of the remaining dimensions

plt.plot(A[:,0])
figure()
plt.imshow(A[:5,:])

for item in np.array([1,3,5,7,9]):
    print item

A = np.random.random((2,10))
for i in A:
    print i

# we can temporarily switch the axis to iterate over with np.rollaxis(), see below
# alternatively you can get the shape of the array and use those as loop variables.

A = np.arange(25).reshape(5,5)

print A.shape  # tuple of dimension sizes, I use this one constantly
print A.ndim   # equal to len(asdf.shape)
print A.dtype  # default on my platform is 64bit int
print A.size   # total number of elements in the array

# .copy() - return a copy of an array, not a view.  also is a method on an ndarray object
A = np.array([1,2,3,4])
B = A.copy()  # or - np.copy(asdf)
B[0] = 99

A, B

# .reshape() - takes an array and a tuple of the new size to change that array to
# and returns a view on the *same array*

A = np.array([1,2,3,4])
print A
print
B = A.reshape((2,2))
print B
print 
B[0,0] = 9
print A

# .astype() - casting method, can use built-in data types, a host of datatypes includined in NumPy, or strings
A = np.array([0,1,2,3,4])
print A.astype(float)
print A.astype(bool) # 0 is False, >0 is True

B = np.array([0.1,0.2,0.3,0.4,0.5])
print B.astype(int) # note the floor 

C = np.array([-4,-3,-2,-1,0,1,2,3,4]).astype('uint16')
print C # note wraparound

# finally, note that astype() returns a casted COPY of the array

# np.roll() - rotate the data along a given axis, returns a copy
A = np.array([1,2,3,4])
B = np.roll(A, 1)

print A, B
print '---------'

A = A.reshape((2,2))
C = np.roll(A, 1, axis=0)
D = np.roll(A, 1, axis=1)
print A
print
print C
print
print D

# .T - transpose method on an ndarray, returns a view on the original array
# unlike the others listed here, this is an attribute, not a method
A = np.arange(25).reshape(5,5)
print A
print
print A.T
print
print A.T[4,0]
A.T[4,0] = 99
print
print A[4,0]

# np.rollaxis() - reorder axes in an array, returning a view
# takes an array, the axis to reorder, and the position to roll before.
# the relative order of the other axes are not changed

A = np.zeros((2,3,4))
B = np.rollaxis(A, 2, 0)
print A.shape, B.shape
print '---------'
B[0,0,0] = 1
print A

# note that we can use rollaxis with iteration to get a slices of an array in a given dimension
# simply roll that axes to the start

A = np.zeros((2,3,4))
for array_slice in np.rollaxis(A,2,0):
    print array_slice.shape

# np.unique() - like set, but for ndarrays:
A = np.array([1,1,3,4,5,0])
print np.unique(A)

# .flatten() - return a copy of an multidimensional array flattened into one dimension
# compart to .flat - an attribute which is a 1d iterator over the whole array
A = np.arange(25).reshape(5,5)
B = A.flatten()

A, B

# .tolist() - convert an ndarray to a list - works for multidimensonial arrays
A = np.arange(25).reshape(5,5)
print A.tolist()

A = np.ones((5,5))
print A
print
print A+A
print
print A*A # 1*1 = 1!
print
print np.dot(A, A) # multiple matrices with np.dot

print A
print
print A*5

# np.sum(), .sum() - calculate the sum of the elements over a given axes, returning a new array
A = np.arange(25).reshape(5,5)
print A
print
print A.sum()
print
print A.sum(axis=0)
print
print A.sum(axis=1)

# np.mean(), .mean() - calcuated the average of elements over a given axis, returning a new array
A = np.arange(25).reshape(5,5)
print A
print
print A.mean()
print
print A.mean(axis=0)
print
print A.mean(axis=1)

# np.argmin(), .argmin(), np.argmax(), .argmax(), - return the index of the mininum or maximum value in the array
# note that these take the axis keyword- otherwise they return the max of the flattened array
# see also .argwhere()

A = np.arange(25)
B = A.reshape(5,5)
print A
print
print B

print np.argmin(A)
print
print np.argmax(B, axis=1)

# np.ceil, np.floor, np.round - rounding operations, element-wise
A = np.random.random((10,10)) # create an array of random numbers between 0 and 1

A_ceil = np.ceil(A)
A_floor = np.floor(A)
A_round = np.round(A)

figure(figsize=(12,16))
subplot(1,4,1)
imshow(A)
subplot(1,4,2)
imshow(A_ceil, vmin=0, vmax=1)
subplot(1,4,3)
imshow(A_floor)
subplot(1,4,4)
imshow(A_round)

# np.dot - dot product.  works for vectors and matrices - for 2d matricies it is the matrix product
A = [[1, 0], [0, 1]]
B = [[4, 1], [2, 2]]
print np.dot(A, B)

# np.exp(), np.log(), np.power(), calculate exponential, log, and power functions on an element by element basis

A = np.arange(10)
print A
print
print np.log(A)
print
print np.power(A, 2) # second argument is the power function -- can be another array!
print 
print np.exp(A+1)

# np.diff - discrete difference
# note: has an axis argument

A = np.arange(25)
A[12:] += 10
plot(A, label='data')
plot(np.diff(A), label='first difference')
legend()

A = np.arange(10,0,-1) # make a reverse list
print A[[1,5,7]]

A = np.arange(100).reshape(10,10)
index = A > 30
print A[index]

print A[index].shape

imshow(A)
figure()
imshow(index)



A = np.random.random((100,100,100)) * np.arange(100) # 100 frames of increasing random numbers
figure()
subplot(1,3,1)
imshow(A[:,:,0], vmax=100)
subplot(1,3,2)
imshow(A[:,:,50], vmax=100)
subplot(1,3,3)
imshow(A[:,:,99], vmax=100)

index = np.zeros((100,100), dtype=bool)
index[50:65,50:60] = True
print A[index,:].shape


figure()
imshow(index)
figure()
plot(A[index,:].mean(axis=0))

# np.random.random() - takes a tuple that is the shape of the array (like np.zeros()) and files with with random floats on the interval [0.0, 1.0)
# note that if you call it with no arguments, it returns a single number, of type 'float'

A = np.random.random(5)
B = np.random.random((5,5))

print A
print
print B

# np.random.seed()
for i in range(10):
    np.random.seed(0)
    print np.random.random()

# np.random.randint() is different- it returns a specified number of random integers in your specified interval
for i in range(10):
    print np.random.randint(0,7) # dice roll! (note that the high value is exclusive)  # randomintegers() is inclusive
    
print
# if you want a collection of them, you use the 'size' keyword
print np.random.randint(0,7,10)

A = np.arange(15)
print A
print
np.random.shuffle(A)
print A

# np.random.normal() - takes mean, std, and # of samples
for i in [10,100,1000,10000,100000]:
    figure()
    plt.hist(np.random.normal(5, 5, i))


# solve for Ax = b  
A = np.array([[1,2], [-1,1]])  
b = np.array([1,1])  
x = np.linalg.solve(A,b)  

print x

A = np.random.randint(0, 10, 25).reshape(5, 5)
print A
e_vals, e_vecs = np.linalg.eig(A)
print
print e_vals
print e_vecs

# example stolen from http://glowingpython.blogspot.com/2011/08/how-to-plot-frequency-spectrum-with.html

figure(figsize=(12,8))
def plotSpectrum(y,Fs):
 """
 Plots a Single-Sided Amplitude Spectrum of y(t)
 """
 n = len(y) # length of the signal
 k = np.arange(n)
 T = n/Fs
 frq = k/T # two sides frequency range
 frq = frq[range(n/2)] # one side frequency range

 Y = np.fft.fft(y)/n # fft computing and normalization
 Y = Y[range(n/2)]
 
 plot(frq,abs(Y),'r') # plotting the spectrum
 xlabel('Freq (Hz)')
 ylabel('|Y(freq)|')

Fs = 150.0;  # sampling rate
Ts = 1.0/Fs; # sampling interval
t = np.arange(0,1,Ts) # time vector

ff = 5;   # frequency of the signal
y = np.sin(2*np.pi*ff*t)

subplot(2,1,1)
plot(t,y)
xlabel('Time')
ylabel('Amplitude')
subplot(2,1,2)
plotSpectrum(y,Fs)

# polynomials have their own convience classes, you pass in the coefficents.
# the degree is infered from the length of the coefficent list.

p = np.polynomial.Polynomial([5,1,1]) 

# we can evaluate at a point or over a range of values
print p(3)
plot(np.arange(0,10,0.01), p(np.arange(0,10,0.01)))
xlabel('x')
ylabel('y')

# let's cheat and use the polynomial function to make some data for us
x_vals = np.arange(0,10,0.01)
p = np.polynomial.Polynomial([5,1,1])
clean_data = p(x_vals)

num_points = clean_data.shape[0]
dirty_data = clean_data + (np.random.random(num_points) * 10 - 5) # now dirty_data is clean data +/- 5 

fit_coefs = np.polynomial.polynomial.polyfit(x_vals, dirty_data, 2) # takes x values, y values and the degree to try and fit the data to
print fit_coefs
fit_p = np.polynomial.Polynomial(fit_coefs) # now build a polynomial with those fit values

# plot them
plot(x_vals, dirty_data, label='dirty data')
plot(x_vals, fit_p(x_vals), lw=2, color='black', label='fit')
legend()

A = np.random.random(100)
A_masked_small = np.ma.MaskedArray(A, A>0.5) # second arguement is the mask, must be the same shape as the array
A_masked_big = np.ma.MaskedArray(A, A<=0.5) # second arguement is the mask, must be the same shape as the array


print 'all average', A.mean()
print 'big average', A_masked_big.mean()
print 'small average', A_masked_small.mean()

figure(figsize=(12,8))

plot(A, label='full')
plot(A_masked_small, label='small')
plot(A_masked_big, label='big')

legend()