import numpy as np

"""
Explanation: NumPy: Vectorized Array Processing in Python
NumPy, short for Numerical Python, is the fundamental package required for high performance scientific computing and data analysis. It is the foundation on which nearly all of the higher-level tools we will use are built. Here are some of the things it provides:

ndarray, a fast and space-efficient multidimensional array providing vectorized arithmetic operations and sophisticated broadcasting capabilities. We use "array" to refer to ndarrays colloquially.
Standard mathematical functions for fast operations on entire arrays of data without having to write loops
Tools for reading / writing array data to disk and working with memory-mapped files
Linear algebra, random number generation, and Fourier transform capabilities
Tools for integrating code written in C, C++, and Fortran

The standard import line for NumPy is the following:
End of explanation
"""


my_list = [1, 2, 3, 4, 5, 6, 7, 8, 9] # equivalently, range(1,10)
my_array = np.array(my_list)

"""
Explanation: Arrays
NumPy arrays work fundamentally differently than the standard Python list, though they can be initialized by wrapping around a Python list, as below.
End of explanation
"""


print my_array[0] # returns the first element of my_array
print my_array[3:5] # returns the fourth and fifth elements of my_array
print my_array[-2] # returns the second-to-last element of my_array

"""
Explanation: Arrays are indexed just like in standard Pyhon. For example,
End of explanation
"""


print my_list + my_list

print my_array + my_array

print my_array * 8

"""
Explanation: However, there are some immediate differences between Python lists and NumPy arrays. For example, consider these lines of code
End of explanation
"""


np.arange(0,10,0.4)

"""
Explanation: It looks as though NumPy arrays handle arithmetic operations differently than Python lists! Stay tuned to learn more about this.
Initializing NumPy arrays
There are four main ways to initialize an array in NumPy. The first is by doing what was shown above, wrapping an NumPy array around an existing Python list. A second common initialization is by using arange, which takes a starting point, an ending point (which is excluded), and an increment value, and generates the array. For example:
End of explanation
"""


np.zeros((2,4)) # note that np.zeros takes a tuple (m,n)

"""
Explanation: This is similar to Python's range function for lists. Another common way to initialize arrays is to use the zeros function, which generates an array of zeros for a given size. For example:
End of explanation
"""


np.ones((3,1)) # don't forget the tuple!

"""
Explanation: Finally, sometimes you want to initialize an array to all 1s, so NumPy provides the ones function:
End of explanation
"""


fidentity = np.eye (5, dtype='float16')
bidentity = np.eye (5, dtype='bool')
sidentity = np.eye (5, dtype='string_')

print fidentity, "\n\n"
print bidentity, "\n\n"
print sidentity

"""
Explanation: The full list of array creation functions are given in the following table:
Array creation functions
array
Convert input data (list, tuple, array, or other sequence type) to an ndarray either by inferring a dtype or explicitly specifying a dtype. Copies the input data by default 
asarray
Convert input to ndarray, but do not copy if the input is already an ndarray 
arange
Like the built-in range but returns an ndarray instead of a list
ones, ones_like
Produce an array of all 1's with the given shape and dtype.ones_like takes another array and produces a ones array of the same shape and dtype.
zeros, zeros_like
Like ones and ones_like but producing arrays of 0's instead
empty, empty_like
Create new arrays by allocating new memory, but do not populate with any values like ones and zeros
eye, identity
Create a square NxN identity matrix (1's on the diagonal and 0's elsewhere)
Be careful using empty, because it stores the ndarray with whatever happens to be stored on the main stack. It's almost always preferable to initialize an "empty" array with zeros, because it is guaranteed what the values will be.
Data types
If necessary, you can specifiy the data type in the ndarray by including a dtype= parameter. The acceptable dtypes are:
int8, uint8 (i1, u1)
Signed and unsigned 8-bit (1 byte) integer types
int16, uint16 (i2, u2)
Signed and unsigned 16-bit integer types
int32, uint32 (i4, u4)
Signed and unsigned 32-bit integer types
int64, uint64 (i8, u8)
Signed and unsigned 64-bit integer types
float16 (f2)
Half-precision floating point
float32 (f4 or f)
Standard single-precision floating point. Compatibile with C float
float64 (f8 or d)
Standard double-precision floating point. Compatibile with C double and Python float object
float128 (f16 or g)
Extended floating point
complex64, complex128, complex256 (c8, c16, c32)
Complex numbers represented by 32, 64, or 128 floats, respectively
bool  (?)
Boolean type storing True and False values
object (O)
Python object type
string_  (S)
Fixed-length string type (1 byte per character). For example, to create a string dtype with length 10, use S10
unicode_ (U)
Fixed-length unicdoe type (number of bytes platform specific). Same specification semantics as string_ (e.g. U10).
End of explanation
"""


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

