# load packages
import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import OneHotEncoder, OrdinalEncoder, MinMaxScaler
from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay, roc_auc_score
from sklearn.model_selection import train_test_split
import statsmodels.api as sm
from scipy.stats import kstest, norm

from imblearn.under_sampling import NearMiss, RandomUnderSampler

# To ignore warinings
import warnings
warnings.filterwarnings('ignore')

# load the training and test dataset
df_train = pd.read_csv('train.csv')
df_test = pd.read_csv('test.csv')

# check whether the training dataset has missing observations
print(df_train.isna().sum())
print("\n")
print(df_test.isna().sum())

# check whether there are duplicated observations
print(df_train.duplicated().sum())
print("\n")
print(df_test.duplicated().sum())

# show the first five rows of df_train
df_train.head()

# print the number of rows and columns in df_train and df_test
print("The training data have {} rows and {} columns.".format(df_train.shape[0], df_train.shape[1]))
print("\nThe test data have {} rows and {} columns.".format(df_test.shape[0], df_test.shape[1]))

# prints information about df_train
df_train.info()

df_train['Region_Code'] = df_train['Region_Code'].astype('int64')
df_test['Region_Code'] = df_test['Region_Code'].astype('int64')

df_train['Policy_Sales_Channel'] = df_train['Policy_Sales_Channel'].astype('int64')
df_test['Policy_Sales_Channel'] = df_test['Policy_Sales_Channel'].astype('int64')

# get a quick overview of the numeric data
df_train.describe().apply(lambda s: s.apply(lambda x: format(x, 'g')))

# get a quick overview of the categorical data
df_train.describe(include='object')

# how many unique values Gender, Vehicle_Age and Vehicle_Damange have?
df_train[['Gender', 'Vehicle_Age', 'Vehicle_Damage']].nunique()

# one-hot encode for Gender and Vehicle_Damage as they have 2 unique values
enc = OneHotEncoder(drop='if_binary', sparse_output=False, dtype=np.int64)

df_train['Gender'] = enc.fit_transform(df_train[['Gender']]) # 1 for male and 0 otherwise
df_test['Gender'] = enc.fit_transform(df_test[['Gender']])

df_train['Vehicle_Damage'] = enc.fit_transform(df_train[['Vehicle_Damage']]) # 1 for Yes and 0 otherwise
df_test['Vehicle_Damage'] = enc.fit_transform(df_test[['Vehicle_Damage']])

# ordinal encoding for Vehicle_Age
enc = OrdinalEncoder(categories=[['< 1 Year', '1-2 Year', '> 2 Years']], dtype=np.int64)

df_train['Vehicle_Age'] = enc.fit_transform(df_train[['Vehicle_Age']])
df_test['Vehicle_Age'] = enc.fit_transform(df_test[['Vehicle_Age']])

# display the frequency of each category to compare the counts across categories.
sns.countplot(data=df_train, x='Gender')

# show the proportion of each category
sns.countplot(data=df_train, x='Gender', stat='percent')

# display the frequency of each category to compare the counts across categories.
sns.countplot(data=df_train, x='Driving_License')

# show the proportion of each category
sns.countplot(data=df_train, x='Driving_License', stat='percent');

df_train['Driving_License'].value_counts(normalize=True)

# display the frequency of each category to compare the counts across categories.
plt.figure(figsize=(10, 6))
sns.countplot(data=df_train, x='Region_Code')
plt.xticks(rotation=45);

# show the proportion of each category
plt.figure(figsize=(10, 6))
sns.countplot(data=df_train, x='Region_Code', stat='percent');
plt.xticks(rotation=45);

# Summary Statistics
region_mode = df_train['Region_Code'].mode()[0]
unique_regions = df_train['Region_Code'].nunique()

print(f"Mode: {region_mode}")
print(f"Number of Unique Regions: {unique_regions}")

# display the frequency of each category to compare the counts across categories.
sns.countplot(data=df_train, x='Previously_Insured')

# show the proportion of each category
sns.countplot(data=df_train, x='Previously_Insured', stat='percent');

df_train['Previously_Insured'].value_counts(normalize=True)

# display the frequency of each category to compare the counts across categories.
sns.countplot(data=df_train, x='Vehicle_Age');

# show the proportion of each category
sns.countplot(data=df_train, x='Vehicle_Age', stat='percent');

df_train['Vehicle_Age'].value_counts(normalize=True)

# display the frequency of each category to compare the counts across categories.
sns.countplot(data=df_train, x='Vehicle_Damage');

# show the proportion of each category
sns.countplot(data=df_train, x='Vehicle_Damage', stat='percent');

