import numpy as np
import matplotlib.pyplot as plt

def load(filename):
    f = open(filename, 'rb')
    return np.loadtxt(f, delimiter=',')

data = load('ex1data2.txt')
X = data[:,0:2]
y = data[:,2:3]
m = len(y)

for i in range(10):
    print('x = [%.0f %.0f], y = %.0f \n' % (X[i,0], X[i,1], y[i, :]))


def featureNormalize(X):
    X_norm = X
    if (len(X) <= 1):
        return np.zeros(X.shape), np.zeros((1, X.shape[1])), np.zeros((1, X.shape[1]))

    mu = np.mean(X, axis=0, keepdims=True)
    # https://stackoverflow.com/questions/27600207/why-does-numpy-std-give-a-different-result-to-matlab-std
    sigma = np.std(X, axis=0, keepdims=True, ddof=1)
    X_norm = (X - mu) / sigma
    return X_norm, mu, sigma

# print('featureNormalize() test')
# def magic(n):
#   n = int(n)
#   if n < 3:
#     raise ValueError("Size must be at least 3")
#   if n % 2 == 1:
#     p = np.arange(1, n+1)
#     return n*np.mod(p[:, None] + p - (n+3)//2, n) + np.mod(p[:, None] + 2*p-2, n) + 1
#   elif n % 4 == 0:
#     J = np.mod(np.arange(1, n+1), 4) // 2
#     K = J[:, None] == J
#     M = np.arange(1, n*n+1, n)[:, None] + np.arange(n)
#     M[K] = n*n + 1 - M[K]
#   else:
#     p = n//2
#     M = magic(p)
#     M = np.block([[M, M+2*p*p], [M+3*p*p, M+p*p]])
#     i = np.arange(p)
#     k = (n-2)//4
#     j = np.concatenate((np.arange(k), np.arange(n-k+1, n)))
#     M[np.ix_(np.concatenate((i, i+p)), j)] = M[np.ix_(np.concatenate((i+p, i)), j)]
#     M[np.ix_([k, k+p], [0, k])] = M[np.ix_([k+p, k], [0, k])]
#   return M

# X_test, mu_test, sigma_test = featureNormalize(np.array([[1], [2], [3]]))
# print(X_test, mu_test, sigma_test)

# X_test, mu_test, sigma_test = featureNormalize(magic(3))
# print(X_test, mu_test, sigma_test)

# X_test, mu_test, sigma_test = featureNormalize(np.concatenate((-np.ones((1, 3)), magic(3))))
# print(X_test, mu_test, sigma_test)
# print('test over')

X, mu, sigma = featureNormalize(X)

print(X)
print(mu)
print(sigma)

X = np.concatenate((np.ones((m, 1)), X), 1)
alpha = 0.1
num_iters = 400

def computeCostMulti(X, y, theta):
    m = len(y)
    error = X @ theta - y;
    return error.T @ error / (2 * m)

# print('computeCostMulti() test')
# X_test = np.array([[2, 1, 3], [7, 1, 9], [1, 8, 1], [3, 7, 4]])
# y_test = np.array([[2], [5], [5], [6]])
# theta_test = np.array([[0.4], [0.6], [0.8]])
# print(computeCostMulti(X_test, y_test, theta_test))
# print('test over')

def gradientDescentMulti(X, y, theta, alpha, num_iters):
    m = len(y)
    J_history = np.zeros((num_iters, 1))
    for iters in range(num_iters):
        theta = theta - alpha * (X.T @ (X @ theta - y)) / m
        J_history[iters] = computeCostMulti(X, y, theta)
    return theta, J_history

# print('gradientDescentMulti() test')
# X_test = np.array([[2, 1, 3], [7, 1, 9], [1, 8, 1], [3, 7, 4]])
# y_test = np.array([[2], [5], [5], [6]])
# print(gradientDescentMulti(X_test, y_test, np.array([[0.1], [-0.2], [0.3]]), 0.01, 10))
# print('test over')

theta = np.zeros((3, 1))
theta, _ =  gradientDescentMulti(X, y, theta, alpha, num_iters)

