# -*- coding: utf-8 -*-
"""***CAP5625_Assigment 3-libaries_LogisticRidgeRegresion.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1jsQiNbecYaXL0vqUE09P5EQeZx1EyEnx

- **Programmer:**
    - **Shaun Pritchard**
    - **Ismael A Lopez** 

## **Assigment 3**
**Brief overview of assignment**

***perform a penalized
(regularized) logistic (multinomial) regression fit using ridge regression, with the model
parameters obtained by batch gradient descent. Your predictions will be based on K=5
continental ancestries (African, European, East Asian, Oceanian, or Native American). Ridge
regression will permit you to provide parameter shrinkage (tuning parameter 𝜆=0) to mitigate
overfitting. The tuning parameter 𝜆 will be chosen using five-fold cross validation, and the best-
fit model parameters will be inferred on the training dataset conditional on an optimal tuning
parameter. This trained model will be used to make predictions on new test data points***

> ***Table of Contents***
<hr>


* [Import Packages](#Import_Packages)
    * [Import packages for manipulating data](#Import_packages_for_manipulating_data)
    * [Import packages for splitting data](#Import_packages_for_splitting_data)
    * [Import packages for modeling data](#Import_packages_for_modeling_data)
    * [Import packages for Scaling and Centering data](#Import_packages_for_Scaling_and_Centering_data)
    * [Import packages for Measuring Model Perormance](#Import_packages_for_Measuring_Model_Perormance)
    
* [Data Processing](#Data_Processing)
    * [Import Data](#Import_data)
    * [Lets change the categorical values](#Lets_change_the_categorical_values)
    * [Create Predictor and Target numpy array](#Create_Predictor_and_Target_numpy_array)
    * [Create a Normalize copy of variables](#Create_a_Normalize_copy_of_variables)
    * [Split Data](#Split_Data:)
* [Regression Model](#Regression_Model)
    * [Define our learning rates:](#Define_our_learning_rates)
    * [Create the Regression Objects](#Create_the_Regression_Objects)
        * [LogisticRegression Library](#LogisticRegression_Library)

> ***Deliverables***
<hr>

* [**Deliverable 6.1**](#Deliverable_6.1)
* [**Deliverable 6.2**](#Deliverable_6.2)
* [**Deliverable 6.3**](#Deliverable_6.3)
* [**Deliverable 6.4**](#Deliverable_6.4)
* [**Deliverable 6. Reason for difference**](#Deliverable_6_Reason_for_difference)

# Import Packages <a class="anchor" id="Import_Packages"></a>

### Import packages for manipulating data  <a class="anchor" id="Import_packages_for_manipulating_data"></a>
"""

# Commented out IPython magic to ensure Python compatibility.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.mlab as mlab
import math
import csv
import random
# %matplotlib inline
from sklearn.preprocessing import LabelEncoder
import seaborn as sns
import math

"""### Import packages for splitting data  <a class="anchor" id="Import_packages_for_splitting_data"></a>"""

from sklearn.model_selection import train_test_split, cross_val_score, KFold, StratifiedKFold
from sklearn.model_selection import GridSearchCV

"""### Import packages for modeling data  <a class="anchor" id="Import_packages_for_modeling_data"></a>"""

# Import models:
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression as linearR_Model, Ridge  as RidgeR_Model
from sklearn.linear_model import RidgeCV

from sklearn.linear_model import ElasticNet
from sklearn.linear_model import ElasticNetCV

from sklearn.linear_model import LogisticRegression


from sklearn.exceptions import ConvergenceWarning
#from sklearn.utils._testing import ignore_warnings
import warnings
warnings.filterwarnings('ignore', category=ConvergenceWarning) # To filter out the Convergence warning
warnings.filterwarnings('ignore', category=UserWarning)
from itertools import product

"""### Import packages for Scaling and Centering data  <a class="anchor" id="Import_packages_for_Scaling_and_Centering_data"></a>"""

from sklearn.preprocessing import StandardScaler