# display the frequency of each category to compare the counts across categories.
plt.figure(figsize=(10, 6))
sns.countplot(data=df_train, x='Policy_Sales_Channel');
plt.xticks(rotation=45);

# show the proportion of each category
plt.figure(figsize=(10, 6))
sns.countplot(data=df_train, x='Policy_Sales_Channel', stat='percent');
plt.xticks(rotation=45);

# Summary Statistics
channel_mode = df_train['Policy_Sales_Channel'].mode()[0]
unique_channels = df_train['Policy_Sales_Channel'].nunique()

print(f"Mode: {channel_mode}")
print(f"Number of Unique Regions: {unique_channels}")

# display the frequency of each category to compare the counts across categories.
sns.countplot(data=df_train, x='Response');

# show the proportion of each category
sns.countplot(data=df_train, x='Response', stat='percent');

# the frequency distribution of a numeric variable
sns.histplot(data=df_train, x='Age');

# combine the distribution shape with summary statistics
sns.violinplot(data=df_train, x='Age')

# check if the data follows a normal distribution.
sm.qqplot(df_train['Age'], line ='45');

# the frequency distribution of a numeric variable
sns.histplot(data=df_train, x='Annual_Premium');

# combine the distribution shape with summary statistics
sns.violinplot(data=df_train, x='Annual_Premium')

# check if the data follows a normal distribution.
sm.qqplot(df_train['Annual_Premium'], line ='45');

# the frequency distribution of a numeric variable
sns.histplot(data=df_train, x='Vintage');

# combine the distribution shape with summary statistics
sns.violinplot(data=df_train, x='Vintage')

# check if the data follows a normal distribution.
sm.qqplot(df_train['Vintage'], line ='45');

plt.figure(figsize=(10, 10))
sns.heatmap(df_train.corr(), annot=True)

sns.boxplot(x='Vehicle_Age', y='Age', data=df_train)
plt.show()

sns.histplot(data=df_train, x='Age', hue='Vehicle_Age', palette=['red', 'blue', 'green']);

sns.countplot(x='Previously_Insured', hue='Vehicle_Damage', data=df_train);

g = sns.FacetGrid(df_train, col="Previously_Insured", row="Vehicle_Damage", margin_titles=True)
g.map(sns.countplot, "Response")
plt.show()

top7_policy_sales_channels = df_train['Policy_Sales_Channel'].value_counts()[:7].index
df_train_top7_policy_sales_channels = df_train[df_train['Policy_Sales_Channel'].isin(top7_policy_sales_channels)]

sns.boxplot(x='Policy_Sales_Channel', y='Age', data=df_train_top7_policy_sales_channels)
plt.show()

# normalize continuous features
scaler = MinMaxScaler()

df_train['Age'] = scaler.fit_transform(df_train[['Age']])
df_test['Age'] = scaler.transform(df_test[['Age']])

df_train['Annual_Premium'] = scaler.fit_transform(df_train[['Annual_Premium']])
df_test['Annual_Premium'] = scaler.transform(df_test[['Annual_Premium']])

df_train['Vintage'] = scaler.fit_transform(df_train[['Vintage']])
df_test['Vintage'] = scaler.transform(df_test[['Vintage']])

# normalize continuous features
df_train['Region_Code'] = scaler.fit_transform(df_train[['Region_Code']])
df_test['Region_Code'] = scaler.fit_transform(df_test[['Region_Code']])

df_train['Vehicle_Age'] = scaler.fit_transform(df_train[['Vehicle_Age']])
df_test['Vehicle_Age'] = scaler.fit_transform(df_test[['Vehicle_Age']])

df_train['Policy_Sales_Channel'] = scaler.fit_transform(df_train[['Policy_Sales_Channel']])
df_test['Policy_Sales_Channel'] = scaler.fit_transform(df_test[['Policy_Sales_Channel']])

y = df_train['Response']
X = df_train.drop(['Response', 'id'], axis=1)

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

