#!/usr/bin/env python
# coding: utf-8

# # Numpy

# $\newcommand{\re}{\mathbb{R}}$
# Numpy provides many useful facilities to perform numerical computations including vectors, matrices and linear algebra.

# In[1]:


import numpy
print(numpy.pi)
print(numpy.sin(numpy.pi/2))


#  It is common practice to import numpy like this.

# In[2]:


import numpy as np


# Note that ```np``` acts as an alias or short-hand for ```numpy```.

# In[3]:


print(np.pi)
print(np.sin(np.pi/2))


# ## 1-d Arrays

# Create an array of zeros

# In[4]:


x = np.zeros(5)
print(x)


# These one dimensional arrays are of type ```ndarray```

# In[5]:


type(x)


# Create an array of ones

# In[6]:


x = np.ones(5)
print(x)


# Create and add two arrays

# In[7]:


x = np.array([1.0, 2.0, 3.0])
y = np.array([4.0, 5.0, 6.0])
z = x + y
print(z)


# Get the size of array

# In[8]:


print(len(x))
print(x.size)


# Get the shape of an array

# In[9]:


print(x.shape)


# We see that these are arrays of reals by default. We can specify the type

# In[10]:


a = np.zeros(5, dtype=int)
print(a)


# ## linspace

# This behaves same way as Matlab's [linspace](https://mathworks.com/help/matlab/ref/linspace.html) function.
# 
# Generate 10 uniformly spaced numbers in [1,10]

# In[11]:


x = np.linspace(1,10,10)
print(x)


# Note that this includes the end points 1 and 10. The output of linspace is an ```ndarray``` of floats.

# In[12]:


type(x)


# `x = linspace(a,b,n)` is such that `x` is an array of `n` elements
# ```
# x[i] = a + i*h,   i=0,1,2,...,n-1,    h = (b-a)/(n-1)
# ```
# so that
# ```
# x[0] = a,   x[-1] = b
# ```

# ## arange

# In[13]:


x = np.arange(1,10)
print(x)
print(type(x))


# For $m < n$, `arange(m,n)` returns
# $$
# m, m+1, \ldots, n-1
# $$
# Note the last value $n$ is not included. 
# 
# We can specify a step size

# In[14]:


x = np.arange(1,10,2)
print(x)


# If the arguments are float, the returned value is also float.

# In[15]:


x = np.arange(1.0,10.0)
print(x)


# In[16]:


x = np.arange(0,1,0.1)
print(x)


# ## Beware of pitfalls - 1
# 
# Create an array of ones.

# In[17]:


x = np.ones(10)
print(x)


# Maybe we want set all elements to zero, so we might try this

# In[18]:


x = 0.0
print(x)


# ```x``` has changed from an array to a scalar !!! The correct way is this.

# In[19]:


x = np.ones(10)
print(x)
x[:] = 0.0
print(x)


# ## Beware of pitfalls - 2

# In[20]:


x = np.ones(5)
y = x
x[:] = 0.0
print(x)
print(y)


# Why did ```y``` change ? This happened because when we do
# ```
# y = x
# ```
# then `y` is just a pointer to `x`, so that changing `x` changes `y` also. If we want ```y``` to be an independent copy of ```x``` then do this

# In[21]:


x = np.ones(5)
y = x.copy()     # or y = np.copy(x)
x[:] = 0.0
print(x)
print(y)


# ## 2-d Arrays

# 2-d arrays can be considered as matrices, though Numpy has a separate [matrix](https://numpy.org/doc/stable/reference/generated/numpy.matrix.html) class.
# 
# Create an array of zeros

# In[22]:


A = np.zeros((5,5))
print(A)


# We can access the elements by two indices

# In[23]:


A[1,1] = 1.0
print(A)


# Create an array of ones

# In[24]:


A = np.ones((2,3))
print(A)


# Create identity matrix

# In[25]:


A = np.eye(5)
print(A)


# Create an array by specifying its elements

# In[26]:


A = np.array([[1.0, 2.0], 
              [3.0, 4.0]])