"""### Import packages for Measuring Model Perormance  <a class="anchor" id="Import_packages_for_Measuring_Model_Perormance"></a>"""

from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score
from sklearn.metrics import make_scorer

"""# Data Processing <a class="anchor" id="Data_Processing"></a>

### Import Data  <a class="anchor" id="Import_data"></a>

***Traing Dataset***
"""

Train_dataset = pd.read_csv ('TrainingData_N183_p10.csv')
Train_dataset.head(3)

# What are the datatypes of each observation:
print(Train_dataset.dtypes)
# Shape of my data
print('The size of our data are: ',Train_dataset.shape)

print('Training Dataset Missing Values: \n',Train_dataset.isnull().sum())

"""***Test Dataset***"""

Test_dataset = pd.read_csv ('TestData_N111_p10.csv')
Test_dataset.head(3)

# What are the datatypes of each observation:
print(Test_dataset.dtypes)
# Shape of my data
print('The size of our data are: ',Test_dataset.shape)

# Are there any null or missing values
print('Test Dataset Missing Values: \n',Test_dataset.isnull().sum())

"""### Lets change the categorical values  <a class="anchor" id="Lets_change_the_categorical_values"></a>

"""

# recode the categories
Training_Class = Train_dataset['Ancestry'].unique().tolist()
Test_Class = Test_dataset['Ancestry'].unique().tolist()
num_features = len(Training_Class)


print("Unique Values for Train Ancestry: ", Training_Class)
print("Unique Values for Test Ancestry: ", Test_Class)

from sklearn.preprocessing import LabelEncoder
le = LabelEncoder()
Train_dataset['Ancestry_Encoded'] = le.fit_transform(Train_dataset.iloc[:,-1:])
Test_dataset['Ancestry_Encoded'] = le.fit_transform(Test_dataset.iloc[:,-1:])

Train_dataset.head(3)

"""### Create Predictor and Target numpy array  <a class="anchor" id="Create_Predictor_and_Target_numpy_array"></a>"""

# Target:
Y_Train= Train_dataset['Ancestry_Encoded'].to_numpy()
Y_Test= Test_dataset['Ancestry_Encoded'].to_numpy()
Y_Train.shape

# Convert the Pandas dataframe to numpy ndarray for computational improvement
X_Train = Train_dataset.iloc[:,:-2].to_numpy()
X_Test = Test_dataset.iloc[:,:-2].to_numpy()

print(type(X_Train), X_Train[:1], "Shape = ", X_Train.shape)

"""### Create a Normalize copy of variables <a class="anchor" id="Create_a_Normalize_copy_of_variables"></a>"""

# Create Standarizing ObjectPackages:
standardization = StandardScaler()

# Strandardize 
n_observations = len(Train_dataset)
variables = Train_dataset.columns


# Standardize the Predictors (X)
X_Train = standardization.fit_transform(X_Train)

# Add a constanct to the predictor matrix
#X_Train = np.column_stack((np.ones(n_observations),X_Train))


# Save the original M and Std of the original data. Used for unstandardize
original_means = standardization.mean_

# we chanced standardization.std_ to standardization.var_**.5
originanal_stds = standardization.var_**.5


print("observations :", n_observations)
print("variables :", variables[:2])
print('original_means :', original_means)
print('originanal_stds :', originanal_stds)

"""### Split Data: <a class="anchor" id="Split_Data:"></a>

#let's first split it into train and test part
X_train, X_out_sample, y_train, y_out_sample = train_test_split(Xst, y_Centered, test_size=0.40, random_state=101) # Training and testing split

X_validation, X_test, y_validation, y_test = train_test_split(X_out_sample, y_out_sample, test_size=0.50, random_state=101) # Validation and test split

# Print Data size
print ("Train dataset sample size: {}".format(len(X_train)))
print ("Validation dataset sample size: {}".format(len(X_validation)))
print ("Test dataset sample size: {}".format( len(X_test)))

# Regression Model <a class="anchor" id="Regression_Model"></a>
<hr>

### Define our learning rates <a class="anchor" id="Define_our_learning_rates"></a>
"""