class LinearLayer:
    """
        This Class implements all functions to be executed by a linear layer
        in a computational graph
        Args:
            input_shape: input shape of Data/Activations
            n_out: number of neurons in layer
            ini_type: initialization type for weight parameters, default is "plain"
                      Opitons are: plain, xavier and he
        Methods:
            forward(A_prev)
            backward(upstream_grad)
            update_params(learning_rate)
    """

    def __init__(self, input_shape, n_out, ini_type="plain"):
        """
        The constructor of the LinearLayer takes the following parameters
        Args:
            input_shape: input shape of Data/Activations
            n_out: number of neurons in layer
            ini_type: initialization type for weight parameters, default is "plain"
        """

        self.m = input_shape[1]  # number of examples in training data
        # `params` store weights and bias in a python dictionary
        self.params = initialize_parameters(input_shape[0], n_out, ini_type)  # initialize weights and bias
        self.Z = np.zeros((self.params['W'].shape[0], input_shape[1]))  # create space for resultant Z output

    def forward(self, A_prev):
        """
        This function performs the forwards propagation using activations from previous layer
        Args:
            A_prev:  Activations/Input Data coming into the layer from previous layer
        """

        self.A_prev = A_prev  # store the Activations/Training Data coming in
        self.Z = np.dot(self.params['W'], self.A_prev) + self.params['b']  # compute the linear function

    def backward(self, upstream_grad):
        """
        This function performs the back propagation using upstream gradients
        Args:
            upstream_grad: gradient coming in from the upper layer to couple with local gradient
        """

        # derivative of Cost w.r.t W
        self.dW = np.dot(upstream_grad, self.A_prev.T)

        # derivative of Cost w.r.t b, sum across rows
        self.db = np.sum(upstream_grad, axis=1, keepdims=True)

        # derivative of Cost w.r.t A_prev
        self.dA_prev = np.dot(self.params['W'].T, upstream_grad)

    def update_params(self, learning_rate=0.1):
        """
        This function performs the gradient descent update
        Args:
            learning_rate: learning rate hyper-param for gradient descent, default 0.1
        """

        self.params['W'] = self.params['W'] - learning_rate * self.dW  # update weights
        self.params['b'] = self.params['b'] - learning_rate * self.db  # update bias(es)

class SigmoidLayer:
    """
    This file implements activation layers
    inline with a computational graph model
    Args:
        shape: shape of input to the layer
    Methods:
        forward(Z)
        backward(upstream_grad)
    """

    def __init__(self, shape):
        """
        The consturctor of the sigmoid/logistic activation layer takes in the following arguments
        Args:
            shape: shape of input to the layer
        """
        self.A = np.zeros(shape)  # create space for the resultant activations

    def forward(self, Z):
        """
        This function performs the forwards propagation step through the activation function
        Args:
            Z: input from previous (linear) layer
        """
        self.A = 1 / (1 + np.exp(-Z))  # compute activations

    def backward(self, upstream_grad):
        """
        This function performs the  back propagation step through the activation function
        Local gradient => derivative of sigmoid => A*(1-A)
        Args:
            upstream_grad: gradient coming into this layer from the layer above
        """
        # couple upstream gradient with local gradient, the result will be sent back to the Linear layer
        self.dZ = upstream_grad * self.A*(1-self.A)

def initialize_parameters(n_in, n_out, ini_type='plain'):
    """
    Helper function to initialize some form of random weights and Zero biases
    Args:
        n_in: size of input layer
        n_out: size of output/number of neurons
        ini_type: set initialization type for weights
    Returns:
        params: a dictionary containing W and b
    """

    params = dict()  # initialize empty dictionary of neural net parameters W and b

    if ini_type == 'plain':
        params['W'] = np.random.randn(n_out, n_in) *0.01  # set weights 'W' to small random gaussian
    elif ini_type == 'xavier':
        params['W'] = np.random.randn(n_out, n_in) / (np.sqrt(n_in))  # set variance of W to 1/n
    elif ini_type == 'he':
        # Good when ReLU used in hidden layers
        # Delving Deep into Rectifiers: Surpassing Human-Level Performance on ImageNet Classification
        # Kaiming He et al. (https://arxiv.org/abs/1502.01852)
        # http: // cs231n.github.io / neural - networks - 2 /  # init
        params['W'] = np.random.randn(n_out, n_in) * np.sqrt(2/n_in)  # set variance of W to 2/n

    params['b'] = np.zeros((n_out, 1))    # set bias 'b' to zeros

    return params

def compute_stable_bce_cost(Y, Z):
    """
    This function computes the "Stable" Binary Cross-Entropy(stable_bce) Cost and returns the Cost and its
    derivative w.r.t Z_last(the last linear node) .
    The Stable Binary Cross-Entropy Cost is defined as:
    => (1/m) * np.sum(max(Z,0) - ZY + log(1+exp(-|Z|)))
    Args:
        Y: labels of data
        Z: Values from the last linear node
    Returns:
        cost: The "Stable" Binary Cross-Entropy Cost result
        dZ_last: gradient of Cost w.r.t Z_last
    """
    m = Y.shape[1]

    cost = (1/m) * np.sum(np.maximum(Z, 0) - Z*Y + np.log(1+ np.exp(- np.abs(Z))))
    dZ_last = (1/m) * ((1/(1+np.exp(- Z))) - Y)  # from Z computes the Sigmoid so P_hat - Y, where P_hat = sigma(Z)

    return cost, dZ_last

