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

# In[39]:


import numpy as np
import time


# First, a motivating example for why we might want a NUMerical PYthon (numpy) library.

# In[40]:


num_elements = 1000000
input_array = np.arange(num_elements)


# In[41]:


start_time = time.time()
return_array = [0] * len(input_array)
for key, val in enumerate(input_array):
    return_array[key] = val * val
print(time.time() - start_time)


# In[42]:


start_time = time.time()
return_array_vectorized = np.power(input_array, 2)
print(time.time() - start_time)


# <b>Array operations</b>

# Array construction

# In[43]:


np.zeros(5)


# In[44]:


np.identity(3)


# In[45]:


array_1 = np.array([[1,2],[3,4],[5,6]])  # An array with arbitrary elements.
array_1


# Array access

# In[46]:


array_1[1,1]  # Access an element


# In[47]:


array_1[:,1]  # Access a row


# Array broadcasting of addition and scalar multiplication

# In[48]:


array_2 = np.add(array_1, 1)  # array_1 + 1
array_2


# In[49]:


np.multiply(array_1, 2) # array_1 * 2


# In[50]:


array_1 / 2.0  # np.multiply(array_1, 1/2.0)


# Array elementwise addition

# In[51]:


array_1 + array_2  # np.add(array_1, array_2)


# Array elementwise multiplication

# In[52]:


np.multiply(array_1, array_1)  # array_1 * array_1


# Dot-product

# In[53]:


array_3 = np.array([[1, 2]])
array_4 = np.array([[3], [4]])
np.dot(array_3, array_4)  # Inner-product


# In[54]:


array_5 = np.dot(array_4, array_3)  # Outer-product
array_5


# Array info

# In[55]:


array_5.sum(axis=0)


# In[56]:


array_5.mean(axis=0)


# In[57]:


array_5.std(axis=0)


# In[58]:


array_5.max(axis=0)


# In[59]:


array_5.shape


# Transpose

# In[60]:


M1 = array_1.T
M1


# Matrix elementwise exponentiation

# In[61]:


np.power(M1, 2)  # M1 ^ 2?


# Matrix multiplication

# In[62]:


np.matmul(M1.T, M1)  # or np.dot(M1.T, M1)


# Matrix multiplication is not commutative

# In[63]:


print(np.matmul(M1, M1.T))
print(np.matmul(M1.T, M1))


# Matrix Inverse

# In[64]:


np.linalg.inv(np.array([[2, 0], [0, 2]]))


# Inverse fails if matrix is singular.

# In[65]:


np.linalg.inv(np.matmul(M1.T, M1))  # This should raise an exception.


# Evaluating linear systems Ax = b for b

# In[66]:


A1 = np.array([[1, 1], [0, 1]])
x1 = np.array([[2], [2]])
b1 = np.matmul(A1, x1)
b1


# Solving linear systems Ax = b for x

# In[67]:


x1_verify = np.linalg.solve(A1, b1)
all(x1 == x1_verify)


# If A is singular then the linear system may be overdetermined (no solution) when solving for x.

# In[68]:


A2 = np.array([[1, 1], [2, 2]])  # Note that the rank of this matrix is 1.
x2 = np.array([[2], [3]])
potential_sol = np.matmul(A2, x2)
potential_sol


# In[69]:


b2 = np.array([[5], [11]])
x2_verify = np.linalg.solve(A2, b2)  # This should raise an exception.
# There is no selection of x to solve this linear system for A2 and b2.


# In[70]:


x2_verify = np.linalg.solve(A2, potential_sol)  # This should raise an exception.
# There happens to exist an x to solve this system, but numpy still fails.


# Another linear system that is underdetermined (many solutions) when solving for x

# In[71]:


A3 = np.array([[1, 0, 0], [0, 1, 1]])
x3 = np.array([[1], [1], [1]])
b3 = np.matmul(A3, x3)
b3


# In[72]:


np.linalg.solve(A3, b3)  # This should raise an exception.