# Define my tuning parameter values 𝜆:

learning_rates_λ = [1.e-04, 1.e-03, 1.e-02, 1.e-01, 1.e+00, 1.e+01, 1.e+02, 1.e+03, 1.e+04]
print(learning_rates_λ)

# learning rate
α =  1e-4

# K-folds
k = 5


# Itterations
n_iters = 10000

"""### Create the Regression Objects <a class="anchor" id="Create_the_Regression_Objects"></a>

**LogisticRegression Library** <a class="anchor" id="LogisticRegression_Library"></a>
"""

# LogisticRegression
from sklearn.linear_model import LogisticRegression
Library_LogisticRegression = LogisticRegression(max_iter = 10000, multi_class='multinomial', solver='lbfgs', penalty='l2', C=1)

"""## **Deliverable 7.1**  <a class="anchor" id="Deliverable_6.1"></a>
<h>

> Deliverable 1: Illustrate the effect of the tuning parameter on the inferred ridge regression coefficients by generating five plots (one for each of the 𝐾=5 ancestry classes) of 10 lines (one for each of the 𝑝=10 features), with the 𝑦-axis as 𝛽̂
𝑗𝑘, 𝑗=1,2,…,10 for the graph of class 𝑘, and 𝑥-axis the corresponding log-scaled tuning parameter value log10(𝜆) that
7
generated the particular 𝛽̂
𝑗𝑘. Label both axes in all five plots. Without the log scaling of the tuning parameter, the plot will look distorted.

**LogisticRegression with Library**
"""

L𝛽_per_λ=[] # set empty list

# Evaluate tuning parameters with LogisticRegression penalty
for tuning_param in learning_rates_λ:
        Library_LogisticRegression = LogisticRegression(max_iter = 10000, multi_class = 'multinomial', solver = 'lbfgs', penalty ='l2', C = tuning_param)
        Library_LogisticRegression.fit(X_Train, Y_Train)
        c = np.array(Library_LogisticRegression.coef_)
       # c = np.append(tuning_param,c)
        L𝛽_per_λ.append(Library_LogisticRegression.coef_)
#        print(c)

L𝛽_per_λ[0]

# Loop throught the betas, by class generated by each lamda
temp_df = []
for l in range(np.array(L𝛽_per_λ).shape[0]):
    for c in range(np.array(L𝛽_per_λ).shape[1]):
        temp_df.append(np.append(L𝛽_per_λ[l][c],(learning_rates_λ[l],c)))

TunnedL𝛽_df=pd.DataFrame(np.array(temp_df))
TunnedL𝛽_df.columns=['PC1', 'PC2', 'PC3', 'PC4', 'PC5', 'PC6', 'PC7', 'PC8', 'PC9', 'PC10', 'Lambda', 'Class']
#TunnedL𝛽_df['Class_Name'] = TunnedL𝛽_df['Class_Name'].apply(lambda x: Training_Class[int(x)])
TunnedL𝛽_df.head(10)

Training_Class

TunnedL𝛽_df[TunnedL𝛽_df.Class.eq(0)]

# Plot tuning parameter on the inferred ridge regression coefficients
sns.set(rc = {'figure.figsize':(15,8)})
for i, c in enumerate(Training_Class):
    sns.set_theme(style="whitegrid")
    sns.set_palette("mako")
    for j in range(1, 1 + X_Train.shape[1]):
        sns.lineplot( x =  TunnedL𝛽_df[TunnedL𝛽_df.Class.eq(i)]['Lambda'], y = TunnedL𝛽_df[TunnedL𝛽_df.Class.eq(i)]['PC{}'.format(j)], palette='mako',   label = 'PC{}'.format(j) )
        sns.set()
    plt.xscale('log')
    plt.legend(bbox_to_anchor=(1.09, 1), loc='upper left')
    plt.xlabel('Log Lambda')
    plt.ylabel('Coefficient Values')
    plt.suptitle('Inferred Ridge Regression Coefficient Tuning Parameters of' + ' ' + c + ' ' + 'Class')
    plt.show()