print('Theta computed from gradient descent: [%f, %f %f]' % (theta[0], theta[1], theta[2]))

x = (np.array([1650, 3]) - mu) / sigma
price = np.concatenate(([[1]], x), 1) @ theta
print('Predicted price of a 1650 sq-ft, 3 br house (using gradient descent): $%f' % price)

alpha = [1, 0.3, 0.1, 0.03, 0.01]
num_iters = range(50, 201)
theta = np.zeros((3, 1))
_, J_history = gradientDescentMulti(X, y, theta, alpha[2], num_iters[0])
plt.plot(range(1, num_iters[0]+1), J_history, '-b', linewidth = 2)
plt.xlabel('Number of iterations')
plt.ylabel('Cost J')
plt.show()

# replace the value of alhpa below best alhpa value you found above
alpha = 0.1
num_iters = 50

theta = np.zeros((3, 1))
theta, _ = gradientDescentMulti(X, y, theta, alpha, num_iters)

print('Theta computed from gradent descent: [%f, %f, %f]' % (theta[0], theta[1], theta[2]))
x = (np.array([1650, 3]) - mu) / sigma
price = np.concatenate(([[1]], x), 1) @ theta
print('Predicted price of a 1650 sq-ft, 3 br house (using gradient descent): %f\n' % price)


data = load('ex1data2.txt')
X = data[:, 0:2]
y = data[:, 2:3]
m = len(y)

ax = plt.axes(projection = '3d')
x0 = X[:, 0:1]
x1 = X[:, 1:2]
x0 = np.linspace(np.min(x0), np.max(x0), 10)
x1 = np.linspace(np.min(x1), np.max(x1), 10)
x0, x1 = np.meshgrid(x0, x1)
X_norm = (np.array([x0,x1]).T - mu) / sigma
X_norm = np.concatenate((np.ones((X_norm.shape[0], X_norm.shape[1], 1)), X_norm), 2)
# z = x0 * theta[1] + x1 * theta[2] + theta[0]
z = (X_norm @ theta).T
ax.plot_surface(x0, x1, z.reshape(x0.shape), color='w', alpha=0.5)
ax.scatter3D(X[:,0:1], X[:,1:2], y, c='b', marker='^')
ax.scatter3D(1650, 3, price, c='r', marker='o')
ax.invert_xaxis()

plt.xlabel('sq-ft')
plt.ylabel('br')
plt.show()


data = load('ex1data2.txt')
X = data[:, 0:2]
y = data[:, 2:3]
m = len(y)

X = np.concatenate((np.ones((m, 1)), X), 1)

def normalEqn(X, y):
    return np.linalg.pinv((X.T @ X)) @ X.T @ y

# print('normalEqn() test')
# X_test = np.array([[2, 1, 3], [7, 1, 9], [1, 8, 1], [3, 7, 4]])
# y_test = np.array([[2], [5], [5], [6]])
# print(normalEqn(X_test, y_test))
# print('test over')

theta = normalEqn(X, y)
print('Theta computed from the normal equations: [%f, %f, %f]' % (theta[0], theta[1], theta[2]))

price = np.array([[1, 1650, 3]]) @ theta
print('Predicted price of a 1650 sq-ft, 3 br house (using norma equations): $%f' % price)

# plot normalEqn
ax = plt.axes(projection = '3d')
x0 = X[:, 1:2]
x1 = X[:, 2:3]
x0 = np.linspace(np.min(x0), np.max(x0), 10)
x1 = np.linspace(np.min(x1), np.max(x1), 10)
x0, x1 = np.meshgrid(x0, x1)
z = x0 * theta[1] + x1 * theta[2] + theta[0]
ax.plot_surface(x0, x1, z, color='w', alpha=0.5)
ax.scatter3D(X[:,1:2], X[:,2:3], y, c='b', marker='^')
ax.scatter3D(1650, 3, price, c='r', marker='o')
ax.invert_xaxis()

plt.xlabel('sq-ft')
plt.ylabel('br')
plt.show()