"""
Explanation: Accessing array elements
NumPy arrays can be accessed in similar fashion to Python list types. For example, given an array
End of explanation
"""


print a[0]
print a[:1]
print a[::-1]

"""
Explanation: we can access elements in the natural sense:
End of explanation
"""


a = np.array([[1., 2.], [3., 4.]])
a[0]

"""
Explanation: There are two ways to access elements for multidimensional arrays. The most natural accessing syntax originates from the following observation:
End of explanation
"""


print a[0][0], a[0][1]

"""
Explanation: In a 2x2 matrix, accessing the first element returns the first row vector of the matrix. Thus, we can treat a[0] as its own array, so we can access any of its elements by using the format:
End of explanation
"""


print a[0,0], a[0,1]

"""
Explanation: This is called the recursive approach. There is another approach, which is unique to NumPy arrays, which lets you overload the arguments in the square brackets:
End of explanation
"""


myarr = np.arange(1.,11.,1.)
print myarr * 2.
print myarr - 6.
print myarr * myarr
print np.log(myarr)

"""
Explanation: Exercise
Take a large array, and use %timeit to determine which accessor approach is more efficient.
The moral of the story is, use the NumPy specific approach!
Operations between ndarrays and scalars
Arrays are important because they enable you to express batch operations on data without writing any for loops. This is usually called vectorization. Any arithmetic operations between equal-size arrays applies the operation elementwise.
Try it!
Run the following code to get a sense of how ndarray operations with scalars work.
End of explanation
"""


print my_array
print my_array[::2] # the sub-array starting at 0 to the end with even index

"""
Explanation: Exercise
Write a function that takes a tuple $(m,n)$ and a number $x$ generates an array of containing the value $x$.
Note that you can extract a sub-array from a main array not only by slicing sections but by choosing iterations. For example,
End of explanation
"""


np.random.rand()

"""
Explanation: Write a function that declares an $n\times n$ "checkerboard" array of 0s and 1s.
Lab: Invertible Matrices
NumPy provides powerful methods for randomly-generated arrays using the np.random class. The simplest random-number generation task is to produce a random number $x$ that is a member of the set $[0,1)$; that is, $0 \leq x < 1$. In NumPy, this is achieved by:
End of explanation
"""


np.random.rand(5)

"""
Explanation: Notice that if you run this cell many times, you get a different result. np.random.rand() takes $n$ arguments which determine the shape of the output. For example, to get a 5-dimensional random vector, we can write
End of explanation
"""


Z = np.array([[0,0,0,0,0,0],
                  [0,0,0,1,0,0],
                  [0,1,0,1,0,0],
                  [0,0,1,1,0,0],
                  [0,0,0,0,0,0],
                  [0,0,0,0,0,0]],dtype='int64')