print(A)


# Create a random array and inspect its shape

# In[27]:


m, n = 2, 3
A = np.random.rand(m,n)
print(A)
print('A.size =',A.size)
print('A.shape =',A.shape)
print('A.shape[0] =',A.shape[0],' (number of rows)')
print('A.shape[1] =',A.shape[1],' (number of columns)')


# To print the elements of an array, we need two for loops, one over rows and another over the columns.

# In[28]:


for i in range(m): # loop over rows
    for j in range(n): # loop over columns
        print(i,j,A[i,j])


# To transpose a 2-d array

# In[29]:


A = np.array([[1,2],[3,4]])
print("A ="); print(A)
B = A.T
print("B ="); print(B)


# ## Diagonal matrix creation
# 
# Create
# 
# $$
# A = \begin{bmatrix}
# 1 & 0 & 0 & 0 \\
# 0 & 2 & 0 & 0 \\
# 0 & 0 & 3 & 0 \\
# 0 & 0 & 0 & 4 \end{bmatrix}
# $$

# In[30]:


a = np.array([1,2,3,4])    # diagonal elements
A = np.diag(a)
print(A)


# Create tri-diagonal matrix

# In[31]:


a = np.array([1,2,3])    # sub-diagonal  : length = n-1
b = np.array([4,5,6,7])  # main diagonal : length = n
c = np.array([-1,-2,-3]) # super-diagonal: length = n-1
A = np.diag(a,-1) + np.diag(b,0) + np.diag(c,+1)
print(A)


# ## empty_like, etc.
# 
# Sometimes we have an existing array and we want to create a new array of the same shape and type, but without initializing its elements.

# In[32]:


A = np.random.rand(2,3)
B = np.empty_like(A)
print(A.shape,B.shape)


# In[33]:


C = np.zeros_like(A)
print(C)


# In[34]:


D = np.ones_like(A)
print(D)


# ## 1-D array is neither row nor column vector

# In[35]:


x = np.array([1,2,3])
print(x.shape, x)
y = x.T
print(y.shape, y)


# Create a row vector

# In[36]:


x = np.array([[1,2,3]]) # row vector
print('x.shape =',x.shape)
print(x)
y = x.T
print('y.shape =',y.shape)
print(y)


# Create a column vector

# In[37]:


x = np.array([[1],[2],[3]]) # column vector
print('x.shape =',x.shape)
print(x)
y = x.T
print('y.shape =',y.shape)
print(y)


# Functions like `zeros` and `ones` can return row/column vector rather than array.

# In[38]:


x = np.ones((3,1)) # column vector
print('x ='); print(x)
y = np.ones((1,3)) # row vector
print('y ='); print(y)


# Row/column vectors are like row/column matrix. We have to use two indices to access the elements of such row/column vectors.
# 
# Here we access elements of a column vector.

# In[39]:


print(x[0][0])
print(x[1][0])
print(x[2][0])


# `x[0]` gives an array for the first row.

# In[40]:


print('First row of x       =',x[0], type(x[0]))
print('Element in first row =',x[0][0], type(x[0][0]))


# ## Accessing portions of arrays
# Array of 10 elements
# 
# | x[0] | x[1] | x[2] | x[3] | x[4] | x[5] | x[6] | x[7] | x[8] | x[9] |
# |:----:|:----:|:----:|:----:|:----:|:----:|:----:|:----:|:----:|:----:|
# | 0    | 1    | 2    | 3    | 4    | 5    | 6    | 7    | 8    | 9    |
# |x[-10]| x[-9]| x[-8]| x[-7]| x[-6]| x[-5]| x[-4]| x[-3]| x[-2]|x[-1] |
# 

# In[41]:


x = np.linspace(0,9,10)
print(x)


# Get elements ```x[2],...,x[5]```

# In[42]:


print(x[2:6])


# Hence ```x[m:n]``` gives the elements ```x[m],x[m+1],...,x[n-1]```.
# 
# Get elements from ```x[5]``` upto the last element