# define training constants
learning_rate = 1
number_of_epochs = 500

np.random.seed(48) # set seed value so that the results are reproduceable
                  # (weights will now be initailzaed to the same pseudo-random numbers, each time)


# Our network architecture has the shape: 
#                   (input)--> [Linear->Sigmoid] -> [Linear->Sigmoid] -->(output)  

#------ LAYER-1 ----- define hidden layer that takes in training data 
Z1 = LinearLayer(input_shape=(10, 2263750), n_out=10, ini_type='plain')
A1 = SigmoidLayer(Z1.Z.shape)

#------ LAYER-2 ----- define output layer that take is values from hidden layer
Z2 = LinearLayer(input_shape=A1.A.shape, n_out=1, ini_type='plain')
A2 = SigmoidLayer(Z2.Z.shape)

train_costs = [] # initially empty list, this will store all the costs after a certian number of epochs
test_costs = []

# Set up the undersampling method
undersampler = RandomUnderSampler()

# Start training
for epoch in range(number_of_epochs):

    # Apply the transformation to the dataset
    X_train_sampled, y_train_sampled = undersampler.fit_resample(X_train, y_train)
    X_train_sampled = np.array(X_train_sampled).T
    y_train_sampled = np.array(y_train_sampled).reshape(-1, 1).T
    
    # ------------------------- forward-prop -------------------------
    Z1.forward(X_train_sampled)
    A1.forward(Z1.Z)
    
    Z2.forward(A1.A)
    A2.forward(Z2.Z)
    
    # ---------------------- Compute Cost ----------------------------
    train_cost, dZ2 = compute_stable_bce_cost(y_train_sampled, Z2.Z)
    train_costs.append(train_cost)
    
    # ------------------------- back-prop ----------------------------
    Z2.backward(dZ2)

    A1.backward(Z2.dA_prev)
    Z1.backward(A1.dZ)

    # ---------------------- Forward pass on the test data -----------
    Z1.forward(np.array(X_test).T)
    A1.forward(Z1.Z)
    
    Z2.forward(A1.A)
    A2.forward(Z2.Z)
    
    test_cost, dZ_last = compute_stable_bce_cost(np.array(y_test).reshape(-1, 1).T, Z2.Z)
    test_costs.append(test_cost)
    
    # print and store Costs every 100 iterations.
    if (epoch % 100) == 0:

        print("Cost at epoch# {}: on the training data - {}, on the test data - {}".format(epoch, train_cost, test_cost))

    # ----------------------- Update weights and bias ----------------
    Z2.update_params(learning_rate=learning_rate)
    Z1.update_params(learning_rate=learning_rate)

plt.plot(train_costs, label='Training Cost')
plt.plot(test_costs, label='Test Cost')
plt.legend()
plt.show();

## accuracy on the training data
# forward pass
Z1.forward(X_train.T)
A1.forward(Z1.Z)
Z2.forward(A1.A)
A2.forward(Z2.Z)

y_train_preds = np.round(A2.A).flatten().astype('int64')
y_train_arr = np.array(y_train)
train_acc = np.sum(y_train_preds == y_train_arr) / len(y_train_arr)

print(f"The accuracy on the training data is: {np.round(100*train_acc, 2)}%.")

# accuracy on the test data
# forward pass
Z1.forward(X_test.T)
A1.forward(Z1.Z)
Z2.forward(A1.A)
A2.forward(Z2.Z)

y_test_preds = np.round(A2.A).flatten().astype('int64')
y_test_arr = np.array(y_test)
test_acc = np.sum(y_test_preds == y_test_arr) / len(y_test_arr)

print(f"The accuracy on the test data is: {np.round(100*test_acc, 2)}%.")

# confusion matrix
# cm = confusion_matrix(y_test_arr, y_test_preds)
ConfusionMatrixDisplay.from_predictions(y_test_arr, y_test_preds, values_format='.0f')
plt.show()

# roc auc score
print(f"ROC AUC Score: {roc_auc_score(y_test_arr, y_test_preds)}")

Z1.forward(np.array(df_test.drop(['id'], axis=1)).T)
A1.forward(Z1.Z)
Z2.forward(A1.A)
A2.forward(Z2.Z)

submission = pd.DataFrame({'id': df_test['id'], 'Response': A2.A.flatten()})
submission.to_csv('insurance_nn.csv', index=False)