%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp

"""
Explanation: <h1 align="center">Scientific Programming in Python</h1>
<h2 align="center">Topic 2: NumPy and Efficient Numerical Programming</h2>

Notebook created by Martín Villanueva - martin.villanueva@usm.cl - DI UTFSM - April 2017.
End of explanation
"""


%time {1+2 for i in range(10000)}

"""
Explanation: Table of Contents

1.- Useful Magics
2.- Basic NumPy Operations
3.- Internals of NumPy
4.- Efficient Programming in NumPy
5.- Vectorized Algorithms in NumPy
6.- Bonus 

<div id='magics' />
Useful magic built-in functions for time-profiling

time: See how long it takes a code to run.
End of explanation
"""


%timeit [1+2 for i in range(10000)]

%timeit -n 100 [1+2 for i in range(10000)]

%%timeit
for i in range(10000):
    1+2

"""
Explanation: timeit: 
See how long a code takes to run averaged over multiple runs.  
It will limit the number of runs depending on how long the script takes to execute.  
Provide an accurate time calculation by reducing the impact of startup or shutdown costs on the time calculation by executing the code repeatedly.
End of explanation
"""


from time import sleep

def foo(): sleep(1)
    
def bar(): sleep(2)
    
def baz(): foo(),bar()

%prun baz()

"""
Explanation: prun: See how long it took each function in a script to run.
End of explanation
"""


# Arrays of zeros: np.zeros(shape)
print("Zeros:")
print( np.zeros((3,3)) )

# Arrays of ones: np.ones(shape)
print("\nOnes:")
print( np.ones((3,3)) )

# Empty array: np.empty(shape)
print("\nEmpty:")
print( np.empty((3,3)) )

# Range of values: np.range(start, stop, step)
print("\nRange:")
print( np.arange(0., 10., 1.) )

# Regular grid: np.linspace(start, end, n_values)
print("\nRegular grid:")
print( np.linspace(0., 1., 9) )

# Random secuences: np.random
print("\nRandom sequences:")
print( np.random.uniform(10, size=6) )

# Array constructor: np.array( python_iterable )
print("\nArray constructor")
print( np.array([2, 3, 5, 10, -1]) )


print( 10*np.random.random((5,5)) )

"""
Explanation: <div id='basic_numpy' />
Basic NumPy operations
The reasons of why you should use NumPy instead of any other iterable object_ in Python are:
* NumPy provides an ndarray structure for storing numerical data __in a contiguous way.
* Also implements fast mathematical operations on ndarrays, that exploit this contiguity.
* Brevity of the syntax for array operations. A language like C or Java would require us to write a loop for a matrix operation as simple as C=A+B.
Creating Arrays
There are several NumPy functions for creating common types of arrays. Below is a list of the most common used:
End of explanation
"""


# first we create two random arrays:
A = np.random.random((5,5))
B = np.random.random((5,5))

# sum
print("Sum:")
print( A+B )

# subtraction
print("\nSubtraction")
print( A-B )

# product
print("\nProduct")
print( A*B )

# matricial product
print("\nMatricial Product")
print( np.dot(A,B) )

# power
print("\n Power")
print( A**2 )

# Some common mathematical functions
print("\n np.exp()")
print( np.exp(A) )
print("\n np.sin()")
print( np.sin(A) )
print("\n np.cos()")
print( np.cos(A))
print("\n np.tan()")
print( np.tan(A) )

"""
Explanation: Basic Mathematical operations
Most of the operations performed in NumPy are handled element-wise, i.e, computing C = A + B will translates into $C[i,j] = A[i,j] + B[i,j]$. (The exception is broadcasting, and will be explained soon).
Below is a list with the most common used mathematical operations. For a comprehensive list se here: NumPy mathematical functions.
End of explanation
"""


# Creating two 2d-arrays
A = np.array( [[1, 2, 3], [2, 3, 5], [1, 9, 6]] )
B = np.array( [[1, 2, 3], [3, 5, 5], [0, 8, 5]] )

print("A > B:")
print( A > B )

print("\nA =< B:")
print( A <= B )

print("\n A==B:")
print( A==B )

print("\n A!=B:")
print( A!=B )

# Creating two 2d boolean arrays
C = A==B
D = A>=B

print("\n A and B:")
print( C & D)
print( np.logical_and(C,D) )

print("\n A or B:")
print( C | D)
print( np.logical_or(C,D) )

print("\n not A:")
print( ~C )
print( np.logical_not(C))

"""
Explanation: Boolean operations
Comparisons in NumPy work exaclty the same way as mathematical operations, i.e, element wise!. Let's see some examples:
End of explanation
"""


# Lets create a random array
A = np.random.random((5,5))

print("Dims: ")
print(A.ndim)

print("\nShape: ")
print(A.shape)

print("\nStrides: ")
print(A.strides)

print("\nData type: ")
print(A.dtype)

"""
Explanation: <div id='internals' />
Internals of NumPy
The numpy.ndarray structure
The ndarray is the NumPy object that let us create $N$-dimensional arrays. It is essentially defined by:
1. A number of dimensions
2. a shape
3. strides
4. data type or dtpe
5. The data buffer.
<img src='data/ndarray.png' style="width: 500px;">
End of explanation
"""