# In[43]:


print(x[5:])


# Get the last element

# In[44]:


print(x[-1])


# Get element ```x[5]``` upto last but one element

# In[45]:


print(x[5:-1])


# Access every alternate element of array

# In[46]:


print(x[0::2])


# In[47]:


print(x[1::2])


# These operations work on multi dimensional arrays also.

# In[48]:


A = np.random.rand(3,4)
print(A)


# In[49]:


print(A[0,:]) # 0'th row


# In[50]:


print(A[:,0]) # 0'th column


# In[51]:


print(A[0:2,0:3]) # print submatrix


# In[52]:


A[0,:] = 0.0 # zero out zeroth row
print(A)


# ## Arithmetic operations on arrays

# Arithmetic operations act element-wise

# In[53]:


x = np.array([1.0, 2.0, 3.0])
y = np.array([4.0, 5.0, 6.0])
print(x*y)  # multiply


# In[54]:


print(x/y)  # divide


# In[55]:


print(y**x) # exponentiation


# In[56]:


A = np.ones((3,3))
print(A*x)


# If ```A``` and ```x``` are arrays, then ```A*x``` does not give matrix-vector product. For that use ```dot```

# In[57]:


print(A.dot(x))


# or equivalently

# In[58]:


print(np.dot(A,x))


# In Python3 versions, we can use `@` to achieve matrix operations

# In[59]:


print(A@x)


# We can of course do matrix-matrix products using ```dot``` or ```@```

# In[60]:


A = np.ones((3,3))
B = 2*A
print('A =\n',A)
print('B =\n',B)
print('A*B =\n',A@B)


# On vectors, `@` performs dot product, i.e., `x@y = dot(x,y)`

# In[61]:


x = np.array([1,1,1])
y = np.array([2,3,4])
print(x@y)


# ## Some matrix/vector functions

# In[62]:


x = np.array([-3,-2,-1,0,1,2,3])
print('min     = ',np.min(x))
print('max     = ',np.max(x))
print('abs min = ',np.abs(x).min()) # np.min(np.abs(x))
print('abs max = ',np.abs(x).max()) # np.max(np.abs(x))
print('sum     = ',np.sum(x))
print('dot     = ',np.dot(x,x)) 
print('dot     = ',np.sum(x*x))


# We can compute vector norms using [numpy.linalg.norm](https://numpy.org/doc/stable/reference/generated/numpy.linalg.norm.html) (also see [scipy.linalg.norm](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.norm.html))

# In[63]:


from numpy.linalg import norm
print(norm(x))        # L2 norm
print(norm(x,2))      # L2 norm
print(norm(x,1))      # L1 norm
print(norm(x,np.inf)) # Linf norm


# or import the `linalg` module and use methods inside it.

# In[64]:


import numpy.linalg as la
print(la.norm(x))        # L2 norm
print(la.norm(x,2))      # L2 norm
print(la.norm(x,1))      # L1 norm
print(la.norm(x,np.inf)) # Linf norm


# These methods work on 2-d arrays which are square.

# In[65]:


A = 2*np.random.rand(3,3) - 1
print(A)
print('min       = ',np.min(A))
print('max       = ',np.max(A))
print('abs min   = ',np.abs(A).min()) # np.min(np.abs(x))
print('abs max   = ',np.abs(A).max()) # np.max(np.abs(x))
print('sum       = ',np.sum(A))
print('Frob norm = ',la.norm(A))        # Frobenius norm
print('L1 norm   = ',la.norm(A,1))      # L1 norm
print('L2 norm   = ',la.norm(A,2))      # L2 norm
print('Linf norm = ',la.norm(A,np.inf)) # Linf norm


# ## Example: Matrix-vector product

# Let
# $$
# x, y \in \re^n, \qquad A \in \re^{n \times n}
# $$
# To compute the matrix-vector product $y=Ax$, we can do it element-wise
# $$
# y_i = \sum_{j=0}^{n-1} A_{ij} x_j, \qquad 0 \le i \le n-1
# $$

# In[66]:


n = 10
x = np.random.rand(n)
A = np.random.rand(n,n)
y = np.zeros(n)
for i in range(n):
    for j in range(n):
        y[i] += A[i,j]*x[j]


# We can verify that our result is correct by computing $\| y - A x\|$

# In[67]:


print(np.linalg.norm(y-A@x))


# We can also compute the product column-wise. Let
# $$
# A_{:,j} = \textrm{j'th column of A}
# $$
# Then the matrix-vector product can also be written as
# $$
# y = \sum_{j=0}^{n-1} A_{:,j} x_j
# $$
# **Warning**: This may have inefficient memory access since by default, numpy arrays have column-major ordering, see below.

# In[68]:


y[:] = 0.0
for j in range(n):
    y += A[:,j]*x[j]

# Now check the result
print(np.linalg.norm(y-A@x))


# ## Example: Matrix-Matrix product

# If $A \in \re^{m\times n}$ and $B \in \re^{n \times p}$ then $C = AB \in \re^{m \times p}$ is given by
# $$
# C_{ij} = \sum_{k=0}^{n-1} A_{ik} B_{kj}, \qquad 0 \le i \le m-1, \quad 0 \le j \le p-1
# $$

# In[69]:


m,n,p = 10,8,6
A = np.random.rand(m,n)
B = np.random.rand(n,p)
C = np.zeros((m,p))
for i in range(m):
    for j in range(p):
        for k in range(n):
            C[i,j] += A[i,k]*B[k,j]


# Let us verify the result is correct by computing the Frobenius norm

# In[70]:


print(np.linalg.norm(C - A@B))


# Another view-point is the following
# $$
# C_{ij} = (\textrm{i'th row of A}) \cdot (\textrm{j'th column of B})
# $$

# In[71]:


for i in range(m):
    for j in range(p):
        C[i,j] = A[i,:] @ B[:,j]

# Now check the result
print(np.linalg.norm(C - A@B))


# ## Outer product of two vectors
# 
# Given $m$-vector $a$ and $n$-vector $b$, both considered as column vectors, their outer product
# 
# $$
# A = a b^\top = \begin{bmatrix} 
# |     &  |    &        & | \\
# a b_0 & a b_1 & \ldots & a b_{n-1} \\
# |     &  |    &        & | 
# \end{bmatrix}
# $$
# 
# is an $m \times n$ matrix with elements
# 
# $$
# A_{ij} = a_j b_j, \qquad 0 \le i \le m-1, \quad 0 \le n \le n-1
# $$

# In[72]:


a = np.array([1,2,3])
b = np.array([1,2])
A = np.outer(a,b)
print(A.shape)
print(A)


# ## Math functions
# 
# Numpy provides standard functions like `sin`, `cos`, `log`, etc. which can act on arrays in an element-by-element manner. This is not the case for functions in `math` module, which can only take scalar arguments.

# In[73]:


x = np.linspace(0.0, 2.0*np.pi, 5)
y = np.sin(x)
print('x =',x)
print('y =',y)


# In[74]:


A = np.random.rand(5,5) # 5x5 matrix
Y = np.sin(A)
print('A =\n',A)
print('Y =\n',Y)


# ## Memory ordering in arrays*
# By default, the ordering is same as in C/C++, the inner-most index is the fastest running one. For example, if we have an array of size (2,3), they are stored in memory in this order
# ```
# a[0,0], a[0,1], a[0,2], a[1,0], a[1,1], a[1,2]
# ```
# i.e.,
# ```
# a[0,0] --> a[0,1] --> a[0,2]
#                          |
#    |----------<----------|
#    |
# a[1,0] --> a[1,1] --> a[1,2]
# ```

# In[75]:


a = np.array([[1,2,3], [4,5,6]])
print('Row contiguous =', a[0,:].data.contiguous)
print('Col contiguous =', a[:,0].data.contiguous)
a.flags


