#########################################################################################
#
# Scientific Python Course
#
# L02 - NumPy
#
#########################################################################################
#
# ARCHER, 2015
#
# http://www.archer.ac.uk
# support@archer.ac.uk
#
# L02 - NumPy
#
# Neelofer Banglawala : nbanglaw@epcc.ed.ac.uk
# Kevin Stratford     : kevin@epcc.ed.ac.uk
#
# EPCC, University of Ediburgh, EPSRC, NERC, CRAY
#
#
# ********************************************************************
# * Reusing this material :                                          *
# *                                                                  *
# * This work is licensed under a Creative Commons Attribution-      *
# * Non-Commercial-ShareAlike 4.0 Internatiional License.            *
# * http://creativecommons.org/licenses/by-nc-sa/4.0/deed.en_US      *
# *                                                                  *
# * This means you are free to copy and redistributhe the material   *
# * and adapt and build on the material under the following terms:   *
# * You must give appropriate credit, provide a link to the license  *
# * and indicate if changes were made. If you adapt or build on the  *
# * material you must distribute your work under the same license as *
# * the original.                                                    *
# *                                                                  *
# ********************************************************************
#
#
# 
#########################################################################################
#  
#   S5 [NumPy] Calculating pi II
#  
#########################################################################################
# 
# # calculate an approximation to pi
# # will understand this code by the 
# # end of the session!
# import numpy as np
# 
# # N : number of iterations
# def calc_pi(N):
#     x = np.random.ranf(N);
#     y = np.random.ranf(N);
#     r = np.sqrt(x*x + y*y);
#     c=r[ r <= 1.0 ]
#     return 4*float((c.size))/float(N)
# 
# # time the results
# pts = 6; N = np.logspace(1,8,num=pts);
# result = np.zeros(pts); count = 0;
# for n in N:
#     result = %timeit -o -n1 calc_pi(n) 
#     result[count] = result.best
#     count += 1
#     
# # and save results to file   
# np.savetxt('calcpi_timings.txt', np.c_[N,results],
#            fmt='%1.4e %1.6e');
# 
# 
#########################################################################################
#  
#   S7 [NumPy] Creating arrays I
#  
#########################################################################################
# 
# # import numpy as alias np
# import numpy as np
# 
# # create a scalar (zero dimensional array)
# a = np.array( 42 ); a
# 
#########################################################################################
#  
#   S7-1 [NumPy] Creating arrays II
#  
#########################################################################################
# 
# # create a 1d array with a list
# a = np.array( [-1,0,1] ); a
# 
# # use functions that create lists e.g. range()
# a = np.array( range(-2,6,2) ); a
# 
# # use arrays to create arrays
# b = np.array( a ); b
# 
# # use numpy functions to create arrays 
# # arange for arrays, range for lists!
# a= np.arange( -2, 6, 2 ); a
# 
# 
#########################################################################################
#  
#   S7-2 [NumPy] Creating arrays III
#  
#########################################################################################
# 
# # between start, stop, sample step points
# a = np.linspace(-10,10,5); 
# a;
# 
# # Ex: can you guess these functions do?
# b = np.zeros(3); print b
# c = np.ones(3); print c
# 
# # Ex++: try these and see what you get 
# h = np.hstack( (a, a, a) ); print h
# o = np.ones_like(a); print  o
# 
#########################################################################################
#  
#   S8 [NumPy] Array characteristics
#  
#########################################################################################
# 
# # array characteristics such as:
# print a
# #print a.ndim  # dimensions
# #print a.shape # shape
# #print a.size  # size
# #print a.dtype # data type
# 
# # can choose data type
# a = np.array( [1,2,3], np.int16 ); a.dtype
# 
# # Ex: query the characteristics of the arrays you have created
# 
#########################################################################################
#  
#   S9 [NumPy] Multi-dimensional arrays I
#  
#########################################################################################
# 
# # multi-dimensional arrays e.g. 2d array or matrix 
# # e.g. list of lists
# mat = np.array( [[1,2,3], [4,5,6]]);
# print mat; print mat.size; mat.shape
# 
# # can create 2d arrays with complex elements (e.g.1 + 2*i)
# lst = [[1, 2, 3], [4,5,6]];
# mat = np.array(lst, complex);
# print mat; print mat.size; print mat.shape
# 
# # join arrays along first axis (0)
# d=np.r_[np.array([1,2,3]), 0, 0, [4,5,6]];
# print d; d.shape
# 
#########################################################################################
#  
#   S9-1 [NumPy] Multi-dimensional arrays II
#  
#########################################################################################
# 
# # join arrays along second axis (1)
# d=np.c_[np.array([1,2,3]), [4,5,6]];
# print d; d.shape
# 
# # Ex: use r_, c_ with nd (n>1) arrays
# 
# # Ex: can you guess the shape of these arrays?
# h = np.array( [1,2,3,4,5,6] ); print h.shape
# i = np.array( [[1,1],[2,2],[3,3],[4,4],[5,5],[6,6]] ); 
# print i.shape
# j = np.array( [[[1],[2],[3],[4],[5],[6]]] ); print j.shape
# k = np.array( [[[[1],[2],[3],[4],[5],[6]]]] ); print k.shape
# 
# 
#########################################################################################
#  
#   S10 [NumPy] Reshaping arrays I
#  
#########################################################################################
# 
# # reshape 1d arrays into nd arrays original matrix unaffected
# mat = np.arange(6); print mat
# print mat.reshape( (3, 2) )
# print mat; print mat.size;
# print mat.shape
# 
# mat = np.resize(mat, (1,mat.size)); 
# mat.shape
# 
# # can also use the shape, this modifies the original array
# a = np.zeros(10); 
# print a
# a.shape = (2,5)
# print a; print a.shape;
# 
# 
#########################################################################################
#  
#   S10-1 [NumPy] Reshaping arrays II
#  
#########################################################################################
# 
# # use flatten() or ravel() to go from nd to 1d 
# # this creates a COPY of the original 
# mat = np.array([[1,2,3],[4,5,6]])
# mat2 = mat.flatten()
# print mat2; print mat2.size; mat2.shape
# 
# # unlike shape, original matrix unaffected
# print mat; print mat.size; mat.shape
# 
# # Ex: split a martix? Change the cuts and axis values
# # need help?: np.split?
# cuts=2;
# np.split(mat, cuts, axis=0)
# 
#########################################################################################
#  
#   S10-2 [NumPy] More array functions I
#  
#########################################################################################
# 
# # use copyto to copy values
# # to array
# a = np.array( [-2,6,2] ); print a
# b = np.ones(3); print b
# np.copyto(b, a); print b
# 
# # Ex: create some nd arrays from the methods we've seen
# # query their characteristics, have a play
# 
# # Ex++: can you guess what these functions do?
# v = np.vstack( (arr2d, arr2d) ); print v; v.ndim;
# c0 = np.concatenate( (arr2d, arr2d), axis=0); c0;
# 
# # Ex++: can you guess what this will do?
# c1 = np.concatenate(( mat, mat ), axis=1); print "c1:", c1;
# 
#########################################################################################
#  
#   S10-3 [NumPy] More array functions II
#  
#########################################################################################
# 
# # Ex++: other functions to explore
# #
# # stack(arrays[, axis])
# # tile(A, reps)
# # repeat(a, repeats[, axis])
# # unique(ar[, return_index, return_inverse, ...])
# # trim_zeros(filt[, trim]), fill(scalar)
# # xv, yv = meshgrid(x,y)
# 
#########################################################################################
#  
#   S11 [NumPy] Accessing arrays I
#  
#########################################################################################
# 
# # basic indexing and slicing we know from lists
# # a[start:stop:step] --> [start, stop every step)
# a = np.arange(8); print a
# print a[0:7:2]
# print a[0::2]
# 
# # negative indices are valid!
# print a[2:-3:2]
# 
#########################################################################################
#  
#   S11-1 [NumPy] Accessing arrays II
#  
#########################################################################################
# 
# # basic indexing of a 2d array :take care of each dimension
# nd = np.arange(12).reshape((4,3)); print nd;
# print nd[(2,2)];
# print nd[2][2];
# 
# # get corner elements 0,2,9,11
# print nd[0:4:3, 0:3:2]
# 
# # Ex: get elements 7,8,10,11 that make up the bottom right corner
# nd = np.arange(12).reshape((4,3)); 
# print nd; nd[2:4, 1:3]
# 
#########################################################################################
#  
#   S11-2 [NumPy] Slices and copies I
#  
#########################################################################################
# 
# # slices are views (like references) 
# # on array, can change elements
# nd[2:4, 1:3] = -1; nd
# 
# # assign slice to a variable to prevent this
# s = nd[2:4, 1:3]; print nd; 
# s = -1; nd
# 
#########################################################################################
#  
#   S11-3 [NumPy] Slices and copies II
#  
#########################################################################################
# 
# # simple assignment creates references, 
# nd = np.arange(12).reshape((4,3))
# md = nd
# md[3] = 1000
# print nd
# 
# # can avoid this by using copy()
# nd = np.arange(12).reshape((4,3))
# md = nd.copy()
# md[3]=999
# print nd
# 
#########################################################################################
#  
#   S11-4 [NumPy] Fancy indexing I
#  
#########################################################################################
# 
# # advanced or fancy indexing lets you do more
# p = np.array([[ 0,  1,  2],[ 3,  4,  5],[ 6,  7,  8],[ 9, 10, 11]]); 
# print p
# 
# rows = [0,0,3,3]; cols = [0,2,0,2];
# print p[rows,cols]
# 
# # Ex: what will this slice look like?
# m = np.array([[0,-1,4,20,99],[-3,-5,6,7,-10]]);
# print m[[0,1,1,1],[1,0,1,4]];
# 
# 
#########################################################################################
#  
#   S11-5 [NumPy] Fancy Indexing II
#  
#########################################################################################
# 
# # can use conditionals in indexing
# # m = np.array([[0,-1,4,20,99],[-3,-5,6,7,-10]]);
# m[ m < 0 ]
# 
# # Ex: can you guess what this does? query: np.sum?
# y = np.array([[0, 1], [1, 1], [2, 2]]);
# rowsum = y.sum(1);
# y[rowsum <= 2, :]
# 
# # Ex: and this? 
# a = np.arange(10);
# mask = np.ones(len(a), dtype=bool);
# mask[[0,2,4]] = False; print mask
# result = a[mask]; result
# 
#########################################################################################
#  
#   S12 [NumPy] Manipulating arrays
#  
#########################################################################################
# 
# # add an element with insert
# a = np.arange(6).reshape([2,3]); print a
# np.append(a,np.ones([2,3]),axis=0)
# 
# # inserting an array of elements
# np.insert(a, 1, -10, axis=0)
# 
# # can use delete, or a boolean mask, to delete array elements
# a = np.arange(10)
# np.delete(a, [0,2,4], axis=0)
# 
#########################################################################################
#  
#   S13 [NumPy] Vectorization
#  
#########################################################################################
# 
# # vectorization allows element-wise operations (no for loop!)
# a = np.arange(10).reshape([2,5]); b = np.arange(10).reshape([2,5]);
# 
# -0.1*a
# 
# a*b
# 
# a/(b+1)  #.astype(float)
# 
#########################################################################################
#  
#   S13-1 [NumPy] Random number generation
#  
#########################################################################################
# 
# # random floats
# a = np.random.ranf(10); a 
# 
# 
# # create random 2d int array 
# a = np.random.randint(0,high=5,size=25).reshape(5,5); 
# print a;
# 
# 
# # generate sample from normal distribution 
# # (mean=0, standard deviation=1)
# s = np.random.standard_normal((5,5)); s;
# 
# # Ex: what other ways are there to generate random numbers?
# # What other distributions can you sample? 
# 
#########################################################################################
#  
#   S14 [NumPy] File IO
#  
#########################################################################################
# 
# # easy way to save data to text file
# pts=5; x=np.arange(pts); y=np.random.random(pts);
# 
# # format specifiers: d = int, f = float, e = exponential
# np.savetxt('savedata.txt', np.c_[x,y],header='DATA', footer='END',
#            fmt='%d %1.4f')
# 
# !cat savedata.txt
# #p=np.loadtxt('savedata.txt')
# 
# # much more flexibility with genfromtext
# p = np.genfromtxt('savedata.txt', skip_header=2,skip_footer=1); p
# 
# # Ex++: what do numpy.save, numpy.load do ?
# 
# 
#########################################################################################
#  
#   S15-2 [NumPy] Polynomials : calculating pi II
#  
#########################################################################################
# 
# # calculate pi using polynomials
# # import Polynomial class 
# from numpy.polynomial import Polynomial as poly;
# num = 100000;
# denominator = np.arange(num); 
# 
# denominator[3::4] *=-1 # every other odd coefficient is -ve
# numerator = np.ones(denominator.size); 
# 
# # avoid dividing by zero, drop first element denominator
# almost=numerator[1:]/denominator[1:]; 
# 
# # make even coefficients zero
# almost[1::2] = 0
# 
# # add back zero coefficient
# coeffs = np.r_[0,almost]; 
# 
# p=poly(coeffs); 4*p(1) # pi approximation
# 
#########################################################################################
#  
#   S16-1 [NumPy] Performance II
#  
#########################################################################################
# 
# # accessing a 2d array
# nd = np.arange(100).reshape((10,10))
# 
# # accessing element of 2d array
# %timeit -n10000000 -r3 nd[5][5]
# %timeit -n10000000 -r3 nd[(5,5)]
# 
# # Ex: multiplying two vectors
# x=np.arange(10E7)
# %timeit -n1 -r10 x*x
# %timeit -n1 -r10 x**2
# 
# # Ex++: from the linear algebra package
# %timeit -n1 -r10 np.dot(x,x) 
# 
#########################################################################################
#  
#   S16-2 [NumPy] Performance III
#  
#########################################################################################
# 
# import numpy as np
# # Ex: range functions and iterating
# # in  for loops
# size = int(1E6)
# 
# %timeit for x in range(size): x ** 2
# 
# # faster than range for very large arrays
# %timeit for x in xrange(size): x ** 2
# 
# %timeit for x in np.arange(size): x ** 2
# 
# %timeit np.arange(size) ** 2
# 
# # Ex: look over the two ways of calculating pi. 
# # Make sure you understand the code. 
# # Time each method, which is faster?
# 
#########################################################################################
#  
#  END
#  
#########################################################################################