# Lets create a random array in C-order
A = np.random.random((5,2))
print("C strides:")
print(A.strides)

# Lets create a random array in F-order
B = np.asfortranarray(A)
print("\nF strides:")
print(B.strides)

"""
Explanation: How an ndarray is stored in memory?
When there is more than one dimension, there are two ways of storing the elements in the memory block:
1. Elements can be stored in row-major order (also known as C-order) or,
2. In column-major order (also known as Fortran-order).
<img src='data/ndarray_storage.png' style="width: 800px;">
What are the strides?
NumPy uses the notion of strides to convert between a multidimensional index and the memory location of the underlying (1D) sequence of elements. 
For example, the mapping between array[i,j] and the corresponding address the byte is:
* offset = array.strides[0] * i1 + array.strides[1] * i2
* address = base + offset
where base is the address of the first byte (array[0,0]).
End of explanation
"""


X = np.random.random((5000,5000))

%timeit X[0,:].sum()

%timeit X[:,0].sum()

"""
Explanation: An interesting result
End of explanation
"""


#Python Lists implementation
def norm_square_list(vector):
    norm = 0
    for v in vector:
        norm += v*v
    return norm

#Naive NumPy implementation
def norm_square_array(vector):
    norm = 0
    for v in vector:
        norm += v*v
    return norm

#Vectorized NumPy implementation
def norm_square_numpy(vector):
    return np.sum(vector * vector)

#Clever NumPy implementation
def norm_square_dot(vector):
    return np.dot(vector, vector)

#Vector to use - dimension 10^6
vector = range(1000000)
npvector = np.array(vector)

#Timing the list implementation
%timeit norm_square_list(vector)

#Timing the naive array implementation
%timeit norm_square_array(npvector)

#Timing the NumPy-vectorized implementation
%timeit norm_square_numpy(npvector)

#Timing the clever NumPy-vectorized implementation
%timeit norm_square_dot(npvector)

"""
Explanation: Why Python is so slow?

 Python is Dynamically Typed rather than Statically Typed. What this means is that at the time the program executes, the interpreter doesn't know the type of the variables that are defined. <img src='data/cint_vs_pyint.png' style="width: 300px;">
Python is interpreted rather than compiled. A smart compiler can look ahead and optimize for repeated or unneeded operations, which can result in speed-ups. 
Python's object model can lead to inefficient memory access. A NumPy array in its simplest form is a Python object build around a C array. That is, it has a pointer to a contiguous data buffer of values. A Python list, on the other hand, has a pointer to a contiguous buffer of pointers, each of which points to a Python object which in turn has references to its data (in this case, integers). <img src='data/array_vs_list.png' style="width: 500px;">

Why Numpy is so fast?

Computations follow the Single Instruction Multiple Data (SIMD) paradigm. So that NumPy can take advantage of vectorized instructions on modern CPUs, like Intel's SSE and AVX, AMD's XOP.
<img src='data/devectorized.png' style="width: 350px;", caption='asdf'>
<img src='data/vectorized.png' style="width: 350px;">
A NumPy array is described by metadata (number of dimensions, shape, data type, strides, and so on) and the data (which is stored in a homogeneous and contiguous blocks of memory).
Array computations can be written very efficiently in a low-level language like C (and a large part of NumPy is actually written in C). Aditionally many internal methods and functions are linked to highly optimized linear algebra libraries like BLAS (Basic Linear Algebra Subprograms) and LAPACK (Linear Algebra PACKage).
Spatial locality in memory access patterns results in significant performance gains, notably thanks to the CPU cache.  Indeed, the cache loads bytes in chunks from RAM to the CPU registers.

Lets see the problem...
Consider the problem of calculating the squared norm of a vector ($\displaystyle || \mathbf{v} ||^2 = \mathbf{v} \cdot \mathbf{v}$) with the following 4 implementations:
End of explanation
"""


def id(x):
    """
    This function returns the memory
    block address of an array.
    """
    return x.__array_interface__['data'][0]

"""
Explanation: As you can see, some interesting things have happened:

The naive NumPy implementation, which iterates over data is actually slower than simply using a list. This is because the array stores a  very low-level representation of the numbers it stores , and this must be converted into Python-compatible objects before being returned to user, causing this extra overhead each time you index an array.
norm_square_numpy is slower than the clever NumPy implementation because two reasons:
There is time spend in allocating memory for storing the temporary result (vector x vector) and
This creates two implied loops, one to do the multiplication and one to do the sum.


The clever implementation uses np.dot() NumPy function which has no need to store intermediate results, and iterates just one time (but at C speed).

<div id='efficient' />
Eficient programming with NumPy
In-place and implicit copy operations
Prefer in-place over implicit-copy operations whenever possible. This will save memory (less work to garbage collector) and performs faster.
End of explanation
"""


a = np.zeros(10); aid = id(a)

# in-place operation
a *= 2; id(a) == aid

# implicit-copy operation
a = a * 2; id(a) == aid

