import numpy as np
import matplotlib.pyplot as plt

np.random.seed(666)
x = 2 * np.random.random(size=100)
y = x*3. + 4. + np.random.normal(size=100)

X = x.reshape(-1,1)

#plt.scatter(x,y)
#plt.show()

def J(theta,X_b,y):
    try:
        return np.sum((y - X_b.dot(theta))**2)/len(X_b)
    except:
        return float('inf')

def dJ(theta,X_b,y):
    res = np.empty(len(theta))
    res[0] = np.sum(X_b.dot(theta) - y)
    for i in range(1,len(theta)):
        res[i] = (X_b.dot(theta) - y).dot(X_b[:,i])
    return res * 2 / len(X_b)

def  gradient_descent(X_b,y , initial_theta,eta,n_iters=1e4,epsilon=1e-8):
    theta = initial_theta
    i_iters = 0
    #theta_history.append(theta)

    while i_iters < n_iters:
        gradient = dJ(theta,X_b,y)
        last_theta = theta
        theta = theta - eta * gradient
        #theta_history.append(theta)
        if (abs(J(theta,X_b,y) - J(last_theta,X_b,y)) < epsilon):
            break
        i_iters +=1
    return theta

X_b = np.hstack([np.ones((len(X),1)), X])

initial_theta = np.zeros(X_b.shape[1])
print(initial_theta)
eta = 0.01
theta = gradient_descent(X_b,y,initial_theta,eta)
print(theta)