'''
numpy : is a mathematical computation library. 
matplotlib : A plotting library 
'''
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# First, size of the training set we want to generate.
observations = 2000

# We generate them randomly, drawing from an uniform distribution. There are 3 arguments of this method (low, high, size).
# The size of x0 and z0 is observations by 1. In this case: 2000 x 1.
x0 = np.random.uniform(low=-10, high=10, size=(observations,1))
z0 = np.random.uniform(low=-10, high=10, size =(observations,1))

# Combine  two dimensions of  x0 , z0 into one input matrix. 
# This is the X matrix from the linear model y = x*w + b.
# column_stack is a Numpy method, which combines two vectors into a matrix.
inputs = np.column_stack((x0,z0))

#display the dimensions of input matrix
inputs.shape

# noise is something which isn't in our hands There are 3 arguments of this method (low, high, size).
noise = np.random.uniform(-1, 1, (observations,1))

# Produce the targets according to the f(x,z) = 3x - 4z + 5 + noise definition.
# In this way, we are basically saying: the weights should be 3 and -4, while the bias is 5.
targets = 3*x0 - 4*z0 + 5 + noise

#display the dimensions of input matrix
targets.shape

'''
To use the 3D plot, the objects should have a certain shape, so we reshape the targets.
The proper method to use is reshape and takes as arguments the dimensions in which we want to fit the object.
'''
targets = targets.reshape(observations,)

# Plotting according to the conventional matplotlib.pyplot syntax
# Declare the figure
fig = plt.figure()

# A method allowing us to create the 3D plot
ax = fig.add_subplot(111, projection='3d')

# Choose the axes.
ax.plot(x0, z0, targets)

# Set labels
ax.set_xlabel('x0')
ax.set_ylabel('z0')
ax.set_zlabel('Our Targets')

# You can fiddle with the azim parameter to plot the data from different angles. Just change the value of azim=100
# to azim = 0 ; azim = 200, or whatever. Check and see what happens.
ax.view_init(azim=100)

# So far we were just describing the plot. This method actually shows the plot. 
plt.show()

''' 
We reshape the targets back to the shape that they were in before plotting.
This reshaping is a side-effect of the 3D plot. Sorry for that.
'''
targets = targets.reshape(observations,1)

'''
We will initialize the weights and biases randomly in some small initial range.
init_range is the variable that will measure that.
You can play around with the initial range, but we don't really encourage you to do so.
High initial ranges may prevent the machine learning algorithm from learning.
'''
init_range = 0.1

# Weights are of size k x m, where k is the number of input variables and m is the number of output variables
# In our case, the weights matrix is 2x1 since there are 2 inputs (x and z) and one output (y)
weights = np.random.uniform(low=-init_range, high=init_range, size=(2, 1))

# Biases are of size 1 since there is only 1 output. The bias is a scalar.
biases = np.random.uniform(low=-init_range, high=init_range, size=1)

#Print the weights to get a sense of how they were initialized.
print (weights)
print (biases)

'''
Set some small learning rate (denoted eta in the lecture). 
0.02 is going to work quite well for our example. Once again, you can play around with it.
it is HIGHLY recommended that you play around with it.
'''
learning_rate = 0.02

'''
We iterate over our training dataset 100 times. That works well with a learning rate of 0.02.
The proper number of iterations is something we will talk about later on, but generally
a lower learning rate would need more iterations, while a higher learning rate would need less iterations
keep in mind that a high learning rate may cause the loss to diverge to infinity, instead of converge to 0.
'''
for i in range (100):
    
    # This is the linear model: y = xw + b equation
    outputs = np.dot(inputs,weights) + biases
    # The deltas are the differences between the outputs and the targets
    # Note that deltas here is a vector 1000 x 1
    deltas = outputs - targets
        
    # We are considering the L2-norm loss, but divided by 2, so it is consistent with the lectures.
    # Moreover, we further divide it by the number of observations.
    # This is simple rescaling by a constant. We explained that this doesn't change the optimization logic,
    # as any function holding the basic property of being lower for better results, and higher for worse results
    # can be a loss function.
    loss = np.sum(deltas ** 2) / 2 / observations
    
    # We print the loss function value at each step so we can observe whether it is decreasing as desired.
    print ("Loss :" , loss)
    
    # Another small trick is to scale the deltas the same way as the loss function
    # In this way our learning rate is independent of the number of samples (observations).
    # Again, this doesn't change anything in principle, it simply makes it easier to pick a single learning rate
    # that can remain the same if we change the number of training samples (observations).
    # You can try solving the problem without rescaling to see how that works for you.
    deltas_scaled = deltas / observations
    
    # Finally, we must apply the gradient descent update rules from the relevant lecture.
    # The weights are 2x1, learning rate is 1x1 (scalar), inputs are 1000x2, and deltas_scaled are 1000x1
    # We must transpose the inputs so that we get an allowed operation.
    weights = weights - learning_rate * np.dot(inputs.T,deltas_scaled)
    biases = biases - learning_rate * np.sum(deltas_scaled)
    
    # The weights are updated in a linear algebraic way (a matrix minus another matrix)
    # The biases, however, are just a single number here, so we must transform the deltas into a scalar.
    # The two lines are both consistent with the gradient descent methodology. 

'''
We print the weights and the biases, so we can see if they have converged to what we wanted.
When declared the targets, following the f(x,z), we knew the weights should be 3 and -4, while the bias: 5.
'''
print (weights, biases)

'''
We print the outputs and the targets in order to see if they have a linear relationship.
Again, that's not needed. Moreover, in later lectures, that would not even be possible.
'''
plt.plot(outputs,targets)
plt.xlabel('outputs')
plt.ylabel('targets')
plt.show()