# To get fortran style ordering, where the outer-most index is the fastest running one, which corresponds to the following layout in memory
# ```
# a[0,0], a[1,0], a[0,1], a[1,1], a[0,2], a[1,2]
# ```
# i.e.,
# ```
# a[0,0]    -- a[0,1]    -- a[0,2]
#    |     |     |      |      |  
#    v     ^     v      ^      v
#    |     |     |      |      |
# a[1,0] --    a[1,1] --    a[1,2]
# ```
# create like this

# In[76]:


b = np.array([[1,2,3], [4,5,6]], order='F')
print('Row contiguous =', b[0,:].data.contiguous)
print('Col contiguous =', b[:,0].data.contiguous)
b.flags


# ## Tensor product array: meshgrid

# In[77]:


x = np.linspace(0,3,4)
y = np.linspace(0,2,3)
X, Y = np.meshgrid(x,y)
print('len(x) = ',len(x))
print('len(y) = ',len(y))
print('shape X= ',X.shape)
print('shape Y= ',Y.shape)
print('x = ', x)
print('y = ', y)
print('X = '); print(X)
print('Y = '); print(Y)


# If $x \in \re^m$ and $y \in \re^n$, then the output is arranged like this
# $$
# X[i,j] = x[j], \qquad
# Y[i,j] = y[i], \qquad 0 \le i \le n-1, \quad 0 \le j \le m-1
# $$
# If we want the following arrangement
# $$
# X[i,j] = x[i], \qquad
# Y[i,j] = y[j], \qquad 0 \le i \le m-1, \quad 0 \le j \le n-1
# $$
# we have to do the following

# In[78]:


X, Y = np.meshgrid(x,y,indexing='ij')
print('len(x) = ',len(x))
print('len(y) = ',len(y))
print('shape X= ',X.shape)
print('shape Y= ',Y.shape)


# The second form is useful when working with finite difference schemes on Cartesian grids, where we want to use i index running along x-axis and j index running along y-axis.

# ## Reshaping arrays

# In[79]:


A = np.array([[1,2,3],
              [4,5,6]])
print(A)


# In[80]:


B = np.reshape(A,2*3,order='C')
print(B)


# In[81]:


A1 = np.reshape(B,(2,3),order='C')
print(A1)


# In[82]:


C = np.reshape(A,2*3,order='F')
print(C)


# In[83]:


A2 = np.reshape(C,(2,3),order='F')
print(A2)


# ## Writing and reading files
# Write an array to file

# In[84]:


A = np.random.rand(3,3)
print(A)
np.savetxt('A.txt',A)


# Check the contents of the file in your terminal
# ```
# cat A.txt
# ```
# We can also print it from within the notebook

# In[85]:


get_ipython().system('cat A.txt')


# Write two 1-D arrays as columns into file

# In[86]:


x = np.array([1.0,2.0,3.0,4.0])
y = np.array([2.0,4.0,6.0,8.0])
np.savetxt('data.txt',np.column_stack([x,y]))
get_ipython().system('cat data.txt')


# We can control the number of decimals saved, and use scientific notation

# In[87]:


np.savetxt('data.txt',np.column_stack([x,y]),fmt='%8.4e')
get_ipython().system('cat data.txt')


# We can read an existing file like this

# In[88]:


d = np.loadtxt('data.txt')
print('Shape d =', d.shape)
x1 = d[:,0]
y1 = d[:,1]
print('x =',x1)
print('y =',y1)


# ## Measuring time
# 
# Loops are slow in Python and it is good to avoid them. The following example show how to time code execution.
# 
# Create some random matrices

# In[91]:


m = 1000
A = np.random.random((m,m))
B = np.random.random((m,m))


# First we add two matrices with an explicit for loop.

# In[96]:


get_ipython().run_cell_magic('timeit', '-r10 -n10', 'C = np.empty_like(A)\nfor i in range(m):\n    for j in range(m):\n        C[i,j] = A[i,j] + B[i,j]\n')


# Now we use the + operator.

# In[97]:


get_ipython().run_cell_magic('timeit', '-r10 -n10', 'C = A + B\n')


# The loop version is significantly slower.