"""# **Deliverable 7.2**  <a class="anchor" id="Deliverable_6.2"></a>
Illustrate the effect of the tuning parameter on the cross validation error by generating a plot with the 𝑦-axis as CV(5) error, and the 𝑥-axis the corresponding log-scaled tuning parameter value log10(𝜆) that generated the particular CV(5) error. Label both axes in the plot. Without the log scaling of the tuning parameter 𝜆, the plots will look distorted.

**CV Elastic Net with Library**
"""

from sklearn.model_selection import GridSearchCV
from sklearn.linear_model import LogisticRegression

#Define the model
Library_LogisticRegression = LogisticRegression(max_iter = 10000, multi_class = 'multinomial', solver = 'lbfgs', penalty ='l2')


# Create the Kfold:
cv_iterator = KFold(n_splits = 5, shuffle=True, random_state=101)

cv_score = cross_val_score(Library_LogisticRegression, X_Train, Y_Train, cv=cv_iterator, scoring='neg_mean_squared_error', n_jobs=1)
print (cv_score)
print ('Cv score: mean %0.3f std %0.3f' % (np.mean(cv_score), np.std(cv_score)))

# define grid
Parm_grid = dict()
Parm_grid['C'] = learning_rates_λ
Parm_grid

# Lets define search
GsearchCV = GridSearchCV(estimator = Library_LogisticRegression, param_grid = Parm_grid, scoring = 'neg_mean_squared_error', n_jobs=1, refit=True, cv=cv_iterator)
GsearchCV.fit(X_Train, Y_Train)

GCV_df = pd.concat([pd.DataFrame(GsearchCV.cv_results_["params"]),pd.DataFrame(GsearchCV.cv_results_["mean_test_score"], columns=["mean_test_score"])],axis=1)
#GCV_df.index=GCV_df['alpha']
GCV_df.rename(columns={"C": "learning_rates_λ"}, inplace=True)

GCV_df

sns.set_theme(style="whitegrid")
sns.set_palette("mako")


plt.plot(GCV_df["learning_rates_λ"] , GCV_df["mean_test_score"])


sns.set_palette("mako")
sns.set()


plt.suptitle('Effect of the uning parameter on the cross validation error log10(lambda)')
plt.xscale('log')
plt.xlabel('Lambda')
plt.ylabel('loss from cross valiadation')
plt.show()

"""# **Deliverable 7.3**  <a class="anchor" id="Deliverable_6.3"></a>
Indicate the value of 𝜆 that generated the smallest CV(5) error

**Smallest CV with Library**
"""

print ('Best: ',GsearchCV.best_params_)
print ('Best CV mean squared error: %0.3f' % np.abs(GsearchCV.best_score_))

GCV_df.sort_values(by=['mean_test_score'], ascending=False)[:1]

# Alternative: sklearn.linear_model.ElasticNetCV
from sklearn.linear_model import LogisticRegressionCV

auto_LR = LogisticRegressionCV(Cs = learning_rates_λ, cv=5, max_iter = 10000, multi_class = 'multinomial', solver = 'lbfgs', penalty ='l2', n_jobs= 1 )
auto_LR.fit(X_Train, Y_Train)
#print ('Best alpha: %0.5f' % auto_LR.alpha_)
print ('Best 𝜆: ' , auto_LR.C_)

"""# **Deliverable 7.4**  <a class="anchor" id="Deliverable_6.4"></a>

Given the optimal 𝜆, retrain your model on the entire dataset of 𝑁=183 observations to obtain an estimate of the (𝑝+1)×𝐾 model parameter matrix as 𝐁̂ and make predictions of the probability for each of the 𝐾=5 classes for the 111 test individuals located in TestData_N111_p10.csv. That is, for class 𝑘, compute 
𝑝𝑘(𝑋;𝐁̂)=exp(𝛽̂0𝑘+Σ𝑋𝑗𝛽̂ 𝑗𝑘𝑝𝑗=1)
/ Σexp(𝛽̂0ℓ+Σ𝑋𝑗𝛽̂ 𝑗ℓ𝑝𝑗=1)

- for each of the 111 test samples 𝑋, and also predict the most probable ancestry label as 
    - 𝑌̂(𝑋)=arg max𝑘∈{1,2,…,𝐾}𝑝𝑘(𝑋;𝐁̂)
- Report all six values (probability for each of the 𝐾=5 classes and the most probable ancestry label) for all 111 test individuals.

**Tunned with best λ with Library**
"""

