ex5

import os
import numpy as np
import scipy.io as sio
import matplotlib.pyplot as plt


# Read the dataset
path = os.path.dirname(os.path.realpath(__file__)) + '\ex5\ex5data1'
data = sio.loadmat(path)

# training set
X = data['X']
y = data['y']
# validation set
Xval = data['Xval']
yval = data['yval']
# test set
Xtest = data['Xtest']
ytest = data['ytest']

# Number of training examples
m = X.shape[0]
# Learning rate
learningRate = 0
# Theta
theta = np.ones(X.shape[1]+1)

# Visualize the data
def plotData():
    plt.figure(figsize=(10,8))
    plt.ylabel('Water flowing out of the dam (y)')
    plt.xlabel('Change in water level (x)')
    plt.plot(X, y, 'rx', label='Training')
    plt.legend()

plotData()

##############################################################
# 1.2. Regularized linear regression cost function
##############################################################
def computeCost(theta, X, y, learningRate):
    theta = np.matrix(theta).T
    X = np.insert(X, 0, 1, axis=1)
    y = np.matrix(y)

    inner = (1 / (2 * m)) * np.sum(np.power((X.dot(theta) - y), 2))
    reg = np.sum(np.power(theta[1:], 2))
    return inner + (learningRate / (2 * m)) * reg

print(computeCost(theta, X, y, learningRate))

##############################################################
# 1.3. Regularized linear regression gradient
##############################################################
def computeGradient(theta, X, y, learningRate):
   theta = np.matrix(theta).T
   X = np.insert(X, 0, 1, axis=1)
   y = np.matrix(y)

   parameters = int(theta.ravel().shape[1])
   grad = np.zeros(parameters)

   error = X.dot(theta) - y

   for i in range(parameters):
       term = np.multiply(error.T, X[:,i])

       if (i == 0):
          grad[i] = np.sum(term) / m
       else:
          grad[i] = (np.sum(term) / m) + ((learningRate / m) * theta[i])

   return grad

print(computeGradient(theta, X, y, learningRate))

##############################################################
# 1.4. Fitting linear regression
##############################################################
import scipy.optimize

def optimizeTheta(theta, X, y, learningRate=0.):
   fit_theta = scipy.optimize.fmin_cg(computeCost, x0=theta,
                                       fprime=computeGradient,
                                       args=(X,y,learningRate),
                                       epsilon=1.49e-12,
                                       maxiter=1000)
   return fit_theta

fit_theta = optimizeTheta(theta, X, y, 0.)
print(computeCost(fit_theta, X, y, 0))
# Plot the data
XOne = np.insert(X, 0, values=np.ones(m), axis=1)
plotData()
plt.plot(X, XOne.dot(fit_theta).flatten())

##############################################################
# 2.1. Learning curves
##############################################################
def plotLearningCurve():
    """
    Loop over first training point, then first 2 training points, then first 3 ...
    and use each training-set-subset to find trained parameters.
    With those parameters, compute the cost on that subset (Jtrain)
    remembering that for Jtrain, lambda = 0 (even if you are using regularization).
    Then, use the trained parameters to compute Jval on the entire validation set
    again forcing lambda = 0 even if using regularization.
    Store the computed errors, error_train and error_val and plot them.
    """
    mym, error_train, error_val = [], [], []
    for i in range(1,13,1):
       x_subset = X[:i]
       y_subset = y[:i]
       fit_theta = optimizeTheta(theta, x_subset, y_subset, learningRate=0)
       error_train.append(computeCost(fit_theta, x_subset, y_subset, learningRate=0.))
       error_val.append(computeCost(fit_theta, Xval, yval, learningRate=0))
       mym.append(i)

    plt.figure(figsize=(10,8))
    plt.ylabel('Error')
    plt.xlabel('Number of training examples')
    plt.title('Learning curve for linear regression')
    plt.plot(mym, error_train, 'b-', label='Train')
    plt.plot(mym, error_val, 'g-', label='Cross Validation')
    plt.legend()

plotLearningCurve()

##############################################################
# 3. Polynomial Regression
##############################################################
def genPolyFeatures(X, p):
   myX = X.copy()
   for i in range(2, p+1):
      myX = np.insert(myX, myX.shape[1], np.power(myX[:,0], i), axis=1)

   return myX

#print(genPolyFeatures(X, 2))
def featureNormalize(X):
    Xnorm = X.copy()
    stored_feature_means = np.mean(Xnorm,axis=0) #column-by-column
    Xnorm[:,:] = Xnorm[:,:] - stored_feature_means[:]
    stored_feature_stds = np.std(Xnorm,axis=0,ddof=1)
    Xnorm[:,:] = Xnorm[:,:] / stored_feature_stds[:]
    return Xnorm, stored_feature_means, stored_feature_stds

##############################################################
# 3.1. Learning polynomial regression
##############################################################
degree = 6
newX = genPolyFeatures(X, degree)
newX_norm, stored_means, stored_stds = featureNormalize(newX)
mytheta = np.ones((newX_norm.shape[1]+1, 1))
fit_theta = optimizeTheta(mytheta, newX_norm, y, 1)
#print(len(fit_theta))

def plotFit(fit_theta,means,stds):
    """
    Function that takes in some learned fit values (on feature-normalized data)
    It sets x-points as a linspace, constructs an appropriate X matrix,
    un-does previous feature normalization, computes the hypothesis values,
    and plots on top of data
    """
    n_points_to_plot = 50
    xvals = np.linspace(-50, 50, n_points_to_plot)

    xmat = np.ones((n_points_to_plot,1)) # shape (50,1)
    xmat = np.insert(xmat,xmat.shape[1],xvals.T,axis=1) # shape (50,2)

    for i in range(2, len(fit_theta)):
        xmat = np.insert(xmat, xmat.shape[1], np.power(xmat[:,1], i), axis=1)
    #xmat_new = np.insert(xmat,0,1,axis=1)
    #print(xmat.shape)
    #This is undoing feature normalization
    xmat[:,1:] = xmat[:,1:] - means[:]
    xmat[:,1:] = xmat[:,1:] / stds[:]
    plotData()
    plt.plot(xvals, xmat.dot(fit_theta).flatten(), 'b--')

plotFit(fit_theta,stored_means,stored_stds)

##############################################################
# 3.2. Adjusting the regularization parameter
##############################################################

##############################################################
# 3.3. Selecting λ using a cross validation set
##############################################################

plt.show()