"""
Explanation: Each additional argument specifies another "axis" of the array output. So if you give two arguments, it produces a matrix; three arguments produces a "cube" matrix; and $n$ arguments produces an $n$-tensor.
The class np.random has many more random array capabilities which you can find here. We will just use np.random.rand for this lab.
Exercise

Generate a three-dimensional random vector whose entries range between 4 and 8.
Generalize this process into a function that produces an $n$-dimensional random vector whose components are elments of $[a,b)$. 

Main lab
Consider the following problem: given an $n\times n$-dimension matrix $A$, what is the probability that $A$ is invertible? One reasonable interpretation of this problem is to ask how frequently a randomly-generated matrix is invertible. Your task is to

write a function that generates a random $n\times n$ matrix whose values fall in the set $[-1,1)$,
write a function to determine if such a matrix is invertible,
write a function to approximate the probability of such an event by considering $N$ matrices and by recording the number of these that are invertible.

Recall that a matrix $A$ is invertible if and only if $\det(A)\ne 0$. To this end, you might want to check out linalg, NumPy's linear algebra class.
Also, because we are dealing with numerical precision, it probably is best not to expect that the determinant will ever actually be 0. Instead, we suggest including an error tolerance.
It turns out that the overwhelming majority of random matrices are invertible. This simulation supports the theory, which you can read about on this StackExchange post. Hopefully, you're beginning to see the power of programming in problem-solving.
Lab: Cellular Automata
(From Nicolas P. Rougier's Numpy tutorial)
We will construct a simulation of John Conway's Game of Life using the skills we have learned about NumPy.

The Game of Life, also known simply as Life, is a cellular automaton devised by the British mathematician John Horton Conway in 1970. It is the best-known example of a cellular automaton. The "game" is actually a zero-player game, meaning that its evolution is determined by its initial state, needing no input from human players. One interacts with the Game of Life by creating an initial configuration and observing how it evolves.
The universe of the Game of Life is an infinite two-dimensional orthogonal grid of square cells, each of which is in one of two possible states, live or dead. Every cell interacts with its eight neighbours, which are the cells that are directly horizontally, vertically, or diagonally adjacent. At each step in time, the following transitions occur:

Any live cell with fewer than two live neighbours dies, as if by needs caused by underpopulation.
Any live cell with more than three live neighbours dies, as if by overcrowding.
Any live cell with two or three live neighbours lives, unchanged, to the next generation.
Any dead cell with exactly three live neighbours becomes a live cell.

The initial pattern constitutes the 'seed' of the system. The first generation is created by applying the above rules simultaneously to every cell in the seed – births and deaths happen simultaneously, and the discrete moment at which this happens is sometimes called a tick. (In other words, each generation is a pure function of the one before.) The rules continue to be applied repeatedly to create further generations.

Getting started
The first thing to do is to create the proper numpy array to hold the cells. This can be done very easily:
End of explanation
"""


Z.dtype

"""
Explanation: You don't have to specify the data type, as NumPy will interpret an array of integers as an int64 data type anyways. Sometimes, though, it's good practice to be specific. You can always check what datatype an array is by running
End of explanation
"""


Z.shape

"""
Explanation: We can also check the shape of the array to make sure it is 6x6:
End of explanation
"""


A = Z[1:5,1:5]

"""
Explanation: You already know how to access elements of Z. Write a statement to obtain the the element in the 3rd row and the 4th column.
Better yet, we can access a subpart of the array using the slice notation. Obtain the subarray of Z containing the rows 2-5 and columns 2-5, and set it equal to an ndarray labelled A.
End of explanation
"""


A[0, 0] = 9
print A, "\n\n"
print Z

"""
Explanation: Be mindful of array pointers! Look at what happens when you run the following code:
End of explanation
"""


print Z.base is None
print A.base is Z

A[0,0] = 0 # put A (and Z) back to normal.

"""
Explanation: We set the value of A[0,0] to 9 and we see immediate change in Z[1,1] because A[0,0] actually corresponds to Z[1,1]. This may seem trivial with such simple arrays, but things can become much more complex (we'll see that later). If in doubt, you can use 
ndarray.base 
to check easily if an array is part of another one:
End of explanation
"""


1 + (2*Z + 3)

"""
Explanation: Counting neighbors
We now need a function to count the neighbours. We could do it the same way as for the python version, but this would make things very slow because of the nested loops. We would prefer to act on the whole array whenever possible, this is called vectorization.
First, you need to know that you can manipulate Z as if (and only as if) it was a regular scalar:
End of explanation
"""


Z + Z[-1:1, -1:1]

"""
Explanation: If you look carefully at the output, you may realize that the ouptut corresponds to the formula above applied individually to each element. Said differently, we have 
(1+(2*Z+3))[i,j] == (1+(2*Z[i,j]+3)) 
for any i,j.
Ok, so far, so good. Now what happens if we add Z with one of its subpart, let's say Z[-1:1,-1:1] ?
End of explanation
"""


Z + 1