Library_LogisticRegression_best= LogisticRegression(max_iter = 10000, multi_class='multinomial', solver='lbfgs', penalty='l2', C= auto_LR.C_[0])
Library_LogisticRegression_best.fit( X_Train, Y_Train )

y_predM_best = Library_LogisticRegression_best.predict(X_Test)
print ("Betas= ", np.mean(Library_LogisticRegression_best.coef_, 0))

yhat = Library_LogisticRegression_best.predict_proba(X_Test)
# summarize the predicted probabilities
print('Predicted Probabilities: %s' % yhat[0])

ŷ_test = Library_LogisticRegression_best.predict_proba(X_Test)
ŷ_test[:3]

Y_class = Library_LogisticRegression_best.predict(X_Test)
Y_class

# Re-lable feature headers and add new class prediction index column
new_colNames = ['{}_Probability'.format(c_name) for c_name in Training_Class] + ['ClassPredInd']
new_colNames

# Implemnt index array of probabilities
i_prob = np.concatenate((ŷ_test, Y_class[:, None]), 1)

# Create New dataframe for probality indeces
df2 = pd.DataFrame(i_prob, columns = new_colNames)
df2

# Concat dependant Ancestory features to dataframe
dep_preds = pd.concat([Test_dataset['Ancestry'], df2], axis = 1)

# Add new 
dep_preds['ClassPredName'] = dep_preds['ClassPredInd'].apply(lambda x: Training_Class[int(x)])

# Validate Probability predictions dataframe
dep_preds.head()

# Slice prediction and set new feature vector column variable
prob_1 = dep_preds.loc[:, 'Ancestry':'NativeAmerican_Probability']

# Unpivot convert dataFrame to long format
prob_2 = pd.melt(prob_1, id_vars = ['Ancestry'], var_name = 'Ancestry_Predictions', value_name = 'Probability')

# Test for true probability
prob_2['Ancestry_Predictions'] = prob_2['Ancestry_Predictions'].apply(lambda x: x.split('Prob')[0])

# Validate dataframe
prob_2.head(5)

# Validate dataframe features
print('Describe Columns:=', prob_2.columns, '\n')
print('Data Index values:=', prob_2.index, '\n')
print('Describe data:=', prob_2.describe(), '\n')

# Plot Probality prediction matrix
sns.set(rc = {'figure.figsize':(15,8)})
sns.set_theme(style="whitegrid")
fig, ax = plt.subplots()
sns.barplot(data = prob_2[prob_2['Ancestry'] != 'Unknown'],color = 'r', x = 'Ancestry', y = 'Probability', hue = 'Ancestry_Predictions', palette = 'mako')
plt.xlabel('Ancestory Classes')
plt.ylabel('Probability')
plt.suptitle('Probabilty of Ancestor classes')
#plt.savefig("Assignment3_Deliverable4.png")
plt.show()

"""# **Deliverable 7.5**  <a class="anchor" id="Deliverable_6.5"></a>
How do the class label probabilities differ for the Mexican and African American samples when compared to the class label probabilities for the unknown samples? Are these class probabilities telling us something about recent history? Explain why these class probabilities are reasonable with respect to knowledge of recent history?

- In comparison to the class label probabilities for the unknown samples, those with unknown ancestry show a probability close to or equal to one while the other classes show a probability close to zero or less than one. African American samples showed similar results. The model assigned high probabilities to the African ancestry class for each of these samples. However, both Native American and European ancestry contribute high probabilities to the Mexican population on average with Native American slightly higher than European.
"""