#-- Reading the (training) data in a data frame
df = pd.read_csv("pm25_train_data.csv",delimiter=";")
df.dropna(inplace=True)
df.head()

#-- Save the explanatory variables in a variable X (and their names in a varaible feature_names), and the target variable in Y
feature_names = list(df.columns[1:])
print("Feature Names: ",feature_names)

# extract target
Y = np.array(df["PM2.5"])

# extract variables
X = np.array(df[feature_names])

#-- Display some statistics on the data by using the describe function on the dataframe
df.describe()

# PM2.5 vs C0
df.plot.scatter("CO","PM2.5")

# plot scatter for each variable
for i in feature_names:
    df.plot.scatter("PM2.5",i,title="PM2.5 vs. "+str(i))

# correlation matrix
df.corr()

#-- Select the explanatory variables for the simple linear regression, then the multiple linear regression, and display the scatter plots
plt.figure(figsize=(10,5))

plt.subplot(131)
plt.scatter(df["PM2.5"],df["CO"])
plt.title("PM2.5 vs. CO")

plt.subplot(132)
plt.scatter(df["PM2.5"],df["NO2"],c="orange")
plt.title("PM2.5 vs. NO2")

plt.subplot(133)
plt.scatter(df["PM2.5"],df["SO2"],c="green")
plt.title("PM2.5 vs. SO2")

plt.show()

#-- Extract the data and creates two X matrices that will be used for the regression (have a look at page 26 to know the form of X):
#---- Xs for simple lin reg and Xm for multiple lin reg
from sklearn.model_selection import train_test_split
df = df[["PM2.5","CO","NO2","SO2"]]


X_train_single, X_test_single, y_train_single, y_test_single = train_test_split(df["CO"], df["PM2.5"], test_size=0.1, random_state=42)
X_train_multi, X_test_multi, y_train_multi, y_test_multi = train_test_split(df[["CO","NO2","SO2"]], df["PM2.5"], test_size=0.1, random_state=42)

# transfer to np arrays
X_train_single = np.array(X_train_single)
X_test_single = np.array(X_test_single)
y_train_multi = np.array(y_train_multi)
y_train_single = np.array(y_train_single)


X_train_multi = np.array(X_train_multi)
X_test_multi = np.array(X_test_multi)
y_test_multi = np.array(y_test_multi)
y_train_multi = np.array(y_train_multi)





#-- Check the size of both matrices
print("X_train_single:\t",len(X_train_single),type(X_train_single))
print("X_test_single:\t",len(X_test_single),type(X_test_single))
print("X_train_multi:\t",len(X_train_multi),type(X_train_multi))
print("X_test_multi:\t",len(X_test_multi),type(X_test_multi))
## Hint: use stack/hstack/vstack

#--- Write the standardisation function to mean-center the X data
from sklearn.linear_model import LinearRegression
def standardisation(X):
    """
    from sklearn import preprocessing
    scaler = preprocessing.StandardScaler().fit(X)
    return(scaler.transform(X))
    """
    scaled = (X-X.mean(axis=0))/X.std(axis=0,ddof=1)
    return(scaled)




    
#-- Test 1 - simple lin reg
X_train_single = standardisation(X_train_single)
print("Single:\n",stats.describe(X_train_single))
#-- Test 2 - multiple lin reg
X_train_multi = standardisation(X_train_multi)
print("\nMulti:\n",stats.describe(X_train_multi))

#-- To compare the results of the gradient descent algorithm, we will first implement an exact solution with the maximum likelihood
#Formulae recall: (X^T X)^-1 X^T Y
def coef_ml(X,Y):
    beta = np.dot(np.dot(np.linalg.inv(np.dot(X.T,X)),X.T),Y)
    return(beta)
    
#-- Test 1  - simple reg
beta_single = coef_ml(X_train_single,y_train_single)
print("beta single: \t",beta_single)
#-- Test 2 - multiple reg

#-- (We can also use the native functions of Scikit-Learn, but they are more parameters that needs to be tuned)
beta_multi = coef_ml(X_train_multi,y_train_multi)
print("beta multiple: \t",coef_ml(X_train_multi,y_train_multi))

# output: n(input)+1

#-- Model
def f(X,beta):
    return(np.dot(X,beta)) 

#-- Test 1 - simple reg
Yhat1 = f(X_train_single,beta_single)
print(Yhat1[:5])

#-- Test 2 - multiple reg
Yhat3 = f(X_train_multi,beta_multi)
print(Yhat3[:5])

display(Image(filename='cost_function.png'))