"""
Explanation: This raises a Value Error but more interestingly, numpy complains about the impossibility of broadcasting the two arrays together. Broadcasting is a very powerful feature of numpy and most of the time, it saves you a lot of hassle. Let's consider for example the following code:
End of explanation
"""


N = np.zeros(Z.shape, dtype=int)

N[1:-1,1:-1] += (Z[ :-2, :-2] + Z[ :-2,1:-1] + Z[ :-2,2:] +
                     Z[1:-1, :-2]                + Z[1:-1,2:] +
                     Z[2:  , :-2] + Z[2:  ,1:-1] + Z[2:  ,2:])

N

"""
Explanation: How can a matrix and a scalar be added together? Well, they can't. But NumPy is smart enough to guess that you actually want to add 1 to each of the element of Z. This concept of broadcasting is quite powerful and it will take you some time before masterizing it fully (if even possible).
However, in the present case (counting neighbours if you remember), we won't use broadcasting (uh?). But we'll use vectorize computation using the following code:
End of explanation
"""


def iterate(Z):
    # Iterate the game of life : naive version
    # Count neighbours
    N = np.zeros(Z.shape, int)
    N[1:-1,1:-1] += (Z[0:-2,0:-2] + Z[0:-2,1:-1] + Z[0:-2,2:] +
                     Z[1:-1,0:-2]                + Z[1:-1,2:] +
                     Z[2:  ,0:-2] + Z[2:  ,1:-1] + Z[2:  ,2:])
    N_ = N.ravel()
    Z_ = Z.ravel()

    # Apply rules
    R1 = np.argwhere( (Z_==1) & (N_ < 2) )
    R2 = np.argwhere( (Z_==1) & (N_ > 3) )
    R3 = np.argwhere( (Z_==1) & ((N_==2) | (N_==3)) )
    R4 = np.argwhere( (Z_==0) & (N_==3) )

    # Set new values
    Z_[R1] = 0
    Z_[R2] = 0
    Z_[R3] = Z_[R3]
    Z_[R4] = 1

    # Make sure borders stay null
    Z[0,:] = Z[-1,:] = Z[:,0] = Z[:,-1] = 0

"""
Explanation: To understand this code, have a look at the figure below:

What we actually did with the above code is to add all the darker blue squares together. Since they have been chosen carefully, the result will be exactly what we expected. If you want to convince yourself, consider a cell in the lighter blue area of the central sub-figure and check what will the result for a given cell.
Cell Iteration
In a first approach, we can write the iterate function using the argwhere method that will give us the indices where a given condition is True.
End of explanation
"""


def iterate_2(Z):
    # Count neighbours
    N = (Z[0:-2,0:-2] + Z[0:-2,1:-1] + Z[0:-2,2:] +
         Z[1:-1,0:-2]                + Z[1:-1,2:] +
         Z[2:  ,0:-2] + Z[2:  ,1:-1] + Z[2:  ,2:])

    # Apply rules
    birth = (N==3) & (Z[1:-1,1:-1]==0)
    survive = ((N==2) | (N==3)) & (Z[1:-1,1:-1]==1)
    Z[...] = 0
    Z[1:-1,1:-1][birth | survive] = 1
    return Z

"""
Explanation: Even if this first version does not use nested loops, it is far from optimal because of the use of the 4 argwhere calls that may be quite slow. We can instead take advantages of numpy features the following way.
End of explanation
"""


import matplotlib.pyplot as plt
from matplotlib import animation
from IPython.display import Image
%matplotlib inline

"""
Explanation: Now, let's throw together a simple visualization for our Game of Life! Don't worry about the matplotlib code, we'll deal with it later.
End of explanation
"""


Z = np.random.randint(0,2,(256,256))
fig = plt.figure()
ax = plt.axes()

im = plt.imshow(Z, cmap='gray', interpolation='bicubic')

def updatefig(*args):
    im.set_array(iterate_2(Z))
    return im,
    
ani = animation.FuncAnimation(fig, updatefig, interval=200, blit=True)

ani.save('demonstration.gif', writer='imagemagick', fps=10)

"""
Explanation: The cell below saves our animated Game of Life as a .gif format on your computer.
End of explanation
"""


Image(url='demonstration.gif')

"""
Explanation: And now, load up the .gif file with:
End of explanation
"""