"""
Explanation: Array computations can involve in-place operations (first example below: the array is modified) or implicit-copy operations (second example: a new array is created).
End of explanation
"""


%%timeit 
a = np.ones(100000000)
a *= 2

%%timeit
a = np.ones(100000000)
b = a * 2

"""
Explanation: Be sure to choose the type of operation you actually need. Implicit-copy operations are slower!
End of explanation
"""


m, n = 1000000, 100
a = np.random.random_sample((m, n))
index = np.arange(0, m, 10)

#fancy indexing - indexing with lists
%timeit a[index,:]

#memory slice - memory views
%timeit a[::10]

"""
Explanation: Efficient memory access
We have basically three alternatives to access arrays without loops:

array slicing
boolean masks
fancy indexing. 

Note. If you find yourself looping over indices to select items on which operation is performed, it can probably be done more efficiently with one of these techniques!
End of explanation
"""


def naive_indexing(vect):
    ret = list()
    for val in vect:
        if val < 0.5: ret.append(val)
    return np.array(ret)

#data to occupy and mask of booleans
vect = np.random.random_sample(1000000)
mask = vect < 0.5
mask

#naive indexing
%timeit naive_indexing(vect)

#mask indexing
%timeit vect[mask]

#'improved' mask indexing
%timeit np.compress(mask, vect)

"""
Explanation: Array slices are implemented as memory views, i.e, refer to the original data buffer of an array, but with different offsets, shapes and strides.
Array views should be used whenever possible, but one needs to be careful about the fact that views refer to the original data buffer.
Fancy indexing is several orders of magnitude slower as it involves copying a large array.

Another useful indexing technique is the mask of booleans. Lets suppose we want to get all the elements on array with value less than 0.5
End of explanation
"""


# array([0,1,2]) + 5
np.arange(3) + 5

# array([[1, 1 ,1], [1, 1, 1], [1, 1, 1]]) + array([0, 1, 2])
np.ones((3,3)) + np.arange(3)

# array([[0], [1], [2]]) + array([0, 1 ,2])
np.arange(3).reshape((3,1)) + np.arange(3)

"""
Explanation: Broadcasting
There is no need to always reshape arrays before operate on them. This useful feature implemented in NumPy arrays is called Broadcasting rules. In the visualization below, the extra memory indicated by the dotted boxes is never allocated, but it can be convenient to think about the operations as if it is.
<img src='data/broadcasting.png' style="width: 600px;">
How it works: The two arrays to be operated must match in at least one dimension. Then the array with less dimensions will logically extended to match the dimensions of the other
End of explanation
"""


def naive_orth(Q,v):
    m,n = Q.shape
    for j in range(n):
        v -= np.dot(Q[:,j],v)*Q[:,j]
    return v
    
def vectorized_orth(Q,v):
    proy = np.dot(Q.transpose(),v)
    # v -= (proy*Q).sum(axis=1)
    v -= np.dot(Q,proy)
    return v

# Let's generate a random unitary matrix
# Q unitary matrix, dimensions 100000 x1000
m,n = 10000,100
A = 10 * np.random.random((m,n))
Q,R = np.linalg.qr(A, mode='reduced')
del R

# v will be the starting vector for orthogonalization
v = np.random.random(m)
v1 = v.copy()
v2 = v.copy()

%timeit naive_orth(Q,v1)

%timeit vectorized_orth(Q,v2)

"""
Explanation: <div id='vectorized' />
Vectorized Algorithms with NumPy
Vectorization: The process of converting a scalar implementation, which processes a single pair of operands at a time, to a vector implementation, which processes one operation on multiple pairs of operands at once.
Example: Gram-Shmidt Orthogonalization
The problem: Given a matrix $Q_{m\times n}$ with ($m>n$), find a vector $v$ orthogonal to the column space of $Q$. 
<img src='data/orthogonalization.jpg' style="width: 400px;">
End of explanation
"""


A = np.random.random((100,100))
b = np.random.random(100)
c = np.random.random(100)

### Matrix power ###
np.linalg.matrix_power(A,3)

### Cholesky decomposition ###
#np.linalg.cholesky(A) #A must be positive definite

### QR decomposition ###
np.linalg.qr(A, mode='reduced')

### SVD decomposition ###
np.linalg.svd(A, full_matrices=False)

### Eigenvectors ###
np.linalg.eig(A)

### Eigevalues ###
np.linalg.eigvals(A)

### Matrix or vector norm ###
np.linalg.norm(A, ord='fro')

### Condition number ###
np.linalg.cond(A, p=np.inf)

### Determinant ###
np.linalg.det(A)

### Linear solver Ax=b ###
np.linalg.solve(A,b)

### Least Squares Ax=b (over-determined) ###
np.linalg.lstsq(A,b)

### Inverse ###
np.linalg.inv(A)

### Pseudo-Inverse ###
np.linalg.pinv(A)

### and many more...
del A,b,c

"""
Explanation: <div id='bonus' />
Bonus: Useful libraries based in NumPy
numpy.linalg (Numpy's submodule)
End of explanation
"""


