e_commerce_churn.py

# -*- coding: utf-8 -*-
"""E Commerce Churn.ipynb

Automatically generated by Colab.

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

**Problem statement** :
Online E-Commerce company XYZ has been observing that customers  switching to competitor banks over the past couple of quarters. As such, this has caused a huge dent in the quarterly revenues and might drastically affect annual revenues for the ongoing financial year, causing stocks to plunge and market cap to reduce by X %.  A team of experts from the fields of economics, product development, technology and data science was assembled to halt this decline.

**Objective** : Is it possible to use a model to predict with reasonable accuracy which customers will leave in the near future? The ability to accurately estimate the timing of their departure would be an added bonus.

**Definition of churn** : A customer having closed all their active accounts with the e commerce paltform is said to have churned. Churn can be defined in other ways as well, based on the context of the problem. A customer not transacting for 6 months or 1 year can also be defined as to have churned, based on the business requirements

From the perspective of a business team/product manager: :

(1) Business goal : Arrest slide in revenues or loss of active E-commerce customers

(2) Identify data source : Transactional systems, event-based logs, Data warehouse (MySQL DBs, Redshift/AWS), Data Lakes, NoSQL DBs

(3) Audit for data quality : De-duplication of events/transactions, Complete or partial absence of data for chunks of time in between, Obscuring PII (personal identifiable information) data

(4) Define business and data-related metrics : Tracking of these metrics over time, probably through some intuitive visualizations

(i) Business metrics : Churn rate (month-on-month, weekly/quarterly), Trend of avg. number of products per customer,
    %age of dormant customers, Other such descriptive metrics

(ii) Data-related metrics : F1-score, Recall, Precision
     
     Recall = TP/(TP + FN)
     
     Precision = TP/(TP + FP)
     
     F1-score = Harmonic mean of Recall and Precision
     
     where, TP = True Positive, FP = False Positive and FN = False Negative

(5) Prediction model output format : Since this is not going to be an online model, it doesn't require deployment. Instead, periodic (monthly/quarterly) model runs could be made and the list of customers, along with their propensity to churn shared with the business (Sales/Marketing) or Product team

(6) Action to be taken based on model's output/insights : Based on the output obtained from Data Science team as above, various business interventions can be made to save the customer from getting churned. Customer-centric E-commerce offers, getting in touch with customers to address grievances etc. Here, also Data Science team can help with basic EDA to highlight different customer groups/segments and the appropriate intervention to be applied against them

__Collaboration with Engineering and DevOps :__  

(1) Application deployment on production servers (In the context of this problem statement, not required)

(2) [DevOps] Monitoring the scale aspects of model performance over time (Again, not required, in this case)

### How to set the target/goal for the metrics?

* Data science-related metrics :
    - Recall : >70%
    - Precision : >70%
    - F1-score : >70%


* Business metrics : Usually, it's top down. But a good practice is to consider it to make at least half the impact of the data science metric. For e.g., If we take Recall target as __70%__ which means correctly identifying 70% of customers who's going to churn in the near future, we can expect that due to business intervention (offers, getting in touch with customers etc.), 50% of the customers can be saved from being churned, which means atleast a __35%__ improvement in Churn Rate

## Show me the code!
"""

from google.colab import drive
drive.mount('/content/drive')

!pip install "numpy<2.0" "xgboost==2.1.4" "scikit-learn==1.3.2"

import numpy
import xgboost
print(numpy.__version__)    # Should be 1.26.4
print(xgboost.__version__)  # Should be 2.1.4

## Import required libraries

import os

import pandas as pd
from pandas.plotting import parallel_coordinates

import numpy as np

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from matplotlib.legend_handler import HandlerPathCollection
import matplotlib.lines as mlines

import seaborn as sns

from sklearn.impute import SimpleImputer
from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer
from sklearn.neighbors import LocalOutlierFactor
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import MinMaxScaler
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score, f1_score, recall_score, confusion_matrix, classification_report
from sklearn.tree import DecisionTreeClassifier
from sklearn.model_selection import cross_val_score
from sklearn.ensemble import RandomForestClassifier, ExtraTreesClassifier
from sklearn.model_selection import RandomizedSearchCV, train_test_split
from sklearn.pipeline import Pipeline
from sklearn.metrics import accuracy_score, classification_report, confusion_matrix, roc_curve, auc
from sklearn.calibration import CalibratedClassifierCV
from sklearn.calibration import CalibrationDisplay
from sklearn.metrics import brier_score_loss
from sklearn.neighbors import KNeighborsClassifier
from sklearn.metrics import precision_recall_curve
from sklearn.compose import ColumnTransformer
from sklearn.inspection import PartialDependenceDisplay


from lightgbm import LGBMClassifier
import xgboost as xgb
from xgboost import XGBClassifier

# Commented out IPython magic to ensure Python compatibility.
# %matplotlib inline

## Get multiple outputs in the same cell
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

## Ignore all warnings
import warnings
warnings.filterwarnings('ignore')
warnings.filterwarnings(action='ignore', category=DeprecationWarning)

## Display all rows and columns of a dataframe instead of a truncated version
from IPython.display import display
pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', None)

## Reading the dataset
# This might be present in S3, or obtained through a query on a database
df_churn = pd.read_csv('/content/drive/MyDrive/Churn Projectpro/E Commerce Dataset.csv')

# Dataframe has 16890 observations and 21 columns
print(df_churn.shape)

df_churn.head(10).T

"""### Basic EDA"""

df_churn.describe() # Describe all numerical columns
df_churn.describe(include = ['O']) # Describe all non-numerical/categorical columns

## Checking number of unique customers in the dataset
df_churn.shape[0], df_churn.CustomerID.nunique()

customer_counts = df_churn['CustomerID'].value_counts()
print(customer_counts)

customer_counts = df_churn['CustomerID'].value_counts()

# Check if the minimum and maximum counts are both 3
min_count = customer_counts.min()
max_count = customer_counts.max()

print(f"Minimum observations per customer: {min_count}")
print(f"Maximum observations per customer: {max_count}")

# List of all columns expected to be consistent for a given CustomerID across observations
# Exclude 'CustomerID' itself and 'Unnamed: 0' as they are identifiers or index-like.
all_static_cols = [col for col in df_churn.columns if col not in ['CustomerID', 'Unnamed: 0']]

# This produces a DataFrame showing the number of unique values for each of these columns per CustomerID
consistency_df = df_churn.groupby('CustomerID')[all_static_cols].nunique()

# Check if all customers have a unique count of exactly 1 across all relevant columns
all_consistent_overall = (consistency_df == 1).all().all()

if all_consistent_overall:
    print(f"\nAll customers have consistent values for all specified static attributes across their 3 observations.")
else:
    print(f"\nSome customers have different values across their observations in one or more of the following columns: {all_static_cols}")
    print("\nRun the next cell to see details per column.")

for col in all_static_cols:
    inconsistent_customers_for_col = consistency_df[consistency_df[col] != 1].index.tolist()
    if inconsistent_customers_for_col:
        print(f"\nCustomers with inconsistent values for '{col}': {len(inconsistent_customers_for_col)} customers")
        print(f"Example CustomerIDs for '{col}': {inconsistent_customers_for_col[:5]} (showing first 5 if many)")
        # Optionally, display the actual inconsistent data for the first customer in that column
        if inconsistent_customers_for_col:
            first_inconsistent_customer = inconsistent_customers_for_col[0]
            print(f"Data for CustomerID {first_inconsistent_customer} in column '{col}':")
            display(df_churn[df_churn['CustomerID'] == first_inconsistent_customer][['CustomerID', col]])

check_customer_data = df_churn[df_churn['CustomerID'] == 50011]
display(check_customer_data)

df_churn_ = df_churn.copy()
df_churn_ = df_churn_.drop_duplicates(subset=['CustomerID'])
# Drop the redundant Churn column and rename the remaining one
#df_churn_ = df_churn_.drop(columns=['Churn_y'])
#df_churn_ = df_churn_.rename(columns={'Churn_x': 'Churn'})
df_churn_.shape[0], df_churn.CustomerID.nunique()

df_churn = df_churn.drop_duplicates(subset=['CustomerID'])
# Drop the redundant Churn column and rename the remaining one
#df_churn = df_churn.drop(columns=['Churn_y'])
df_churn.shape[0], df_churn.CustomerID.nunique()

# The two-step merge process is the robust way to ensure everything works correctly.

# 1. Aggregate ONLY numeric columns using median
# Any customer who had ALL 3 rows as NaN will have a NaN here.
# df_agg_numeric = df_churn_.groupby('CustomerID').median(numeric_only=True).reset_index()

# 2. Aggregate categorical columns using 'first'
#df_agg_categorical = df_churn_.groupby('CustomerID')[['PreferredLoginDevice', 'PreferredPaymentMode',
#                                                    'Gender', 'PreferedOrderCat', 'MaritalStatus', 'Churn']].first().reset_index()

# 3. Merge the results
#df_churn_ = pd.merge(df_agg_numeric, df_agg_categorical, on='CustomerID', how='left')

# Drop the redundant Churn column and rename the remaining one
#df_churn_ = df_churn_.drop(columns=['Churn_y'])
#df_churn_ = df_churn_.rename(columns={'Churn_x': 'Churn'})

# Check how many NaNs are in your final 5630-row DataFrame:
print("\nNumber of NaNs remaining in final DataFrame:")
print(df_churn_.isnull().sum())

# Check how many NaNs are in your final 5630-row DataFrame:
print("\nNumber of NaNs remaining in final DataFrame:")
print(df_churn.isnull().sum())

df_tier = df_churn_.groupby(['CityTier']).agg({'CustomerID':'count', 'Churn':'mean'}
                                  ).reset_index().sort_values(by='CustomerID', ascending=False)

df_tier

df_churn_.CityTier.value_counts(normalize=True)

"""#### Conclusion

* Discard Unnamed: 0
* Discard CustomerID as well, since it doesn't convey any extra info. Each row pertains to a unique customer  
* Based on the above, columns/features can be segregated into non-essential, numerical, categorical and target variables
In general, CustomerID is a very useful feature on the basis of which we can calculate a lot of user-centric features. Here, the dataset is not sufficient to calculate any extra customer features
"""

## Separating out different columns into various categories as defined above
## is CityTier, (NumberOfDeviceRegistered) also categorical feature
target_var = ['Churn']
cols_to_remove = ['Unnamed: 0', 'CustomerID']
num_feats = ['Tenure',	'CityTier',	'WarehouseToHome', 'HourSpendOnApp', 'NumberOfDeviceRegistered', 'SatisfactionScore', 'NumberOfAddress', 'Complain', 'OrderAmountHikeFromlastYear', 'CouponUsed', 'OrderCount', 'DaySinceLastOrder', 'CashbackAmount']
cat_feats = ['PreferredLoginDevice', 'PreferredPaymentMode', 'Gender', 'PreferedOrderCat', 'MaritalStatus']

"""Among these, NumberOfDeviceRegistered could also be categorical and	Complain is binary categorical variable and SatisfactionScore probably ordinal.


"""

## Separating out target variable and removing the non-essential columns
y = df_churn_[target_var].values
df_churn_.drop(cols_to_remove, axis=1, inplace=True)

"""### Questioning the data :

 - Why do values in the origional dataset appear three times. An uplaod error?
 - No date/time column. A lot of useful features can be built using date/time columns
 - When was the data snapshot taken? There are certain customer features like : CashbackAmount, Tenure, OrderCount etc., which will have different values across time
 - Are all these values/features pertaining to the same single date or spread across multiple dates?
 - How frequently are customer features updated?
 - Will it be possible to have the values of these features over a period of time as opposed to a single, snapshot date?
 - Some customers have more than one address, do they order for others or for a company?
 - Customer transaction patterns can also help us ascertain whether the customer has actually churned or not. For example, a customer might transact daily/weekly vs a customer who transacts annuallly

 Here, the objective is to understand the data and distill the problem statement and the stated goal further. In the process, if more data/context can be obtained, that adds to the end result of the model performance

### Data Cleaning
"""

# Calculating the percentage of the missing values in every column
percent_missing = (df_churn_.isnull().sum() * 100) / len(df_churn_)

fig = plt.figure(figsize=(12, 6))
_ = percent_missing.plot(kind='bar')
_ = plt.title('Percentage of Missing Values per Column')
_ = plt.xlabel('Column Name')
_ = plt.ylabel('Percentage Missing')
_ = plt.xticks(rotation=45, ha='right')
plt.tight_layout()
plt.show()

print("Checking all columns for negative values:")

# Iterate through all columns in the DataFrame
for col in df_churn_.columns:
    # We only want to check columns that are numeric types (int or float)
    if np.issubdtype(df_churn_[col].dtype, np.number):
        # Check if any value is less than zero
        count_negative = (df_churn_[col] < 0).sum()

        if count_negative > 0:
            print(f"Column '{col}' has {count_negative} negative value(s).")
        # else:
        #     print(f"Column '{col}' is clean (no negatives).")

print("\nCheck complete.")

# Before imputation
print("Value counts for 'PreferredPaymentMode' before imputation:")
display(df_churn_['PreferredPaymentMode'].value_counts(dropna=False))

# Fill missing values in 'PreferredPaymentMode' with 'Unknown'
df_churn_['PreferredPaymentMode'].fillna('Unknown', inplace=True)

# Verify that there are no more NaNs in this column
print(f"\nNumber of NaN values in 'PreferredPaymentMode' after imputation: {df_churn_['PreferredPaymentMode'].isnull().sum()}")

df_churn_['PreferredPaymentMode'] = df_churn_['PreferredPaymentMode'].replace({
    'CC': 'Credit Card',
    'COD': 'Cash on Delivery'
})

df_churn['PreferredPaymentMode'] = df_churn['PreferredPaymentMode'].replace({
    'CC': 'Credit Card',
    'COD': 'Cash on Delivery'
})

# Example 1: PreferredPaymentMode
print("Value counts for 'PreferredPaymentMode':")
display(df_churn_['PreferredPaymentMode'].value_counts())

plt.figure(figsize=(8, 5))
sns.countplot(data=df_churn_, x='PreferredPaymentMode', palette='viridis')
_ = plt.title('Distribution of PreferredPaymentMode')
_ = plt.xlabel('PreferredPaymentMode')
_ = plt.ylabel('Count')
_ = plt.xticks(rotation=45, ha='right') # Added rotation here
plt.show()

# Before imputation
print("Value counts for 'MaritalStatus' before imputation:")
display(df_churn_['MaritalStatus'].value_counts(dropna=False))

# Fill missing values in 'MaritalStatus' with 'Unknown'
df_churn_['MaritalStatus'].fillna('Unknown', inplace=True)

# After imputation
print("\nValue counts for 'MaritalStatus' after imputation:")
display(df_churn_['MaritalStatus'].value_counts(dropna=False))

# Verify that there are no more NaNs in this column
print(f"\nNumber of NaN values in 'MaritalStatus' after imputation: {df_churn_['MaritalStatus'].isnull().sum()}")

# Example 1: MaritalStatus
print("Value counts for 'MaritalStatus':")
display(df_churn_['MaritalStatus'].value_counts())

plt.figure(figsize=(8, 5))
sns.countplot(data=df_churn_, x='MaritalStatus', palette='viridis')
_ = plt.title('Distribution Marital Status')
_ = plt.xlabel('MaritalStatus')
_ = plt.ylabel('Count')
_ = plt.xticks(rotation=45, ha='right') # Added rotation here
plt.show()

rows_with_zeros = df_churn_[df_churn_['PreferredLoginDevice'] == '0']
display(rows_with_zeros)

# Before replacement
print("Value counts for 'PreferredLoginDevice' before standardization:")
display(df_churn_['PreferredLoginDevice'].value_counts(dropna=False))

# Define the mapping
PreferredLoginDevice_mapping = {
    '0': 'Unknown'
}

# Apply the replacement
df_churn_['PreferredLoginDevice'] = df_churn_['PreferredLoginDevice'].replace(PreferredLoginDevice_mapping)

# After replacement
print("\nValue counts for 'PreferredLoginDevice' after standardization:")
display(df_churn_['PreferredLoginDevice'].value_counts(dropna=False))

# PreferredLoginDevice
print("Value counts for 'PreferredLoginDevice':")
display(df_churn_['PreferredLoginDevice'].value_counts())

plt.figure(figsize=(8, 5))
sns.countplot(data=df_churn_, x='PreferredLoginDevice', palette='viridis')
_ = plt.title('Distribution PreferredLoginDevice')
_ = plt.xlabel('Preferred PreferredLoginDevice')
_ = plt.ylabel('Count')
_ = plt.xticks(rotation=45, ha='right') # Added rotation here
plt.show()

# df_churn contains missing values
df_churn['PreferredLoginDevice'] = df_churn['PreferredLoginDevice'].replace("0", np.nan)

# df_churn contains missing values
print("Value counts for 'PreferredLoginDevice':")
display(df_churn['PreferredLoginDevice'].value_counts())

plt.figure(figsize=(8, 5))
sns.countplot(data=df_churn, x='PreferredLoginDevice', palette='viridis')
_ = plt.title('Distribution PreferredLoginDevice')
_ = plt.xlabel('Preferred PreferredLoginDevice')
_ = plt.ylabel('Count')
_ = plt.xticks(rotation=45, ha='right') # Added rotation here
plt.show()

# Before replacement
print("Value counts for 'Gender' before treatment:")
display(df_churn_['Gender'].value_counts(dropna=False))

# Define the mapping
gender_mapping = {
    'm': 'Male',
    'f': 'Female'
}

# Apply the replacement
df_churn_['Gender'] = df_churn_['Gender'].replace(gender_mapping)

# After replacement
print("\nValue counts for 'Gender' after treatment:")
display(df_churn_['Gender'].value_counts(dropna=False))

# df_churn chainging m to Male and f to Female
print("Value counts for 'Gender' before treatment")
display(df_churn['Gender'].value_counts(dropna=False))

# Define the mapping
gender_mapping = {
    'm': 'Male',
    'f': 'Female'
}

# Apply the replacement
df_churn['Gender'] = df_churn['Gender'].replace(gender_mapping)

# After replacement
print("\nValue counts for 'Gender' after treatment:")
display(df_churn['Gender'].value_counts(dropna=False))

rows_with_unreasonable_values = df_churn_[df_churn_['SatisfactionScore'] == 589314.0]
display(rows_with_unreasonable_values)

count_unreasonable_SatisfactionScore = (df_churn_['SatisfactionScore'] == -1).sum()
print(f"Number of rows where 'SatisfactionScore' is 589314.0: {count_unreasonable_SatisfactionScore}")

# Example 1: SatisfactionScore
print("Value counts for 'SatisfactionScore':")
display(df_churn_['SatisfactionScore'].value_counts())

plt.figure(figsize=(8, 5))
sns.countplot(data=df_churn_, x='SatisfactionScore', palette='viridis')
_ = plt.title('Distribution SatisfactionScore')
_ = plt.xlabel('SatisfactionScore')
_ = plt.ylabel('Count')
_ = plt.xticks(rotation=45, ha='right') # Added rotation here
plt.show()

df_churn['SatisfactionScore'] = df_churn['SatisfactionScore'].replace(589314.0, np.nan)

# Example 1: SatisfactionScore
print("Value counts for 'SatisfactionScore':")
display(df_churn['SatisfactionScore'].value_counts())

plt.figure(figsize=(8, 5))
sns.countplot(data=df_churn, x='SatisfactionScore', palette='viridis')
_ = plt.title('Distribution SatisfactionScore')
_ = plt.xlabel('SatisfactionScore')
_ = plt.ylabel('Count')
_ = plt.xticks(rotation=45, ha='right') # Added rotation here
plt.show()

CategoricalFeatures = ['PreferredLoginDevice', 'PreferredPaymentMode', 'Gender', 'PreferedOrderCat', 'MaritalStatus',
                       "SatisfactionScore"]

def ValueCounts():
  for i in CategoricalFeatures:
    print(df_churn_[i].value_counts())

ValueCounts()

# Count the total number of entries where the value in 'Tenure' is less than 0
count_negative_values = (df_churn_['Tenure'] < 0).sum()

print(f"Total number of negative values in df_churn_['Tenure']: {count_negative_values}")

rows_with_negative_tenure = df_churn_[df_churn_['Tenure'] == -10000]
display(rows_with_negative_tenure)

count_negative_tenure = (df_churn_['Tenure'] == -10000).sum()
print(f"Number of rows where 'Tenure' is -10000: {count_negative_tenure}")

count_nan_tenure = df_churn_['Tenure'].isnull().sum()
print(f"Number of NaN values in 'Tenure': {count_nan_tenure}")

df_churn_['Tenure'] = df_churn_['Tenure'].replace(-10000, np.nan)

# Let's verify the change in the new DataFrame
count_negative_tenure_new_df = (df_churn_['Tenure'] == -10000).sum()
count_nan_tenure_new_df = df_churn_['Tenure'].isnull().sum()

print(f"In the new DataFrame, 'df_churn_':")
print(f"  Number of -10000 values in 'Tenure': {count_negative_tenure_new_df}")
print(f"  Number of NaN values in 'Tenure' after replacement: {count_nan_tenure_new_df}")

# Also show a sample of the updated column
display(df_churn_[df_churn_['Tenure'].isnull()].head())

# Counting Nan values in df_churn

count_nan_tenure_in_df_churn = df_churn['Tenure'].isnull().sum()
print(f"Number of NaN values in 'Tenure': {count_nan_tenure_in_df_churn}")

# Veriyfying that -10000 has changed to NaN in df_chrun

df_churn['Tenure'] = df_churn['Tenure'].replace(-10000, np.nan)

# Let's verify the change in the new DataFrame
count_negative_tenure_new_df = (df_churn['Tenure'] == -10000).sum()
count_nan_tenure_new_df = df_churn['Tenure'].isnull().sum()

print(f"In the new DataFrame, 'df_churn':")
print(f"  Number of -10000 values in 'Tenure': {count_negative_tenure_new_df}")
print(f"  Number of NaN values in 'Tenure' after replacement: {count_nan_tenure_new_df}")

# Also show a sample of the updated column
display(df_churn[df_churn['Tenure'].isnull()].head())

"""* The missing values of categorical features have been filled with Unkown and also a new category been generated for the missing values
* The gender feature included m and f. m has bee assigned to Male and f has been assigned to Female.
* Tenure had among NaN also -10000 which are now treated as NaN values.
* For the second churn df (churn_df) the missing values or unreasonable values have been defined as NaN

### Separating out train-test-valid sets

Since this is the only data available to us, we keep aside a holdout/test set to evaluate our model at the very end in order to estimate our chosen model's performance on unseen data / new data.

A validation set is also created which we'll use in our baseline models to evaluate and tune our models
"""

#Data Frames
#1. df_churn original dataframe contains duplicate entries (NOTE: DUPLICATES REMOVED BEFOR THE SPLIT)
#2. df_churn_ duplicate entries removed missing values in categorical features are categorized as unknown and
#unreasonable numbers in numeric features are defined as NaN
#3. df_churn --> df_val_m --> df_train_m --> df_test_m --> y_train_m --> y_val_m (NOTE: df_churn dubplicates removed, contains missing values)
#4. df_churn_ --> df_val --> df_train --> df_test --> y_train --> y_val --> y_test
#5. df_train, df_test, df_val -> df_train_II, df_test_II -> df_test_II, df_val_II (II = Imputed, No multicollinearity, encoded)
# df_train_m, df_val_m, df_test_m --> df_train_m_X, df_val_m_X, df_test_m_X
#6. df_train_II, df_X_test_II, df_X_val_II -> sc_en_X_train_II, sc_en_X_test_II, sc_en_X_val_II (sclaed and encoded)
#   df_train_II, df_X_test_II, df_X_val_II --> df_train_II_NoOutliers, df_test_II_NoOutliers, df_val_II_NoOutliers
#  sc_en_X_train_II, sc_en_X_test_II, sc_en_X_val_II --> sc_en_X_train_II_NoOutliers, sc_en_X_test_II_NoOutliers, sc_en_X_val_II_NoOutliers
#7. X_train, X_val, X_test -> df_train_II, df_val_II, df_test_II (Imputed, No multicollinearity, encoded)
#8. X_train_m, X_val_m, X_test_m -> df_train_m, df_val_m, df_test_m

## Keeping aside a test/holdout set
df_train_val, df_test, y_train_val, y_test = train_test_split(df_churn_, y.ravel(), test_size = 0.1, random_state = 42)

## Splitting into train and validation set
df_train, df_val, y_train, y_val = train_test_split(df_train_val, y_train_val, test_size = 0.12, random_state = 42)

df_train.shape, df_val.shape, df_test.shape, y_train.shape, y_val.shape, y_test.shape
np.mean(y_train), np.mean(y_val), np.mean(y_test)

## Applying on df_churn that has missing values
## Keeping aside a test/holdout set
df_train_val_m, df_test_m, y_train_val_m, y_test_m = train_test_split(df_churn, y.ravel(), test_size = 0.1, random_state = 42)

## Splitting into train and validation set
df_train_m, df_val_m, y_train_m, y_val_m = train_test_split(df_train_val_m, y_train_val_m, test_size = 0.12, random_state = 42)

# 'df' is your cleaned dataframe
#df_churn.to_csv('df_churn.csv', index=False)

#from google.colab import files
#files.download('df_churn.csv')

df_train_m_X = df_train_m.copy()
df_val_m_X = df_val_m.copy()
df_test_m_X = df_test_m.copy()

df_train_m_X.columns

sns.set(style="whitegrid")

# Create a count plot for the 'Churn' variable
plt.figure(figsize=(6, 4)) # Optional: Adjusts the size of the plot
sns.countplot(x="Churn", data=df_train) # 'x' maps to the categorical column
plt.title("Distribution of Churn df_train (No Churn vs Churn)") # Add a descriptive title
plt.xlabel("Churn Status") # Optional: label for the x-axis
plt.ylabel("Count") # Optional: label for the y-axis
plt.show() # Displays the plot

sns.set(style="whitegrid")

# Create a count plot for the 'Churn' variable
plt.figure(figsize=(6, 4)) # Optional: Adjusts the size of the plot
sns.countplot(x="Churn", data=df_test) # 'x' maps to the categorical column
plt.title("Distribution of Churn df_test (No Churn vs Churn)") # Add a descriptive title
plt.xlabel("Churn Status") # Optional: label for the x-axis
plt.ylabel("Count") # Optional: label for the y-axis
plt.show() # Displays the plot

sns.set(style="whitegrid")

# Create a count plot for the 'Churn' variable
plt.figure(figsize=(6, 4)) # Optional: Adjusts the size of the plot
sns.countplot(x="Churn", data=df_val) # 'x' maps to the categorical column
plt.title("Distribution of Churn df_val (No Churn vs Churn)") # Add a descriptive title
plt.xlabel("Churn Status") # Optional: label for the x-axis
plt.ylabel("Count") # Optional: label for the y-axis
plt.show() # Displays the plot

sns.set(style="whitegrid")

# Create a count plot for the 'Churn' variable
plt.figure(figsize=(6, 4)) # Optional: Adjusts the size of the plot
sns.countplot(x="Churn", data=df_train_m) # 'x' maps to the categorical column
plt.title("Distribution of Churn df_train_m (No Churn vs Churn)") # Add a descriptive title
plt.xlabel("Churn Status") # Optional: label for the x-axis
plt.ylabel("Count") # Optional: label for the y-axis
plt.show() # Displays the plot

sns.set(style="whitegrid")

# Create a count plot for the 'Churn' variable
plt.figure(figsize=(6, 4)) # Optional: Adjusts the size of the plot
sns.countplot(x="Churn", data=df_test_m) # 'x' maps to the categorical column
plt.title("Distribution of Churn df_test_m (No Churn vs Churn)") # Add a descriptive title
plt.xlabel("Churn Status") # Optional: label for the x-axis
plt.ylabel("Count") # Optional: label for the y-axis
plt.show() # Displays the plot

sns.set(style="whitegrid")

# Create a count plot for the 'Churn' variable
plt.figure(figsize=(6, 4)) # Optional: Adjusts the size of the plot
sns.countplot(x="Churn", data=df_val_m) # 'x' maps to the categorical column
plt.title("Distribution of Churn df_val_m (No Churn vs Churn)") # Add a descriptive title
plt.xlabel("Churn Status") # Optional: label for the x-axis
plt.ylabel("Count") # Optional: label for the y-axis
plt.show() # Displays the plot

"""### Univariate plots of numerical variables in training set"""

num_list = ['Tenure', 'WarehouseToHome', 'HourSpendOnApp', 'NumberOfDeviceRegistered',
            'NumberOfAddress', 'OrderAmountHikeFromlastYear','CouponUsed', 'OrderCount',
            'DaySinceLastOrder', 'CashbackAmount']
def numeric_features_visuals():
  for visual in num_list:
    plt.figure()
    sns.set(style = "whitegrid")
    sns.boxplot(y = df_train[visual])
  return plt.show()

print(numeric_features_visuals())

num_list = ['Tenure', 'WarehouseToHome', 'HourSpendOnApp', 'NumberOfDeviceRegistered',
            'NumberOfAddress', 'OrderAmountHikeFromlastYear','CouponUsed', 'OrderCount',
            'DaySinceLastOrder', 'CashbackAmount']
def numeric_features_visuals():
  for visual in num_list:
    plt.figure()
    #sns.set(style = "whitegrid")
    sns.violinplot(y = df_train[visual])
  return plt.show()

print(numeric_features_visuals())

"""Violin plot suggests that NumberOfDeviceRegistered and HoursSpendOnApp could be treated as categorical features as well?"""

num_list = ['Tenure', 'WarehouseToHome', 'HourSpendOnApp', 'NumberOfDeviceRegistered',
            'NumberOfAddress', 'OrderAmountHikeFromlastYear','CouponUsed', 'OrderCount',
            'DaySinceLastOrder', 'CashbackAmount']
def numeric_features_visuals():
  for visual in num_list:
    plt.figure()
    sns.set(style = 'ticks')
    sns.distplot(df_train[visual], hist=True, kde=False)
  return plt.show()

print(numeric_features_visuals())

num_list = ['Tenure', 'WarehouseToHome', 'HourSpendOnApp', 'NumberOfDeviceRegistered',
            'NumberOfAddress', 'OrderAmountHikeFromlastYear','CouponUsed', 'OrderCount',
            'DaySinceLastOrder', 'CashbackAmount']
def numeric_features_visuals():
  for visual in num_list:
    plt.figure()
    #sns.set(style = 'ticks')
    sns.kdeplot(df_train[visual])
  return plt.show()

print(numeric_features_visuals())

num_list = ['Tenure', 'WarehouseToHome', 'HourSpendOnApp', 'NumberOfDeviceRegistered',
            'NumberOfAddress', 'OrderAmountHikeFromlastYear','CouponUsed', 'OrderCount',
            'DaySinceLastOrder', 'CashbackAmount']
def numeric_features_visuals():
  for visual in num_list:
    plt.figure()
    #sns.set(style = 'ticks')
    sns.histplot(df_train[visual])
  return plt.show()

print(numeric_features_visuals())

"""* HoursSpendOnApp and NumberOfDeviceRegistered might not be strong predictors because they are only 3 or 6 numbers (could be treated as categorical or even ordinal).
* CouponUsed and OrderCount might also not be strong predictrs due to fewer numbers in some categories.

### Bivariate plots of categorical variables in training set
"""

print("Unique categories in 'PreferredLoginDevice':")
display(df_train['PreferredLoginDevice'].unique())

print("\nCategories and their counts in 'PreferredPaymentMode':")
display(df_train['PreferredPaymentMode'].unique())

print("\nCategories and their counts in 'Gender':")
display(df_train['Gender'].unique())


print("\nCategories and their counts in 'PreferedOrderCat':")
display(df_train['PreferedOrderCat'].unique())


print("\nCategories and their counts in 'MaritalStatus':")
display(df_train['MaritalStatus'].unique())

print("\nCategories and their counts in 'SatisfactionScore':")
display(df_train['SatisfactionScore'].unique())

print("\nCategories and their counts in 'Complain':")
display(df_train['Complain'].unique())

print("\nCategories and their counts in 'CityTier':")
display(df_train['CityTier'].unique())

categorical_for_churn_viz = ['PreferredLoginDevice', 'PreferredPaymentMode', 'Gender', 'PreferedOrderCat',
            'MaritalStatus', 'SatisfactionScore', 'Complain', 'CityTier']
def categorical_vs_target_visuals():
  for visual in categorical_for_churn_viz:
    # Calculate counts for each category and Churn status
    temp_df = df_train.groupby(visual)['Churn'].value_counts().rename('count').reset_index()

    # Filter for 'Not Churned' (Churn == 0) and sort by their counts
    churn_0_df = temp_df[temp_df['Churn'].astype(int) == 0].sort_values(by='count', ascending=False)

    # Get the ordered list of categories to use in the plot
    order_list = churn_0_df[visual].tolist()

    plt.figure(figsize=(10, 6))
    sns.set(style = "whitegrid")
    sns.countplot(data=df_train, x=visual, hue='Churn', palette='viridis', order=order_list)
    plt.title(f'Distribution of {visual} by Churn (Ordered by Not Churned Count)')
    plt.xlabel(visual)
    plt.ylabel('Count')
    plt.xticks(rotation=45, ha='right')
    plt.legend(title='Churn', labels=['Not Churned', 'Churned'])
  return plt.show()

print(categorical_vs_target_visuals())

def categorical_vs_target_visuals_percentage_ordered():
    categorical_for_churn_viz = ['PreferredLoginDevice', 'PreferredPaymentMode', 'Gender', 'PreferedOrderCat',
                'MaritalStatus', 'SatisfactionScore', 'Complain', 'CityTier']

    for visual in categorical_for_churn_viz:
        fig, ax = plt.subplots(figsize=(10, 6))

        # Calculate churn percentages within each category
        temp_df = df_train.groupby(visual)['Churn'].value_counts(normalize=True).rename('percentage').reset_index()
        temp_df['percentage'] = temp_df['percentage'] * 100 # Convert to actual percentage

        # --- NEW STEP: Determine the order based on 'Not Churned' percentage (assuming 0 means Not Churned) ---
        # Filter for 'Not Churned' (Churn == 0) and sort in descending order of percentage
        # Use .astype(int) for Churn if it's currently stored as float or object
        churn_0_df = temp_df[temp_df['Churn'].astype(int) == 0].sort_values(by='percentage', ascending=False)

        # Get the ordered list of categories to use in the plot
        order_list = churn_0_df[visual].tolist()

        # Plot the bar chart using the calculated percentages and the new 'order' parameter
        sns.barplot(data=temp_df, x=visual, y='percentage', hue='Churn', palette='viridis', ax=ax, order=order_list)

        plt.title(f'Churn Percentage Distribution within {visual} (Ordered by Not Churned %)')
        plt.xlabel(visual)
        plt.ylabel('Percentage (%)')
        plt.xticks(rotation=45, ha='right')

        # --- Fix Legend ---
        handles, labels = ax.get_legend_handles_labels()
        # Ensure labels are correct if Churn values are 0/1 or string '0'/'1'
        ax.legend(handles, ['Not Churned', 'Churned'], title='Churn', bbox_to_anchor=(1.05, 1), loc='upper left')

        plt.tight_layout()

    return plt.show()

# Run the new function
print(categorical_vs_target_visuals_percentage_ordered())

# Assuming 'df_train' is your DataFrame and 'Churn' is 0/1 (0 for Not Churned, 1 for Churned)

def categorical_vs_target_visuals_percentage_ordered():
    categorical_for_churn_viz = ['HourSpendOnApp', 'NumberOfDeviceRegistered']

    for visual in categorical_for_churn_viz:
        fig, ax = plt.subplots(figsize=(10, 6))

        # Calculate churn percentages within each category
        temp_df = df_train.groupby(visual)['Churn'].value_counts(normalize=True).rename('percentage').reset_index()
        temp_df['percentage'] = temp_df['percentage'] * 100 # Convert to actual percentage

        # --- NEW STEP: Determine the order based on 'Not Churned' percentage (assuming 0 means Not Churned) ---
        # Filter for 'Not Churned' (Churn == 0) and sort in descending order of percentage
        # Use .astype(int) for Churn if it's currently stored as float or object
        churn_0_df = temp_df[temp_df['Churn'].astype(int) == 0].sort_values(by='percentage', ascending=False)

        # Get the ordered list of categories to use in the plot
        order_list = churn_0_df[visual].tolist()

        # Plot the bar chart using the calculated percentages and the new 'order' parameter
        sns.barplot(data=temp_df, x=visual, y='percentage', hue='Churn', palette='viridis', ax=ax, order=order_list)

        plt.title(f'Churn Percentage Distribution within {visual} (Ordered by Not Churned %)')
        plt.xlabel(visual)
        plt.ylabel('Percentage (%)')
        plt.xticks(rotation=45, ha='right')

        # --- Fix Legend ---
        handles, labels = ax.get_legend_handles_labels()
        # Ensure labels are correct if Churn values are 0/1 or string '0'/'1'
        ax.legend(handles, ['Not Churned', 'Churned'], title='Churn', bbox_to_anchor=(1.05, 1), loc='upper left')

        plt.tight_layout()

    return plt.show()

# Run the new function
print(categorical_vs_target_visuals_percentage_ordered())

"""Conlcusion

- HourSpendOnApp does not look like an informative feature. It might one of the least informative features in the dataset. There is no clear pattern.

- Regarding gender feature, both female and male the number of churn is the same for both gender.

- The other features show good pattern. For instance, the more somone is satisfied lessl likely that someone is going to churn. Also if someone complained then the churn likelyhood increases.

### Bivariate plots of numeric features
"""

num_list = ['Tenure', 'WarehouseToHome', 'HourSpendOnApp', 'NumberOfDeviceRegistered',
            'NumberOfAddress', 'OrderAmountHikeFromlastYear','CouponUsed', 'OrderCount',
            'DaySinceLastOrder', 'CashbackAmount']
def numeric_features_visuals():
  for visual in num_list:
    plt.figure(figsize=(8, 6)) # Adjust figure size for better readability
    sns.set(style = "whitegrid")
    sns.boxplot(x='Churn', y=visual, data=df_train, palette='viridis')
    plt.title(f'Distribution of {visual} by Churn')
    plt.xlabel('Churn (0 = No Churn, 1 = Churn)')
    plt.ylabel(visual)
    plt.grid(axis='y', linestyle='--', alpha=0.7)
  return plt.show()

print(numeric_features_visuals())

"""Cconlusion
- Regarding tenure, the longer somone is registered as a member, less likely is the churn.
- The more devices someone registred, the more likely someone is going churn.
- The more adresses someone has, the more likely someone is going to churn.
- The likelihood of churn is higher if the day since the previous order is lower.
- The higher the cash back amount, less likely someone is going to churn.
"""

# continue from here

"""### Missing values and outlier treatment

#### Outliers

* Can be observed from univariate plots of different features

* Outliers can either be logically improbable (as per the feature definition) or just an extreme value as compared to the feature distribution

* As part of outlier treatment, the particular row containing the outlier can be removed from the training set, provided they do not form a significant chunk of the dataset (< 0.5-1%)

* In cases where the value of outlier is logically faulty, e.g. negative Age or Tenure < 0, the particular record can be replaced with mean of the feature or the nearest among min/max logical value of the feature

Outliers in numerical features can be of a very high/low value, lying in the top 1% or bottom 1% of the distribution or values which are not possible as per the feature definition.

Outliers in categorical features are usually levels with a very low frequency/no. of samples as compared to other categorical levels.

##### Is outlier treatment always required ?

No, Not all ML algorithms are sensitive to outliers. Algorithms like linear/logistic regression are sensitive to outliers.

Tree algorithms, kNN, clustering algorithms etc. are in general, robust to outliers

Outliers affect metrics such as mean, std. deviation

#### Missing values

df_train
"""

# Calculating the percentage of the missing values in every column
percent_missing = (df_train.isnull().sum() * 100) / len(df_train)

fig = plt.figure(figsize=(12, 6))
_ = percent_missing.plot(kind='bar')
_ = plt.title('Percentage of Missing Values per Column')
_ = plt.xlabel('Column Name')
_ = plt.ylabel('Percentage Missing')
_ = plt.xticks(rotation=45, ha='right')
plt.tight_layout()
plt.show()

"""There are a few missing values which can also be observed from df.describe() commands. A couple of things which can be done in such cases :
 - If the column/feature has too many missing values, it can be dropped as it might not add much relevance to the data
 - If there a few missing values, the column/feature can be imputed with its summary statistics (mean/median/mode) and/or numbers like 0, -1 etc. which add value depending on the data and context.
"""

print(df_train.isnull().sum()/df_train.shape[0])

# Extract the 'WarehouseToHome' column into new, separate DataFrames
df_train_WarehouseToHome = df_train[['WarehouseToHome']].copy()
df_test_WarehouseToHome = df_test[['WarehouseToHome']].copy()
df_val_WarehouseToHome = df_val[['WarehouseToHome']].copy()

# The original df_train, df_test, and df_val remain untouched

# Generating a funciton for the specific column to asign them to df_train df_test and df_val
#def impute_missing_values(df_train_1, df_test_1, df_val_1):
 #   df_train_WarehouseToHome = df_train_1[['WarehouseToHome']]
 #   df_test_WarehouseToHome = df_test_1[['WarehouseToHome']]
 #   df_val_WarehouseToHome = df_val_1[['WarehouseToHome']]
 #   return df_train_1, df_test_1, df_val_1

#df_train_WarehouseToHome, df_test_WarehouseToHome, df_val_WarehouseToHome = impute_missing_values(df_train, df_test, df_val)

print(df_train.isnull().sum()/df_train.shape[0])
()
print(df_test.isnull().sum()/df_test.shape[0])
()
print(df_val.isnull().sum()/df_val.shape[0])

# Initialize the imputer once
imp = SimpleImputer(missing_values=np.nan, strategy='median')

# Impute missing values in the training set
df_train_WarehouseToHome['WarehouseToHome'] = imp.fit_transform(df_train_WarehouseToHome[['WarehouseToHome']])

# Impute missing values in the test set using the SAME fitted imputer
df_test_WarehouseToHome['WarehouseToHome'] = imp.transform(df_test_WarehouseToHome[['WarehouseToHome']])

# Impute missing values in the validation set using the SAME fitted imputer
df_val_WarehouseToHome['WarehouseToHome'] = imp.transform(df_val_WarehouseToHome[['WarehouseToHome']])

# Display samples of the imputed columns
print("df_train_WarehouseToHome after imputation:")
display(df_train_WarehouseToHome['WarehouseToHome'].head())
print("\ndf_test_WarehouseToHome after imputation:")
display(df_test_WarehouseToHome['WarehouseToHome'].head())
print("\ndf_val_WarehouseToHome after imputation:")
display(df_val_WarehouseToHome['WarehouseToHome'].head())

# Copying the dataframe so the original does not change
df_train_II = df_train.copy()
df_test_II = df_test.copy()
df_val_II = df_val.copy()

# Numeric Features to be imputed
# Check this again if Number of NumberOfAddress has missing values
numerical_features = ['Tenure', 'WarehouseToHome', 'HourSpendOnApp', 'NumberOfDeviceRegistered',
                      'OrderAmountHikeFromlastYear', 'CouponUsed' ,'OrderCount', 'DaySinceLastOrder']

# Select the relevant numerical data for the imputer
features_for_imputation_train = df_train_II[numerical_features].copy()
features_for_imputation_test = df_test_II[numerical_features].copy()
features_for_imputation_val = df_val_II[numerical_features].copy()


# Initialize a single imputer instance
imputer = IterativeImputer(sample_posterior=True, random_state=42)

# Fit the imputer ONLY on the training data and transform it
II_train_df = imputer.fit_transform(features_for_imputation_train)

# Transform the test and validation sets using the SAME fitted imputer
II_test_df = imputer.transform(features_for_imputation_test)
II_val_df = imputer.transform(features_for_imputation_val)

# Assign the imputed values back to the original DataFrame using the correct index
# The output is a NumPy array, so we create a temporary DataFrame to align indices/columns
df_train_II[numerical_features] = pd.DataFrame(II_train_df,
                                               columns=numerical_features,
                                               index=df_train_II.index)

df_test_II[numerical_features] = pd.DataFrame(II_test_df,
                                              columns=numerical_features,
                                              index=df_test_II.index)

df_val_II[numerical_features] = pd.DataFrame(II_val_df,
                                              columns=numerical_features,
                                              index=df_val_II.index)

# Check the missing values on the CORRECT DataFrame name
print(df_train_II[numerical_features].isnull().sum())
print(df_test_II[numerical_features].isnull().sum())
print(df_val_II[numerical_features].isnull().sum())

df_train_II.columns

df_train_II['WarehouseToHome'] = pd.to_numeric(df_train_II['WarehouseToHome'], errors='coerce')
df_test_II['WarehouseToHome'] = pd.to_numeric(df_test_II['WarehouseToHome'], errors='coerce')
df_val_II['WarehouseToHome'] = pd.to_numeric(df_val_II['WarehouseToHome'], errors='coerce')


plt.figure(figsize=(15, 10)) # Adjust overall figure size as needed
grid = gridspec.GridSpec(3, 2)

# Plot 1: df_train_WarehouseToHome
plt.subplot(grid[0, 0]) # Top-left subplot
sns.kdeplot(df_train.WarehouseToHome.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_train_WarehouseToHome.WarehouseToHome, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_train Before and After Missing Value Treatment')
plt.xlabel('WarehouseToHome')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 2: df_test_WarehouseToHome
plt.subplot(grid[0, 1]) # Top-right subplot
sns.kdeplot(df_test.WarehouseToHome.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_test_WarehouseToHome.WarehouseToHome, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_test Before and After Missing Value Treatment')
plt.xlabel('WarehouseToHome')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 3: df_val_WarehouseToHome (now in the bottom-left, same width as top plots)
plt.subplot(grid[1, 0]) # Bottom-left subplot
sns.kdeplot(df_val.WarehouseToHome.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_val_WarehouseToHome.WarehouseToHome, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Before and After Missing Value Treatment')
plt.xlabel('WarehouseToHome')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)


# Iterative Imputation
# Plot 4: df_val_WarehouseToHome (now in the bottom-left, same width as top plots)
plt.subplot(grid[1, 1]) # Bottom-left subplot
sns.kdeplot(df_train.WarehouseToHome.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_train_II.WarehouseToHome, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Before and After Missing Value Iterative Imputation')
plt.xlabel('WarehouseToHome')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 5: df_val_WarehouseToHome (now in the bottom-left, same width as top plots)
plt.subplot(grid[2, 0]) # Bottom-left subplot
sns.kdeplot(df_test.WarehouseToHome.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_test_II.WarehouseToHome, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_test Before and After Missing Value Iterative Imputation')
plt.xlabel('WarehouseToHome')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)


# Plot 5: df_val_WarehouseToHome (now in the bottom-left, same width as top plots)
plt.subplot(grid[2, 1]) # Bottom-left subplot
sns.kdeplot(df_val.WarehouseToHome.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_val_II.WarehouseToHome, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Before and After Missing Value Iterative Imputation')
plt.xlabel('WarehouseToHome')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)


plt.tight_layout() # Adjust layout to prevent overlapping titles/labels
plt.show()

"""Tenure"""

# Extract the 'Tenure' column into new, separate DataFrames
df_train_Tenure = df_train[['Tenure']]
df_test_Tenure = df_test[['Tenure']]
df_val_Tenure = df_val[['Tenure']]

# The original df_train, df_test, and df_val remain untouched

# Initialize the imputer once
imp = SimpleImputer(missing_values=np.nan, strategy='median')

# Impute missing values in the training set
df_train_Tenure['Tenure'] = imp.fit_transform(df_train_Tenure[['Tenure']])

# Impute missing values in the test set using the SAME fitted imputer
df_test_Tenure['Tenure'] = imp.transform(df_test_Tenure[['Tenure']])

# Impute missing values in the validation set using the SAME fitted imputer
df_val_Tenure['Tenure'] = imp.transform(df_val_Tenure[['Tenure']])

# Display samples of the imputed columns
print("df_train_Tenure after imputation:")
display(df_train_Tenure['Tenure'].head())
print("\ndf_test_Tenure after imputation:")
display(df_test_Tenure['Tenure'].head())
print("\ndf_val_Tenure after imputation:")
display(df_val_Tenure['Tenure'].head())

plt.figure(figsize=(15, 10)) # Adjust overall figure size as needed
grid = gridspec.GridSpec(3, 2)

# Plot 1: df_train_Tenure
plt.subplot(grid[0, 0]) # Top-left subplot
sns.kdeplot(df_churn_.Tenure.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_train_Tenure.Tenure, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_train Tenure Before and After Missing Value Treatment')
plt.xlabel('Tenure')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 2: df_test_Tenure
plt.subplot(grid[0, 1]) # Top-right subplot
sns.kdeplot(df_churn_.Tenure.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_test_Tenure.Tenure, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_test Tenure Before and After Missing Value Treatment')
plt.xlabel('Tenure')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 3: df_val_Tenure (now in the bottom-left, same width as top plots)
plt.subplot(grid[1, 0]) # Bottom-left subplot
sns.kdeplot(df_churn_.Tenure.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_val_Tenure.Tenure, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Tenure Before and After Missing Value Treatment')
plt.xlabel('Tenure')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Iterative Imputation
# Plot 4: df_val_Tenure (now in the bottom-left, same width as top plots)
plt.subplot(grid[1, 1]) # Bottom-left subplot
sns.kdeplot(df_train.Tenure.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_train_II.Tenure, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Before and After Missing Value Iterative Imputation')
plt.xlabel('Tenure')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 5: df_val_Tenure (now in the bottom-left, same width as top plots)
plt.subplot(grid[2, 0]) # Bottom-left subplot
sns.kdeplot(df_test.Tenure.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_test_II.Tenure, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_test Before and After Missing Value Iterative Imputation')
plt.xlabel('Tenure')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)


# Plot 5: df_val_Tenure (now in the bottom-left, same width as top plots)
plt.subplot(grid[2, 1]) # Bottom-left subplot
sns.kdeplot(df_val.Tenure.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_val_II.Tenure, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Before and After Missing Value Iterative Imputation')
plt.xlabel('Tenure')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

plt.tight_layout() # Adjust layout to prevent overlapping titles/labels
plt.show()

"""HourSpendOnApp"""

# Extract the 'HourSpendOnApp' column into new, separate DataFrames
df_train_HourSpendOnApp = df_train[['HourSpendOnApp']]
df_test_HourSpendOnApp = df_test[['HourSpendOnApp']]
df_val_HourSpendOnApp = df_val[['HourSpendOnApp']]

# The original df_train, df_test, and df_val remain untouched

# Initialize the imputer once
imp = SimpleImputer(missing_values=np.nan, strategy= 'most_frequent')

# Impute missing values in the training set
df_train_HourSpendOnApp['HourSpendOnApp'] = imp.fit_transform(df_train_HourSpendOnApp[['HourSpendOnApp']])

# Impute missing values in the test set using the SAME fitted imputer
df_test_HourSpendOnApp['HourSpendOnApp'] = imp.transform(df_test_HourSpendOnApp[['HourSpendOnApp']])

# Impute missing values in the validation set using the SAME fitted imputer
df_val_HourSpendOnApp['HourSpendOnApp'] = imp.transform(df_val_HourSpendOnApp[['HourSpendOnApp']])

# Display samples of the imputed columns
print("df_train_HourSpendOnApp after imputation:")
display(df_train_HourSpendOnApp['HourSpendOnApp'].head())
print("\ndf_test_HourSpendOnApp after imputation:")
display(df_test_HourSpendOnApp['HourSpendOnApp'].head())
print("\ndf_val_HourSpendOnApp after imputation:")
display(df_val_HourSpendOnApp['HourSpendOnApp'].head())

plt.figure(figsize=(15, 10)) # Adjust overall figure size as needed
grid = gridspec.GridSpec(3, 2)

# Plot 1: df_train_HourSpendOnApp
plt.subplot(grid[0, 0]) # Top-left subplot
sns.kdeplot(df_churn_.HourSpendOnApp.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_train_HourSpendOnApp.HourSpendOnApp, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_train HourSpendOnApp Before and After Missing Value Treatment')
plt.xlabel('HourSpendOnApp')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 2: df_test_HourSpendOnApp
plt.subplot(grid[0, 1]) # Top-right subplot
sns.kdeplot(df_churn_.HourSpendOnApp.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_test_HourSpendOnApp.HourSpendOnApp, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_test HourSpendOnApp Before and After Missing Value Treatment')
plt.xlabel('HourSpendOnApp')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 3: df_val_HourSpendOnApp (now in the bottom-left, same width as top plots)
plt.subplot(grid[1, 0]) # Bottom-left subplot
sns.kdeplot(df_churn_.HourSpendOnApp.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_val_HourSpendOnApp.HourSpendOnApp, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val HourSpendOnApp Before and After Missing Value Treatment')
plt.xlabel('HourSpendOnApp')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Iterative Imputation
# Plot 4: df_val_HourSpendOnApp (now in the bottom-left, same width as top plots)
plt.subplot(grid[1, 1]) # Bottom-left subplot
sns.kdeplot(df_train.HourSpendOnApp.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_train_II.HourSpendOnApp, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Before and After Missing Value Iterative Imputation')
plt.xlabel('HourSpendOnApp')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 5: df_val_HourSpendOnApp (now in the bottom-left, same width as top plots)
plt.subplot(grid[2, 0]) # Bottom-left subplot
sns.kdeplot(df_test.HourSpendOnApp.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_test_II.HourSpendOnApp, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_test Before and After Missing Value Iterative Imputation')
plt.xlabel('HourSpendOnApp')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)


# Plot 5: df_val_HourSpendOnApp (now in the bottom-left, same width as top plots)
plt.subplot(grid[2, 1]) # Bottom-left subplot
sns.kdeplot(df_val.HourSpendOnApp.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_val_II.HourSpendOnApp, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Before and After Missing Value Iterative Imputation')
plt.xlabel('HourSpendOnApp')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

plt.tight_layout() # Adjust layout to prevent overlapping titles/labels
plt.show()

"""NumberOfDeviceRegistered"""

# Extract the 'NumberOfDeviceRegistered' column into new, separate DataFrames
df_train_NumberOfDeviceRegistered = df_train[['NumberOfDeviceRegistered']]
df_test_NumberOfDeviceRegistered = df_test[['NumberOfDeviceRegistered']]
df_val_NumberOfDeviceRegistered = df_val[['NumberOfDeviceRegistered']]

# The original df_train, df_test, and df_val remain untouched

# Initialize the imputer once
imp = SimpleImputer(missing_values=np.nan, strategy= 'most_frequent')

# Impute missing values in the training set
df_train_NumberOfDeviceRegistered['NumberOfDeviceRegistered'] = imp.fit_transform(df_train_NumberOfDeviceRegistered[['NumberOfDeviceRegistered']])

# Impute missing values in the test set using the SAME fitted imputer
df_test_NumberOfDeviceRegistered['NumberOfDeviceRegistered'] = imp.transform(df_test_NumberOfDeviceRegistered[['NumberOfDeviceRegistered']])

# Impute missing values in the validation set using the SAME fitted imputer
df_val_NumberOfDeviceRegistered['NumberOfDeviceRegistered'] = imp.transform(df_val_NumberOfDeviceRegistered[['NumberOfDeviceRegistered']])

# Display samples of the imputed columns
print("df_train_NumberOfDeviceRegistered after imputation:")
display(df_train_NumberOfDeviceRegistered['NumberOfDeviceRegistered'].head())
print("\ndf_test_NumberOfDeviceRegistered after imputation:")
display(df_test_NumberOfDeviceRegistered['NumberOfDeviceRegistered'].head())
print("\ndf_val_NumberOfDeviceRegistered after imputation:")
display(df_val_NumberOfDeviceRegistered['NumberOfDeviceRegistered'].head())

plt.figure(figsize=(15, 10)) # Adjust overall figure size as needed
grid = gridspec.GridSpec(3, 2)

# Plot 1: df_train_NumberOfDeviceRegistered
plt.subplot(grid[0, 0]) # Top-left subplot
sns.kdeplot(df_churn_.NumberOfDeviceRegistered.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_train_NumberOfDeviceRegistered.NumberOfDeviceRegistered, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_train NumberOfDeviceRegistered Before and After Missing Value Treatment')
plt.xlabel('NumberOfDeviceRegistered')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 2: df_test_NumberOfDeviceRegistered
plt.subplot(grid[0, 1]) # Top-right subplot
sns.kdeplot(df_churn_.NumberOfDeviceRegistered.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_test_NumberOfDeviceRegistered.NumberOfDeviceRegistered, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_test NumberOfDeviceRegistered Before and After Missing Value Treatment')
plt.xlabel('NumberOfDeviceRegistered')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 3: df_val_NumberOfDeviceRegistered (now in the bottom-left, same width as top plots)
plt.subplot(grid[1, 0]) # Bottom-left subplot
sns.kdeplot(df_churn_.NumberOfDeviceRegistered.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_val_NumberOfDeviceRegistered.NumberOfDeviceRegistered, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val NumberOfDeviceRegistered Before and After Missing Value Treatment')
plt.xlabel('NumberOfDeviceRegistered')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Iterative Imputation
# Plot 4: df_val_NumberOfDeviceRegistered (now in the bottom-left, same width as top plots)
plt.subplot(grid[1, 1]) # Bottom-left subplot
sns.kdeplot(df_train.NumberOfDeviceRegistered.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_train_II.NumberOfDeviceRegistered, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Before and After Missing Value Iterative Imputation')
plt.xlabel('NumberOfDeviceRegistered')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 5: df_val_NumberOfDeviceRegistered (now in the bottom-left, same width as top plots)
plt.subplot(grid[2, 0]) # Bottom-left subplot
sns.kdeplot(df_test.NumberOfDeviceRegistered.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_test_II.NumberOfDeviceRegistered, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_test Before and After Missing Value Iterative Imputation')
plt.xlabel('NumberOfDeviceRegistered')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)


# Plot 5: df_val_NumberOfDeviceRegistered (now in the bottom-left, same width as top plots)
plt.subplot(grid[2, 1]) # Bottom-left subplot
sns.kdeplot(df_val.NumberOfDeviceRegistered.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_val_II.NumberOfDeviceRegistered, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Before and After Missing Value Iterative Imputation')
plt.xlabel('NumberOfDeviceRegistered')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

plt.tight_layout() # Adjust layout to prevent overlapping titles/labels
plt.show()

"""OrderAmountHikeFromlastYear"""

# Extract the 'OrderAmountHikeFromlastYear' column into new, separate DataFrames
df_train_OrderAmountHikeFromlastYear = df_train[['OrderAmountHikeFromlastYear']]
df_test_OrderAmountHikeFromlastYear = df_test[['OrderAmountHikeFromlastYear']]
df_val_OrderAmountHikeFromlastYear = df_val[['OrderAmountHikeFromlastYear']]

# The original df_train, df_test, and df_val remain untouched

# Initialize the imputer once
imp = SimpleImputer(missing_values=np.nan, strategy= 'most_frequent')

# Impute missing values in the training set
df_train_OrderAmountHikeFromlastYear['OrderAmountHikeFromlastYear'] = imp.fit_transform(df_train_OrderAmountHikeFromlastYear[['OrderAmountHikeFromlastYear']])

# Impute missing values in the test set using the SAME fitted imputer
df_test_OrderAmountHikeFromlastYear['OrderAmountHikeFromlastYear'] = imp.transform(df_test_OrderAmountHikeFromlastYear[['OrderAmountHikeFromlastYear']])

# Impute missing values in the validation set using the SAME fitted imputer
df_val_OrderAmountHikeFromlastYear['OrderAmountHikeFromlastYear'] = imp.transform(df_val_OrderAmountHikeFromlastYear[['OrderAmountHikeFromlastYear']])

# Display samples of the imputed columns
print("df_train_OrderAmountHikeFromlastYear after imputation:")
display(df_train_OrderAmountHikeFromlastYear['OrderAmountHikeFromlastYear'].head())
print("\ndf_test_OrderAmountHikeFromlastYear after imputation:")
display(df_test_OrderAmountHikeFromlastYear['OrderAmountHikeFromlastYear'].head())
print("\ndf_val_OrderAmountHikeFromlastYear after imputation:")
display(df_val_OrderAmountHikeFromlastYear['OrderAmountHikeFromlastYear'].head())

plt.figure(figsize=(15, 10)) # Adjust overall figure size as needed
grid = gridspec.GridSpec(3, 2)

# Plot 1: df_train_OrderAmountHikeFromlastYear
plt.subplot(grid[0, 0]) # Top-left subplot
sns.kdeplot(df_churn_.OrderAmountHikeFromlastYear.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_train_OrderAmountHikeFromlastYear.OrderAmountHikeFromlastYear, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_train OrderAmountHikeFromlastYear Before and After Missing Value Treatment')
plt.xlabel('OrderAmountHikeFromlastYear')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 2: df_test_OrderAmountHikeFromlastYear
plt.subplot(grid[0, 1]) # Top-right subplot
sns.kdeplot(df_churn_.OrderAmountHikeFromlastYear.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_test_OrderAmountHikeFromlastYear.OrderAmountHikeFromlastYear, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_test OrderAmountHikeFromlastYear Before and After Missing Value Treatment')
plt.xlabel('OrderAmountHikeFromlastYear')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 3: df_val_OrderAmountHikeFromlastYear (now in the bottom-left, same width as top plots)
plt.subplot(grid[1, 0]) # Bottom-left subplot
sns.kdeplot(df_churn_.OrderAmountHikeFromlastYear.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_val_OrderAmountHikeFromlastYear.OrderAmountHikeFromlastYear, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val OrderAmountHikeFromlastYear Before and After Missing Value Treatment')
plt.xlabel('OrderAmountHikeFromlastYear')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Iterative Imputation
# Plot 4: df_val_OrderAmountHikeFromlastYear (now in the bottom-left, same width as top plots)
plt.subplot(grid[1, 1]) # Bottom-left subplot
sns.kdeplot(df_train.OrderAmountHikeFromlastYear.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_train_II.OrderAmountHikeFromlastYear, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Before and After Missing Value Iterative Imputation')
plt.xlabel('OrderAmountHikeFromlastYear')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 5: df_val_OrderAmountHikeFromlastYear (now in the bottom-left, same width as top plots)
plt.subplot(grid[2, 0]) # Bottom-left subplot
sns.kdeplot(df_test.OrderAmountHikeFromlastYear.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_test_II.OrderAmountHikeFromlastYear, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_test Before and After Missing Value Iterative Imputation')
plt.xlabel('OrderAmountHikeFromlastYear')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)


# Plot 5: df_val_OrderAmountHikeFromlastYear (now in the bottom-left, same width as top plots)
plt.subplot(grid[2, 1]) # Bottom-left subplot
sns.kdeplot(df_val.OrderAmountHikeFromlastYear.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_val_II.OrderAmountHikeFromlastYear, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Before and After Missing Value Iterative Imputation')
plt.xlabel('OrderAmountHikeFromlastYear')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

plt.tight_layout() # Adjust layout to prevent overlapping titles/labels
plt.show()

"""CouponUsed"""

# Extract the 'CouponUsed' column into new, separate DataFrames
df_train_CouponUsed = df_train[['CouponUsed']]
df_test_CouponUsed = df_test[['CouponUsed']]
df_val_CouponUsed = df_val[['CouponUsed']]

# The original df_train, df_test, and df_val remain untouched

# Initialize the imputer once
imp = SimpleImputer(missing_values=np.nan, strategy= 'median')

# Impute missing values in the training set
df_train_CouponUsed['CouponUsed'] = imp.fit_transform(df_train_CouponUsed[['CouponUsed']])

# Impute missing values in the test set using the SAME fitted imputer
df_test_CouponUsed['CouponUsed'] = imp.transform(df_test_CouponUsed[['CouponUsed']])

# Impute missing values in the validation set using the SAME fitted imputer
df_val_CouponUsed['CouponUsed'] = imp.transform(df_val_CouponUsed[['CouponUsed']])

# Display samples of the imputed columns
print("df_train_CouponUsed after imputation:")
display(df_train_CouponUsed['CouponUsed'].head())
print("\ndf_test_CouponUsed after imputation:")
display(df_test_CouponUsed['CouponUsed'].head())
print("\ndf_val_CouponUsed after imputation:")
display(df_val_CouponUsed['CouponUsed'].head())

plt.figure(figsize=(15, 10)) # Adjust overall figure size as needed
grid = gridspec.GridSpec(3, 2) # Create a 2x2 grid

# Plot 1: df_train_CouponUsed
plt.subplot(grid[0, 0]) # Top-left subplot
sns.kdeplot(df_churn_.CouponUsed.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_train_CouponUsed.CouponUsed, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_train CouponUsed Before and After Missing Value Treatment')
plt.xlabel('CouponUsed')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 2: df_test_CouponUsed
plt.subplot(grid[0, 1]) # Top-right subplot
sns.kdeplot(df_churn_.CouponUsed.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_test_CouponUsed.CouponUsed, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_test CouponUsed Before and After Missing Value Treatment')
plt.xlabel('CouponUsed')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 3: df_val_CouponUsed (now in the bottom-left, same width as top plots)
plt.subplot(grid[1, 0]) # Bottom-left subplot
sns.kdeplot(df_churn_.CouponUsed.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_val_CouponUsed.CouponUsed, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val CouponUsed Before and After Missing Value Treatment')
plt.xlabel('CouponUsed')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Iterative Imputation
# Plot 4: df_val_CouponUsed (now in the bottom-left, same width as top plots)
plt.subplot(grid[1, 1]) # Bottom-left subplot
sns.kdeplot(df_train.CouponUsed.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_train_II.CouponUsed, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Before and After Missing Value Iterative Imputation')
plt.xlabel('CouponUsed')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 5: df_val_CouponUsed (now in the bottom-left, same width as top plots)
plt.subplot(grid[2, 0]) # Bottom-left subplot
sns.kdeplot(df_test.CouponUsed.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_test_II.CouponUsed, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_test Before and After Missing Value Iterative Imputation')
plt.xlabel('CouponUsed')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)


# Plot 5: df_val_CouponUsed (now in the bottom-left, same width as top plots)
plt.subplot(grid[2, 1]) # Bottom-left subplot
sns.kdeplot(df_val.CouponUsed.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_val_II.CouponUsed, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Before and After Missing Value Iterative Imputation')
plt.xlabel('CouponUsed')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

plt.tight_layout() # Adjust layout to prevent overlapping titles/labels
plt.show()

"""OrderCount"""

# Extract the 'OrderCount' column into new, separate DataFrames
df_train_OrderCount = df_train[['OrderCount']]
df_test_OrderCount = df_test[['OrderCount']]
df_val_OrderCount = df_val[['OrderCount']]

# The original df_train, df_test, and df_val remain untouched

# Initialize the imputer once
imp = SimpleImputer(missing_values=np.nan, strategy= 'median')

# Impute missing values in the training set
df_train_OrderCount['OrderCount'] = imp.fit_transform(df_train_OrderCount[['OrderCount']])

# Impute missing values in the test set using the SAME fitted imputer
df_test_OrderCount['OrderCount'] = imp.transform(df_test_OrderCount[['OrderCount']])

# Impute missing values in the validation set using the SAME fitted imputer
df_val_OrderCount['OrderCount'] = imp.transform(df_val_OrderCount[['OrderCount']])

# Display samples of the imputed columns
print("df_train_OrderCount after imputation:")
display(df_train_OrderCount['OrderCount'].head())
print("\ndf_test_OrderCount after imputation:")
display(df_test_OrderCount['OrderCount'].head())
print("\ndf_val_OrderCount after imputation:")
display(df_val_OrderCount['OrderCount'].head())

plt.figure(figsize=(15, 10)) # Adjust overall figure size as needed
grid = gridspec.GridSpec(3, 2) # Create a 2x2 grid

# Plot 1: df_train_OrderCount
plt.subplot(grid[0, 0]) # Top-left subplot
sns.kdeplot(df_churn_.OrderCount.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_train_OrderCount.OrderCount, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_train OrderCount Before and After Missing Value Treatment')
plt.xlabel('OrderCount')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 2: df_test_OrderCount
plt.subplot(grid[0, 1]) # Top-right subplot
sns.kdeplot(df_churn_.OrderCount.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_test_OrderCount.OrderCount, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_test OrderCount Before and After Missing Value Treatment')
plt.xlabel('OrderCount')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 3: df_val_OrderCount (now in the bottom-left, same width as top plots)
plt.subplot(grid[1, 0]) # Bottom-left subplot
sns.kdeplot(df_churn_.OrderCount.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_val_OrderCount.OrderCount, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val OrderCount Before and After Missing Value Treatment')
plt.xlabel('OrderCount')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Iterative Imputation
# Plot 4: df_val_OrderCount (now in the bottom-left, same width as top plots)
plt.subplot(grid[1, 1]) # Bottom-left subplot
sns.kdeplot(df_train.OrderCount.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_train_II.OrderCount, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Before and After Missing Value Iterative Imputation')
plt.xlabel('OrderCount')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 5: df_val_OrderCount (now in the bottom-left, same width as top plots)
plt.subplot(grid[2, 0]) # Bottom-left subplot
sns.kdeplot(df_test.OrderCount.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_test_II.OrderCount, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_test Before and After Missing Value Iterative Imputation')
plt.xlabel('OrderCount')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)


# Plot 5: df_val_OrderCount (now in the bottom-left, same width as top plots)
plt.subplot(grid[2, 1]) # Bottom-left subplot
sns.kdeplot(df_val.OrderCount.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_val_II.OrderCount, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Before and After Missing Value Iterative Imputation')
plt.xlabel('OrderCount')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

plt.tight_layout() # Adjust layout to prevent overlapping titles/labels
plt.show()

"""DaySinceLastOrder"""

# Extract the 'DaySinceLastOrder' column into new, separate DataFrames
df_train_DaySinceLastOrder = df_train[['DaySinceLastOrder']]
df_test_DaySinceLastOrder = df_test[['DaySinceLastOrder']]
df_val_DaySinceLastOrder = df_val[['DaySinceLastOrder']]

# The original df_train, df_test, and df_val remain untouched

# Initialize the imputer once
imp = SimpleImputer(missing_values=np.nan, strategy= 'median')

# Impute missing values in the training set
df_train_DaySinceLastOrder['DaySinceLastOrder'] = imp.fit_transform(df_train_DaySinceLastOrder[['DaySinceLastOrder']])

# Impute missing values in the test set using the SAME fitted imputer
df_test_DaySinceLastOrder['DaySinceLastOrder'] = imp.transform(df_test_DaySinceLastOrder[['DaySinceLastOrder']])

# Impute missing values in the validation set using the SAME fitted imputer
df_val_DaySinceLastOrder['DaySinceLastOrder'] = imp.transform(df_val_DaySinceLastOrder[['DaySinceLastOrder']])

# Display samples of the imputed columns
print("df_train_DaySinceLastOrder after imputation:")
display(df_train_DaySinceLastOrder['DaySinceLastOrder'].head())
print("\ndf_test_DaySinceLastOrder after imputation:")
display(df_test_DaySinceLastOrder['DaySinceLastOrder'].head())
print("\ndf_val_DaySinceLastOrder after imputation:")
display(df_val_DaySinceLastOrder['DaySinceLastOrder'].head())

plt.figure(figsize=(15, 10)) # Adjust overall figure size as needed
grid = gridspec.GridSpec(3, 2) # Create a 2x2 grid

# Plot 1: df_train_DaySinceLastOrder
plt.subplot(grid[0, 0]) # Top-left subplot
sns.kdeplot(df_churn_.DaySinceLastOrder.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_train_DaySinceLastOrder.DaySinceLastOrder, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_train DaySinceLastOrder Before and After Missing Value Treatment')
plt.xlabel('DaySinceLastOrder')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 2: df_test_DaySinceLastOrder
plt.subplot(grid[0, 1]) # Top-right subplot
sns.kdeplot(df_churn_.DaySinceLastOrder.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_test_DaySinceLastOrder.DaySinceLastOrder, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_test DaySinceLastOrder Before and After Missing Value Treatment')
plt.xlabel('DaySinceLastOrder')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 3: df_val_DaySinceLastOrder (now in the bottom-left, same width as top plots)
plt.subplot(grid[1, 0]) # Bottom-left subplot
sns.kdeplot(df_churn_.DaySinceLastOrder.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_val_DaySinceLastOrder.DaySinceLastOrder, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val DaySinceLastOrder Before and After Missing Value Treatment')
plt.xlabel('DaySinceLastOrder')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Iterative Imputation
# Plot 4: df_val_DaySinceLastOrder (now in the bottom-left, same width as top plots)
plt.subplot(grid[1, 1]) # Bottom-left subplot
sns.kdeplot(df_train.DaySinceLastOrder.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_train_II.DaySinceLastOrder, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Before and After Missing Value Iterative Imputation')
plt.xlabel('DaySinceLastOrder')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

# Plot 5: df_val_DaySinceLastOrder (now in the bottom-left, same width as top plots)
plt.subplot(grid[2, 0]) # Bottom-left subplot
sns.kdeplot(df_test.DaySinceLastOrder.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_test_II.DaySinceLastOrder, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_test Before and After Missing Value Iterative Imputation')
plt.xlabel('DaySinceLastOrder')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)


# Plot 5: df_val_DaySinceLastOrder (now in the bottom-left, same width as top plots)
plt.subplot(grid[2, 1]) # Bottom-left subplot
sns.kdeplot(df_val.DaySinceLastOrder.dropna(), label='Original (before treatment)', fill=True)
sns.kdeplot(df_val_II.DaySinceLastOrder, label='After treatment (imputed)', fill=True)
plt.title('Distribution of df_val Before and After Missing Value Iterative Imputation')
plt.xlabel('DaySinceLastOrder')
plt.ylabel('Density')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)

plt.tight_layout() # Adjust layout to prevent overlapping titles/labels
plt.show()

"""Conclustion:

Initially, the missing values were imputed using a simple imputer, but this was not good enough. There were shifts upwards and sideways. After iterative imputation was applied, the training, test and validation data fit better.

### Categorical variable encoding

As a rule of thumb, we can consider using :

 1. Label Encoding ---> Binary categorical variables and Ordinal variables
 2. One-Hot Encoding ---> Non-ordinal categorical variables with low to mid cardinality (< 5-10 levels)
 3. Target encoding ---> Categorical variables with > 10 levels

PreferredLoginDevice	PreferredPaymentMode	Gender	PreferedOrderCat	MaritalStatus

* For Gender, a simple Label encoding should be fine.
* For PreferredLoginDevice, PreferredPaymentMode, PreferedOrderCat and MaritalStatus since there are more than 2 levels, OneHotEncoding should do the trick

### Seeing Categories in a Categorical Variable
"""

df_train_II.columns

df_train_m.columns

print("Unique categories in 'PreferredLoginDevice':")
display(df_train['PreferredLoginDevice'].unique())

print("\nCategories and their counts in 'PreferredPaymentMode':")
display(df_train['PreferredPaymentMode'].unique())

print("\nCategories and their counts in 'Gender':")
display(df_train['Gender'].unique())


print("\nCategories and their counts in 'PreferedOrderCat':")
display(df_train['PreferedOrderCat'].unique())


print("\nCategories and their counts in 'MaritalStatus':")
display(df_train['MaritalStatus'].unique())

print("\nCategories and their counts in 'SatisfactionScore':")
display(df_train['SatisfactionScore'].unique())

print("\nCategories and their counts in 'Complain':")
display(df_train['Complain'].unique())

print("\nCategories and their counts in 'CityTier':")
display(df_train['CityTier'].unique())

"""#### Label Encoding for binary variables"""

df_train_m.columns

# Display unique values for numerical features
print("--- Unique values for Numerical Features ---")
for feature in num_feats:
    print(f"\nUnique values for '{feature}':")
    display(df_train_II[feature].unique())

# Display unique values and their counts for categorical features
print("\n--- Unique values and counts for Categorical Features ---")
for feature in cat_feats:
    print(f"\nUnique values for '{feature}':")
    display(df_train_II[feature].unique())
    print(f"Value counts for '{feature}':")
    display(df_train_II[feature].value_counts())

## The non-sklearn method to prevent data leakage!
# 1. Define the exact mapping (Male=1, Female=0)
gender_mapping = {'Male': 1, 'Female': 0}

# 2. Apply this fixed map to all datasets
df_train_II['Gender_cat'] = df_train_II['Gender'].map(gender_mapping)
df_val_II['Gender_cat'] = df_val_II['Gender'].map(gender_mapping)
df_test_II['Gender_cat'] = df_test_II['Gender'].map(gender_mapping)

## The non-sklearn method to prevent data leakage appliyed on df_train_m, df_test_m df_val_m!
# 1. Define the exact mapping (Male=1, Female=0)
gender_mapping = {'Male': 1, 'Female': 0}

# 2. Apply this fixed map to all datasets
df_train_m['Gender_cat'] = df_train_m['Gender'].map(gender_mapping)
df_val_m['Gender_cat'] = df_val_m['Gender'].map(gender_mapping)
df_test_m['Gender_cat'] = df_test_m['Gender'].map(gender_mapping)

"""#### One-Hot encoding for categorical variables with multiple levels"""

df_train_II.columns

# applying one hot encoding on categorical features in df_train_II, df_test_II, df_val_II

# Applying one-hot encoding on MaritalStatus
df_train_II = pd.get_dummies(df_train_II, columns=["MaritalStatus"], dtype=int)
df_test_II = pd.get_dummies(df_test_II, columns=["MaritalStatus"], dtype=int)
df_val_II = pd.get_dummies(df_val_II, columns=["MaritalStatus"], dtype=int)

# Applying one-hot encoding on PreferredPaymentMode

df_train_II = pd.get_dummies(df_train_II, columns=["PreferredPaymentMode"], dtype=int)
df_test_II = pd.get_dummies(df_test_II, columns=["PreferredPaymentMode"], dtype=int)
df_val_II = pd.get_dummies(df_val_II, columns=["PreferredPaymentMode"], dtype=int)

# Applying one-hot encoding on PreferedOrderCat

df_train_II = pd.get_dummies(df_train_II, columns=["PreferedOrderCat"], dtype=int)
df_test_II = pd.get_dummies(df_test_II, columns=["PreferedOrderCat"], dtype=int)
df_val_II = pd.get_dummies(df_val_II, columns=["PreferedOrderCat"], dtype=int)

# Applying one-hot encoding on PreferredLoginDevice

df_train_II = pd.get_dummies(df_train_II, columns=["PreferredLoginDevice"], dtype=int)
df_test_II = pd.get_dummies(df_test_II, columns=["PreferredLoginDevice"], dtype=int)
df_val_II = pd.get_dummies(df_val_II, columns=["PreferredLoginDevice"], dtype=int)

# Applying one-hot encoding on CityTier

df_train_II = pd.get_dummies(df_train_II, columns=["CityTier"], dtype=int)
df_test_II = pd.get_dummies(df_test_II, columns=["CityTier"], dtype=int)
df_val_II = pd.get_dummies(df_val_II, columns=["CityTier"], dtype=int)

# applying one hot encoding on categorical features in df_train_m, df_test_m, df_val_m

# Applying one-hot encoding on MaritalStatus
df_train_m = pd.get_dummies(df_train_m, columns=["MaritalStatus"], dummy_na=True, dtype=int)
df_test_m = pd.get_dummies(df_test_m, columns=["MaritalStatus"], dummy_na=True, dtype=int)
df_val_m = pd.get_dummies(df_val_m, columns=["MaritalStatus"], dummy_na=True, dtype=int)

# Applying one-hot encoding on PreferredPaymentMode

df_train_m = pd.get_dummies(df_train_m, columns=["PreferredPaymentMode"], dummy_na=True, dtype=int)
df_test_m = pd.get_dummies(df_test_m, columns=["PreferredPaymentMode"], dummy_na=True,dtype=int)
df_val_m = pd.get_dummies(df_val_m, columns=["PreferredPaymentMode"], dummy_na=True ,dtype=int)

# Applying one-hot encoding on PreferedOrderCat

df_train_m = pd.get_dummies(df_train_m, columns=["PreferedOrderCat"], dummy_na=True ,dtype=int)
df_test_m = pd.get_dummies(df_test_m, columns=["PreferedOrderCat"],dummy_na=True ,dtype=int)
df_val_m = pd.get_dummies(df_val_m, columns=["PreferedOrderCat"], dummy_na=True,dtype=int)

# Applying one-hot encoding on PreferredLoginDevice

df_train_m = pd.get_dummies(df_train_m, columns=["PreferredLoginDevice"], dummy_na=True,dtype=int)
df_test_m = pd.get_dummies(df_test_m, columns=["PreferredLoginDevice"], dummy_na=True,dtype=int)
df_val_m = pd.get_dummies(df_val_m, columns=["PreferredLoginDevice"], dummy_na=True,dtype=int)

# Applying one-hot encoding on CityTier

df_train_m = pd.get_dummies(df_train_m, columns=["CityTier"], dummy_na=True,dtype=int)
df_test_m = pd.get_dummies(df_test_m, columns=["CityTier"], dummy_na=True,dtype=int)
df_val_m = pd.get_dummies(df_val_m, columns=["CityTier"], dummy_na=True,dtype=int)

"""####  Removing the original feature and Unnamed: 0 and CustomerID"""

## Drop the Gender column in df_train_II, df_val_II, df_test_II
df_train_II.drop(['Gender'], axis = 1, inplace=True)
df_val_II.drop(['Gender'], axis = 1, inplace=True)
df_test_II.drop(['Gender'], axis = 1, inplace=True)

## Drop the Gender column in df_train_m, df_test_m, df_val_m
df_train_m.drop(['Gender'], axis = 1, inplace=True)
df_val_m.drop(['Gender'], axis = 1, inplace=True)
df_test_m.drop(['Gender'], axis = 1, inplace=True)

df_train_m = df_train_m.drop(columns=['Unnamed: 0', "CustomerID"])
df_test_m = df_test_m.drop(columns=['Unnamed: 0', "CustomerID"])
df_val_m = df_val_m.drop(columns=['Unnamed: 0', "CustomerID"])

"""#### _Summarize_ : How to handle unknown categorical levels/values in unseen data in production?

- Use LabelEncoding, OneHotEncoding on training set and then save the mapping and apply on the test set. For missing values, use 0, -1 etc.

### Bivariate analysis
"""

## Check linear correlation (rho) between individual features and the target variable
## Imputred training set
corr_train_II = df_train_II.corr()
corr_train_II.shape

sns.heatmap(corr_train_II, cmap = 'coolwarm')

# If not, you'd calculate it from your main data like this:
# corr_train_II = df.corr()

# Define the target feature name (adjust if your target name is different)
target_feature = 'Churn'

# Assuming the correlation matrix includes the target_feature as a column/index.
# Get all feature names except the target
all_features = [col for col in corr_train_II.columns if col != target_feature]
num_features = len(all_features) # Should be 38 based on your (39, 39) shape

# Randomly select 9 features for each plot to reach 10 features total (9 + target)
# The specific selection logic below is an example; adjust as needed for your specific feature groups.

# Ensure we have at least 27 other features to select 3 distinct groups of 9
if num_features >= 27:
    # Example: Selecting the first 9, next 9, and next 9 features for 3 plots
    features_1 = all_features[:9]
    features_2 = all_features[9:18]
    features_3 = all_features[18:27]
    # For the 4th plot, include the remaining 11 features
    features_4 = all_features[27:]
else:
    # Handle cases with fewer features if necessary (e.g., sample features randomly with replacement)
    # The example code assumes > 27 features are available.
    features_1 = np.random.choice(all_features, 9, replace=False).tolist()
    features_2 = np.random.choice([f for f in all_features if f not in features_1], 9, replace=False).tolist()
    features_3 = np.random.choice([f for f in all_features if f not in features_1 + features_2], 9, replace=False).tolist()
    features_4 = np.random.choice([f for f in all_features if f not in features_1 + features_2 + features_3], 9, replace=False).tolist()


# Combine each feature list with the target feature
plot_features_list = [
    features_1 + [target_feature],
    features_2 + [target_feature],
    features_3 + [target_feature],
    features_4 + [target_feature]
]

# Create a 2x2 grid of subplots
fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(16, 14))
# Flatten the axes array for easier iteration
axs = axs.flatten()

# Iterate through the feature lists and plot the heatmaps
for i, features_to_plot in enumerate(plot_features_list):
    # Select the subset from the main correlation dataframe
    subset_corr = corr_train_II.loc[features_to_plot, features_to_plot]

    # Plot the heatmap on the corresponding axis
    sns.heatmap(
        subset_corr,
        annot=False,
        cmap='coolwarm',
        fmt=".2f",
        linewidths=0.5,
        ax=axs[i],
        square=True
    )

    # Set titles for each subplot
    axs[i].set_title(f'Correlation Matrix Subset {i+1} (N={len(features_to_plot)})')

# Adjust layout to prevent overlap and display the plot
plt.tight_layout()
plt.show()

## Check linear correlation (rho) between individual features and the target variable
## Non Imputed training set
corr_train_m = df_train_m.corr()
corr_train_m

sns.heatmap(corr_train_m, cmap = 'coolwarm')

# Define the target feature name (adjust if your target name is different)
target_feature = 'Churn'

# Assuming the correlation matrix includes the target_feature as a column/index.
# Get all feature names except the target
all_features = [col for col in corr_train_m.columns if col != target_feature]
num_features = len(all_features) # Should be 38 based on your (39, 39) shape

# Randomly select 9 features for each plot to reach 10 features total (9 + target)
# The specific selection logic below is an example; adjust as needed for your specific feature groups.

# Ensure we have at least 27 other features to select 3 distinct groups of 9
if num_features >= 27:
    # Example: Selecting the first 9, next 9, and next 9 features for 3 plots
    features_1 = all_features[:9]
    features_2 = all_features[9:18]
    features_3 = all_features[18:27]
    # For the 4th plot, include the remaining 11 features
    features_4 = all_features[27:]
else:
    # Handle cases with fewer features if necessary (e.g., sample features randomly with replacement)
    # The example code assumes > 27 features are available.
    features_1 = np.random.choice(all_features, 9, replace=False).tolist()
    features_2 = np.random.choice([f for f in all_features if f not in features_1], 9, replace=False).tolist()
    features_3 = np.random.choice([f for f in all_features if f not in features_1 + features_2], 9, replace=False).tolist()
    features_4 = np.random.choice([f for f in all_features if f not in features_1 + features_2 + features_3], 9, replace=False).tolist()


# Combine each feature list with the target feature
plot_features_list = [
    features_1 + [target_feature],
    features_2 + [target_feature],
    features_3 + [target_feature],
    features_4 + [target_feature]
]

# Create a 2x2 grid of subplots
fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(16, 14))
# Flatten the axes array for easier iteration
axs = axs.flatten()

# Iterate through the feature lists and plot the heatmaps
for i, features_to_plot in enumerate(plot_features_list):
    # Select the subset from the main correlation dataframe
    subset_corr = corr_train_m.loc[features_to_plot, features_to_plot]

    # Plot the heatmap on the corresponding axis
    sns.heatmap(
        subset_corr,
        annot=False,
        cmap='coolwarm',
        fmt=".2f",
        linewidths=0.5,
        ax=axs[i],
        square=True
    )

    # Set titles for each subplot
    axs[i].set_title(f'Correlation Matrix Subset {i+1} (N={len(features_to_plot)})')

# Adjust layout to prevent overlap and display the plot
plt.tight_layout()
plt.show()

from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.tools.tools import add_constant

# compute the vif for all given features
def compute_vif(considered_features):

    X = df_train_II[considered_features]
    # the calculation of variance inflation requires a constant
    X['intercept'] = 1

    # create dataframe to store vif values
    vif = pd.DataFrame()
    vif["Variable"] = X.columns
    vif["VIF"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
    vif = vif[vif['Variable']!='intercept']
    return vif

# show the unique values in df_train_II

print(df_train_II.columns)

# Assuming your DataFrame is named 'df_train_II'
features_to_exclude = ['Churn', "MaritalStatus_Unknown", "PreferredPaymentMode_Unknown", "PreferedOrderCat_Others",
                       "PreferredLoginDevice_Unknown", "CityTier_2", "PreferredLoginDevice_Mobile Phone",
                       "PreferedOrderCat_Laptop & Accessory", "PreferredPaymentMode_Debit Card", "MaritalStatus_Married",
                       "CityTier_3"]

# Get all columns and filter out the features that are IN the exclude list
considered_features = [col for col in df_train_II.columns if col not in features_to_exclude]

# Print the result
print(considered_features)

compute_vif(considered_features).sort_values('VIF', ascending=False)

# removing the selected featrues from df_train, df_val and df_test

features_to_exclude = ["MaritalStatus_Unknown", "PreferredPaymentMode_Unknown", "PreferedOrderCat_Others",
                       "PreferredLoginDevice_Unknown", "CityTier_2", "PreferredLoginDevice_Mobile Phone",
                       "PreferedOrderCat_Laptop & Accessory", "PreferredPaymentMode_Debit Card", "MaritalStatus_Married",
                       "CityTier_3"]

df_train_II = df_train_II.drop(columns=features_to_exclude)
df_val_II = df_val_II.drop(columns=features_to_exclude)
df_test_II = df_test_II.drop(columns=features_to_exclude)

print(df_train_II.shape)
print(df_val_II.shape)
print(df_test_II.shape)

"""None of the features are highly correlated with the target variable. But some of them have slight linear associations with the target variable.

* Most of the features have a linear relationship with the target feature Churn

* There is a strong correlation among the independent variables. First, one of the categorical varialbes have been removed to avoid dummy trap (VIF = 1/1-R^2), then based on VIF value iteratively features have been selected and removed one by one.

### Some basic feature engineering NOTHING DONE YET

### Feature scaling and normalization

Different methods :

1. Feature transformations - Using log, log10, sqrt, pow
2. MinMaxScaler - Brings all feature values between 0 and 1
3. StandardScaler - Mean normalization. Feature values are an estimate of their z-score


* Why is scaling and normalization required ?


* How do we normalize unseen data?

#### Feature transformations
"""

### Demo-ing feature transformations
sns.distplot(df_train_II.Tenure, hist=False)

sns.distplot(np.log(df_train_II.Tenure), hist=False)

"""#### StandardScaler"""

from sklearn.preprocessing import StandardScaler
sc = StandardScaler()

df_train_II.columns

"""Scaling only continuous variables"""

cont_vars = ["Tenure", "WarehouseToHome", "HourSpendOnApp", "NumberOfAddress","OrderAmountHikeFromlastYear",
             "CouponUsed", "OrderCount", "DaySinceLastOrder", "CashbackAmount"]

cat_vars = ["NumberOfDeviceRegistered" ,"SatisfactionScore", "Complain", "Gender_cat", "MaritalStatus_Divorced", "MaritalStatus_Single",
            "PreferredPaymentMode_CC", "PreferredPaymentMode_COD", "PreferredPaymentMode_Cash on Delivery",
            'PreferredPaymentMode_Credit Card', 'PreferredPaymentMode_E wallet',
            'PreferredPaymentMode_UPI', 'PreferedOrderCat_Fashion',
            'PreferedOrderCat_Grocery', 'PreferedOrderCat_Mobile',
            'PreferedOrderCat_Mobile Phone', 'PreferredLoginDevice_Computer',
            'PreferredLoginDevice_Phone', 'CityTier_1']

## Scaling only continuous columns
cols_to_scale = cont_vars

sc_X_train_II = sc.fit_transform(df_train_II[cols_to_scale])

## Converting from array to dataframe and naming the respective features/columns
sc_X_train_II = pd.DataFrame(data = sc_X_train_II, columns = cols_to_scale)
sc_X_train_II.shape
sc_X_train_II.head()

df_train_II.shape

## Mapping learnt on the continuous features
sc_map = {'mean':sc.mean_, 'std':np.sqrt(sc.var_)}
sc_map

## Scaling validation and test sets by transforming the mapping obtained through the training set
sc_X_val_II = sc.transform(df_val_II[cols_to_scale])
sc_X_test_II = sc.transform(df_test_II[cols_to_scale])

## Converting val and test arrays to dataframes for re-usability
sc_X_val_II = pd.DataFrame(data = sc_X_val_II, columns = cols_to_scale)
sc_X_test_II = pd.DataFrame(data = sc_X_test_II, columns = cols_to_scale)

print(sc_X_test_II.shape)
print(df_test_II.shape)
print(sc_X_val_II.shape)
print(df_val_II.shape)

# 1. Identify your categorical feature names (you need this list)
cat_vars = ["NumberOfDeviceRegistered","SatisfactionScore", "Complain", "Gender_cat", "MaritalStatus_Divorced", "MaritalStatus_Single",
            "PreferredPaymentMode_Cash on Delivery",
            'PreferredPaymentMode_Credit Card', 'PreferredPaymentMode_E wallet',
            'PreferredPaymentMode_UPI', 'PreferedOrderCat_Fashion',
            'PreferedOrderCat_Grocery', 'PreferedOrderCat_Mobile',
            'PreferedOrderCat_Mobile Phone', 'PreferredLoginDevice_Computer',
            'PreferredLoginDevice_Phone', 'CityTier_1']

# 0. Define your variable types, ensuring 'Churn' is NOT in the feature lists

# 'Churn' is the name of your target variable
target_var = 'Churn'
# 'num_vars' and 'cat_vars' only contain predictor columns (X)

# 1. Separate Features (X) and Target (y) *first* from your raw data
X_train_II = df_train_II.drop(columns=[target_var])
y_train_II = df_train_II[target_var]

# Now, both X_train_II and y_train_II have the consistent shape (4458, ...)

# 2. Extract only those specific columns into a temporary dataframe
#    We use X_train_II (features only) as the source, not df_train_II
X_train_cat_df = X_train_II[cat_vars]

# 3. Concatenate the scaled numerical DF and the encoded categorical DF
#    You need a new variable name for the final combined dataset.

# --- Use the robust index reset fix from before ---
sc_X_train_II_reset = sc_X_train_II.reset_index(drop=True)
X_train_cat_df_reset = X_train_cat_df.reset_index(drop=True)

sc_en_X_train_II = pd.concat([sc_X_train_II_reset, X_train_cat_df_reset], axis=1)


# Verify the shape and columns
print("Shape of combined features data (X):", sc_en_X_train_II.shape)
print("Shape of target data (y):", y_train_II.shape)

# Use sc_en_X_train_II for your outlier analysis/KNN modeling

# 0. Define your variable types, ensuring 'Churn' is NOT in the feature lists

#  'Churn' is the name of your target variable
target_var = 'Churn'
# Assuming 'num_vars' and 'cat_vars' only contain predictor columns (X)

# 1. Separate Features (X) and Target (y) *first* from your raw data
X_test_II = df_test_II.drop(columns=[target_var])
y_test_II = df_test_II[target_var]

# Now, both X_test_II and y_test_II have the consistent shape (4458, ...)

# 2. Extract only those specific columns into a temporary dataframe
#    We use X_test_II (features only) as the source, not df_test_II
X_test_cat_df = X_test_II[cat_vars]

# 3. Concatenate the scaled numerical DF and the encoded categorical DF
#    You need a new variable name for the final combined dataset.

# --- Use the robust index reset fix from before ---
sc_X_test_II_reset = sc_X_test_II.reset_index(drop=True)
X_test_cat_df_reset = X_test_cat_df.reset_index(drop=True)

sc_en_X_test_II = pd.concat([sc_X_test_II_reset, X_test_cat_df_reset], axis=1)


# Verify the shape and columns
print("Shape of combined features data (X):", sc_en_X_test_II.shape)
print("Shape of target data (y):", y_test_II.shape)

# Use sc_en_X_test_II for your outlier analysis/KNN modeling

# 0. Define your variable types, ensuring 'Churn' is NOT in the feature lists

# 'Churn' is the name of your target variable
target_var = 'Churn'
# Assuming 'num_vars' and 'cat_vars' only contain predictor columns (X)

# 1. Separate Features (X) and Target (y) *first* from your raw data
X_val_II = df_val_II.drop(columns=[target_var])
y_val_II = df_val_II[target_var]

# Now, both X_val_II and y_val_II have the consistent shape (4458, ...)

# 2. Extract only those specific columns into a temporary dataframe
#    We use X_val_II (features only) as the source, not df_val_II
X_val_cat_df = X_val_II[cat_vars]

# 3. Concatenate the scaled numerical DF and the encoded categorical DF
#    You need a new variable name for the final combined dataset.

# --- Use the robust index reset fix from before ---
sc_X_val_II_reset = sc_X_val_II.reset_index(drop=True)
X_val_cat_df_reset = X_val_cat_df.reset_index(drop=True)

sc_en_X_val_II = pd.concat([sc_X_val_II_reset, X_val_cat_df_reset], axis=1)


# Verify the shape and columns
print("Shape of combined features data (X):", sc_en_X_val_II.shape)
print("Shape of target data (y):", y_val_II.shape)

# Use sc_en_X_val_II for your outlier analysis/KNN modeling

# checking the shape of the dataframes after preprocess. Target Churn is exluced from scaled dataframes
sc_en_X_val_II.shape
df_val_II.shape
sc_en_X_train_II.shape
df_train_II.shape
sc_en_X_test_II.shape
df_test_II.shape

"""#### Local Outlier Factor"""

# Redefine the classifier with novelty=True
clf_novelty = LocalOutlierFactor(
    n_neighbors= 19, #53
    contamination="auto",
    algorithm='auto',
    metric='minkowski',
    p=2,
    leaf_size = 35,
    novelty=True, # <-- CHANGE THIS
    n_jobs=None,
    metric_params=None
)

# Fit ONLY on the training data
clf_novelty.fit(sc_en_X_train_II)

# Get the scores for training data (like X_scores)
train_scores = clf_novelty.decision_function(sc_en_X_train_II)

# apply the exact same model instance to score the test/validation data SAFELY
# You must have already preprocessed (scaled/encoded) your test/val data identically
val_scores = clf_novelty.decision_function(sc_en_X_val_II)
test_scores = clf_novelty.decision_function(sc_en_X_test_II) # <- This prevents data leakage

# Add the scores and prediction labels as new columns to the original DataFrame

# FIX 1: Change X_scores to train_scores
df_train_II['outlier_score'] = train_scores

# You may also need to define 'y_pred' if you haven't yet,
# likely by using clf_novelty.predict(sc_en_X_train_II)

print(f"Combined DataFrame shape: {df_train_II.shape}")

# 1. Ensure this line runs successfully in your environment:
y_pred = clf_novelty.predict(sc_en_X_train_II)

# 2. Then run these lines:
df_train_II['outlier_score'] = train_scores
df_train_II['is_outlier_pred'] = y_pred

# 3. Print or display the head of your DataFrame to confirm both columns exist:
(df_train_II.head())

# Generate scores and predictions for the VALIDATION set
val_scores = clf_novelty.decision_function(sc_en_X_val_II)
val_predictions = clf_novelty.predict(sc_en_X_val_II)

# Generate scores and predictions for the TEST set
test_scores = clf_novelty.decision_function(sc_en_X_test_II)
test_predictions = clf_novelty.predict(sc_en_X_test_II)

# Apply to Validation DF
df_val_II['outlier_score'] = val_scores
df_val_II['is_outlier_pred'] = val_predictions
print(f"Combined Validation DataFrame shape: {df_val_II.shape}")

# Apply to Test DF
df_test_II['outlier_score'] = test_scores
df_test_II['is_outlier_pred'] = test_predictions
print(f"Combined Test DataFrame shape: {df_test_II.shape}")

# Assuming 'sc_en_X_train_II' is a pandas DataFrame or numpy array with only 2 dimensions (2 features)

def update_legend_marker_size(handle, orig):
    "Customize size of the legend marker"
    handle.update_from(orig)
    handle.set_sizes([20])

# --- Adjustments Needed ---

# 1. Convert DataFrame to NumPy array if it isn't already, so slicing ([:, 0]) works easily
if hasattr(sc_en_X_train_II, 'values'):
    X_plot = sc_en_X_train_II.values
else:
    X_plot = sc_en_X_train_II

# 2. Ensure train_scores is available from your previous code execution
#    Make sure the train_scores variable is populated from the clf.negative_outlier_factor_ attribute

# 3. Handle the removed 'n_errors' variable

# --- Plotting ---

plt.figure(figsize=(8, 6))

plt.scatter(X_plot[:, 0], X_plot[:, 1], color="k", s=3.0, label="Data points")

# plot circles with radius proportional to the outlier scores
# Radius calculation remains the same
radius = (train_scores.max() - train_scores) / (train_scores.max() - train_scores.min())
scatter = plt.scatter(
    X_plot[:, 0], # Use X_plot here
    X_plot[:, 1], # Use X_plot here
    s=1000 * radius,
    edgecolors="r",
    facecolors="none",
    label="Outlier scores",
)
plt.axis("tight")
# The limits below might need adjustment based on the scale of your *standardized* data
# If your data is standardized, limits around -3 to 3 are often reasonable
plt.xlim((X_plot[:, 0].min() - 0.5, X_plot[:, 0].max() + 0.5))
plt.ylim((X_plot[:, 1].min() - 0.5, X_plot[:, 1].max() + 0.5))

# Remove the reference to n_errors since it doesn't exist without ground truth
plt.xlabel("Data points are colored black, outliers circled in red.")

plt.legend(
    handler_map={scatter: HandlerPathCollection(update_func=update_legend_marker_size)}
)
plt.title("Local Outlier Factor (LOF) Visualization")
plt.show()

# Define the parameters we want to experiment with
k_values_to_test = [1, 20, 35 ,53, 100]
contamination_value = 'auto' # Keep this consistent for comparison

plt.figure(figsize=(16, 8))

# Loop through each k value
for i, k in enumerate(k_values_to_test):
    # 1. Instantiate the model with the current k value
    clf = LocalOutlierFactor(
        n_neighbors=k,
        contamination=contamination_value,
        novelty=False
    )

    # 2. Fit the model and get the scores
    # We don't need y_pred here, just the scores
    clf.fit(sc_en_X_train_II)
    X_scores = clf.negative_outlier_factor_

    # 3. Visualize the distribution of scores for this K value
    ax = plt.subplot(1, len(k_values_to_test), i + 1)
    ax.hist(X_scores, bins=30, density=True, color='skyblue', edgecolor='black')
    ax.set_title(f'k = {k}')
    ax.set_xlabel('Negative Outlier Factor Score')
    if i == 0:
        ax.set_ylabel('Density of Observations')

plt.tight_layout()
plt.show()

# Create a DataFrame with all data points having a score -1
potential_anomalies_train_df_II = df_train_II[df_train_II['is_outlier_pred'] == -1].copy()

print(f"Number of potential anomalies predicted -1: {len(potential_anomalies_train_df_II)}")

(potential_anomalies_train_df_II.sort_values(by='is_outlier_pred', ascending=True).head(50))

# 1. Temporarily add the 'is_outlier_pred' column to the *processed* dataframe
# This works if the row order hasn't changed since you created sc_en_X_train_II
#sc_en_X_train_II['is_outlier_pred'] = df_train_II['is_outlier_pred'].values


# 2. Use boolean masking to filter the rows where the prediction is -1
#anomalies_in_processed_data = sc_en_X_train_II[sc_en_X_train_II['is_outlier_pred'] == -1]


# 3. View the results
#(f"Number of anomalies found in sc_en_X_train_II: {len(anomalies_in_processed_data)}")
#("\nDataFrame containing only the anomaly rows (features scaled and encoded):")
#(anomalies_in_processed_data) # Display the first few anomaly rows

#  original boolean column creation
df_train_II['is_manual_outlier'] = df_train_II['is_outlier_pred'] ==-1

# Define the custom color palette: True maps to Red, False maps to Blue
custom_palette = {True: 'red', False: 'blue'}

# --- Sampling Logic ---

# 1. Separate the outliers (True) and non-outliers (False)
df_outliers = df_train_II[df_train_II['is_manual_outlier'] == True]
df_non_outliers = df_train_II[df_train_II['is_manual_outlier'] == False]

# 2. Sample the non-outliers (e.g., take a sample of 500 points)
# Replace '500' with your desired sample size, or use 'frac=0.1' for a percentage
sample_size = 1000
df_non_outliers_sampled = df_non_outliers.sample(n=sample_size, random_state=42) # Using random_state for reproducibility

# 3. Combine the full outlier set with the sampled non-outlier set
df_sampled_plot = pd.concat([df_outliers, df_non_outliers_sampled])

# --- Plotting Code (uses the new combined DataFrame) ---
plt.figure(figsize=(10, 6))

sns.scatterplot(
    data=df_sampled_plot,       # Use the new SAMPLED DataFrame
    x='WarehouseToHome',
    y='Tenure',
    hue='is_manual_outlier',
    palette=custom_palette,
    s=50,
    alpha=0.7
)

plt.title(f"Sampled Data Points: All outliers (Red) + {sample_size} non-outliers (Blue)")
plt.xlabel("WarehouseToHome")
plt.ylabel("Tenure")
plt.legend(title="predicted outlier red")
plt.show()

valid_scores = {1, 2, 3, 4, 5}

# Replace invalid values with -1 (unknown)
df_train_II.loc[
    ~df_train_II['SatisfactionScore'].isin(valid_scores),
    'SatisfactionScore'
] = -1

# --- Setup ---
df_train_II['is_manual_outlier'] = df_train_II['is_outlier_pred'] == -1
df_outliers = df_train_II[df_train_II['is_manual_outlier'] == True].copy()
df_non_outliers = df_train_II[df_train_II['is_manual_outlier'] == False].copy()

# Sample the non-outliers (using the sample size from your code)
sample_size = 50
df_non_outliers_sampled = df_non_outliers.sample(n=sample_size, random_state=42)

# Combine the full outlier set with the sampled non-outlier set
df_sampled_plot = pd.concat([df_outliers, df_non_outliers_sampled])


# IMPORTANT: Map the values BEFORE plotting so the legend handles them correctly
df_sampled_plot['is_manual_outlier'] = df_sampled_plot['is_manual_outlier'].map({True: 'Outlier', False: 'Non-outlier'})


# Define the columns you want to include in the parallel coordinates plot.
# You must list all relevant numerical features.
# Replace the list below with the actual names of your other numerical columns.
columns_to_plot = df_train_II.columns.tolist()

# Ensure the columns you want to plot are numeric (except the target class column)
# If they are not, you may need to cast them: df_sampled_plot[col] = pd.to_numeric(df_sampled_plot[col])


# --- Parallel Coordinates Plotting ---

plt.figure(figsize=(12, 8))

# 1. Update the color list to match alphabetical order: [Non-outlier, Outlier]
colors = ['lightgrey', 'red']

parallel_coordinates(
    frame=df_sampled_plot[columns_to_plot], # Pass only the relevant columns
    class_column='is_manual_outlier',      # The column to use for colors
    cols=[col for col in columns_to_plot if col != 'is_manual_outlier'], # Explicitly define axes order
    color=colors,  # Use the defined list of colors
    sort_labels=True
)
# Automatisch generierte Legende entfernen
if ax.get_legend():
    ax.get_legend().remove()

# --- MANUAL LEGEND FIX ---
# Manuelle Linien für die Legende definieren
# 3. Ensure the manual legend handles match this order
outlier_line = mlines.Line2D([], [], color='red', label='Outlier')
non_outlier_line = mlines.Line2D([], [], color='lightgrey', label='Non-outlier')
plt.legend(handles=[outlier_line, non_outlier_line], title="Classification")

plt.title(f"Parallel Coordinates Plot (Sampled Data: All outliers + {sample_size} non-outliers)")
plt.xlabel("Features")
plt.ylabel("Features Values")
plt.xticks(rotation=90)
plt.show()

"""What is striking, based on the coordinate parallel plot, is that individuals considered outliers have a longer warehouse-to-home distance, more addresses and a higher number of coupons used and orders placed.  """

df_test_II.sample(10)

print(f"The outlier threshold (offset) is: {clf_novelty.offset_}")

# Get scores for the test outliers only
test_outlier_scores = test_scores[test_predictions == -1]

print(f"Min score: {test_outlier_scores.min()}")
print(f"Max score: {test_outlier_scores.max()}")
print(f"Number of 'borderline' outliers (between -1.5 and -1.6): {((test_outlier_scores < -1.5) & (test_outlier_scores > -1.6)).sum()}")

# Create a DataFrame with all data points having a prediction of -1

potential_anomalies_test_df_II = df_test_II[df_test_II['is_outlier_pred'] == -1].copy()

print(f"Number of potential anomalies predicted -1: {len(potential_anomalies_test_df_II)}")

potential_anomalies_test_df_II

# Create a DataFrame with all data points having a prediction of -1

potential_anomalies_val_df_II = df_val_II[df_val_II['is_outlier_pred'] == -1].copy()

print(f"Number of potential anomalies predicted -1: {len(potential_anomalies_val_df_II)}")

potential_anomalies_val_df_II

# 2. Then run these lines:
sc_en_X_train_II['outlier_score'] = train_scores
sc_en_X_train_II['is_outlier_pred'] = y_pred

# 3. Print or display the head of your DataFrame to confirm both columns exist:
(sc_en_X_train_II.head(10))
outliers_df = sc_en_X_train_II[sc_en_X_train_II['is_outlier_pred'] == -1]
outliers_df

# 2. Then run these lines:
sc_en_X_test_II['outlier_score'] = test_scores
sc_en_X_test_II['is_outlier_pred'] = test_predictions

# 3. Print or display the head of your DataFrame to confirm both columns exist:

outliers_test_df = sc_en_X_test_II[sc_en_X_test_II['is_outlier_pred'] == -1]
outliers_test_df

# 2. Then run these lines:
sc_en_X_val_II['outlier_score'] = val_scores
sc_en_X_val_II['is_outlier_pred'] = val_predictions

# 3. Print or display the head of your DataFrame to confirm both columns exist:

outliers_val_df = sc_en_X_val_II[sc_en_X_val_II['is_outlier_pred'] == -1]
outliers_val_df

# Correct Way: Use boolean filtering to select only the non-outliers REMOVING OUTLIERS FROM THE DFS
sc_en_X_val_II_NoOutliers = sc_en_X_val_II[sc_en_X_val_II['is_outlier_pred'] != -1]
sc_en_X_test_II_NoOutliers = sc_en_X_test_II[sc_en_X_test_II['is_outlier_pred'] != -1]
sc_en_X_train_II_NoOutliers = sc_en_X_train_II[sc_en_X_train_II['is_outlier_pred'] != -1]

df_train_II_NoOutliers = df_train_II[df_train_II['is_outlier_pred'] != -1]
df_test_II_NoOutliers = df_test_II[df_test_II['is_outlier_pred'] != -1]
df_val_II_NoOutliers = df_val_II[df_val_II['is_outlier_pred'] != -1]

# Assuming you have original y dataframes aligned with the df_II dataframes:

# Apply the same boolean mask to the target variables (y_train, y_val, y_test)
y_train_NoOutliers = y_train[df_train_II['is_outlier_pred'] != -1]
y_val_NoOutliers = y_val[df_val_II['is_outlier_pred'] != -1]
y_test_NoOutliers = y_test[df_test_II['is_outlier_pred'] != -1]

sc_en_X_train_II.columns

# Use the corrected variable name
columns_to_drop = ['outlier_score', 'is_outlier_pred', "is_manual_outlier"]

# Drop columns from all relevant DataFrames
df_train_II = df_train_II.drop(columns_to_drop, axis=1)
df_train_II_NoOutliers = df_train_II_NoOutliers.drop(columns_to_drop, axis=1)

# Use the corrected variable name
columns_to_drop = ['outlier_score', 'is_outlier_pred']

# Drop columns from all relevant DataFrames
df_test_II = df_test_II.drop(columns_to_drop, axis=1)
df_test_II_NoOutliers = df_test_II_NoOutliers.drop(columns_to_drop, axis=1)

df_val_II = df_val_II.drop(columns_to_drop, axis=1)
df_val_II_NoOutliers = df_val_II_NoOutliers.drop(columns_to_drop, axis=1)

"""Conclusion:
- In this section, new dfs have been generated excluding the outliers from the imputed dataset.

### Inital fits on imputed dataframe
"""

# Obtaining class weights based on the class samples imbalance ratio
_, num_samples = np.unique(y_train, return_counts = True)
weights = np.max(num_samples)/num_samples
weights
num_samples

weights_dict = dict()
class_labels = [0,1]
for a,b in zip(class_labels,weights):
    weights_dict[a] = b

weights_dict

# Assuming 'weights' is a list or array like [1.0, 5.0243...]
weights_dict = dict()
class_labels = [0,1]
for a,b in zip(class_labels,weights):
    weights_dict[a] = b

weights_dict = {0: np.float64(1.0), 1: np.float64(5.024324324324325)}

# ADD THIS LINE:
imbalance_ratio = weights_dict[1] / weights_dict[0] # Calculate the ratio of the positive class weight to the negative class weight
# imbalance_ratio will be np.float64(5.024324324324325)

## Defining model
lr = LogisticRegression(C = 1.0, penalty = 'l2', class_weight = weights_dict, n_jobs = -1)

## Fitting model
lr.fit(sc_en_X_train_II, y_train)

## Fitted model parameters
sc_en_X_train_II.columns
lr.coef_
lr.intercept_

## Training metrics
roc_auc_score(y_train, lr.predict(sc_en_X_train_II))
recall_score(y_train, lr.predict(sc_en_X_train_II))
confusion_matrix(y_train, lr.predict(sc_en_X_train_II))
print(classification_report(y_train, lr.predict(sc_en_X_train_II)))

## Validation metrics
roc_auc_score(y_val, lr.predict(sc_en_X_val_II))
recall_score(y_val, lr.predict(sc_en_X_val_II))
lr.intercept_
confusion_matrix(y_val, lr.predict(sc_en_X_val_II))
print(classification_report(y_val, lr.predict(sc_en_X_val_II)))

# For binary classification, lr.coef_ returns a 2D array, so we take the first row [0]
coefficients = lr.coef_[0] # <-- Changed from 'model' to 'lr'

# 2. Get the feature names
feature_names = sc_en_X_train_II.columns

# 3. Create a pandas DataFrame
feature_importance_df = pd.DataFrame({
    'Feature': feature_names,
    'Coefficient': coefficients,
    'Absolute_Importance': np.abs(coefficients) # Absolute value for ranking by magnitude
})

# 4. Sort the DataFrame by importance (highest magnitude first)
feature_importance_df = feature_importance_df.sort_values(
    by='Absolute_Importance',
    ascending=False
)

# Display the result

(feature_importance_df)

"""Tenure
- exp(-1.3688) yields approximately 0.254

For every unit increase in tenure (e.g., one month or one year, depending on your scaling), the odds of churning are reduced by roughly 74.6% (\(1-0.254\)), holding all other variables constant.

Complain
- exp(0.5714) yields approximately 1.77

A customer who has complained has 1.77 times (or 77% higher) the odds of churning compared to a customer who has not complained, assuming all other factors are equal.

DecisionTreeClassifier
"""

clf = DecisionTreeClassifier(criterion = 'entropy', class_weight = weights_dict, max_depth = 4, max_features = None,
                             min_samples_split = 25, min_samples_leaf = 15)

# Start with your base dataframe copy
X_train = df_train_II.copy()

# Define the columns you want to remove
columns_to_drop = [

    'Churn'  # Assuming 'Churn' is the name of your target variable column in the full dataframe
]

# Drop the columns to create the final feature set (X_train)
X_train = X_train.drop(columns=columns_to_drop)
# 'errors="ignore"' ensures the code doesn't crash if one of the columns is already missing

# Check the shape of the new X_train to confirm the columns are gone
X_train.shape

# Start with your base dataframe copy
X_val = df_val_II.copy()

# Define the columns you want to remove
columns_to_drop = [

    'Churn'  # Assuming 'Churn' is the name of your target variable column in the full dataframe
]

# Drop the columns to create the final feature set (X_train)
X_val = X_val.drop(columns=columns_to_drop)
# 'errors="ignore"' ensures the code doesn't crash if one of the columns is already missing

# Check the shape of the new X_train to confirm the columns are gone
X_val.shape

# Start with your base dataframe copy
X_test = df_test_II.copy()

# Define the columns you want to remove
columns_to_drop = [

    'Churn'  # Assuming 'Churn' is the name of your target variable column in the full dataframe
]

# Drop the columns to create the final feature set (X_train)
X_test = X_test.drop(columns=columns_to_drop)
# 'errors="ignore"' ensures the code doesn't crash if one of the columns is already missing

# Check the shape of the new X_train to confirm the columns are gone
X_test.shape

if X_val.columns.equals(X_test.columns):
    print("Columns successfully align by name and order. You are ready to model.")
else:
    print("Column names or order still do not match! Review the previous set difference code.")

X_train.shape, y_train.shape
X_val.shape, y_val.shape
X_test.shape, y_test.shape

X_train.columns

X_train.columns

clf.fit(X_train, y_train)

## Training metrics
roc_auc_score(y_train, clf.predict(X_train))
recall_score(y_train, clf.predict(X_train))
confusion_matrix(y_train, clf.predict(X_train))
print(classification_report(y_train, clf.predict(X_train)))

## Validation metrics
roc_auc_score(y_val, clf.predict(X_val))
recall_score(y_val, clf.predict(X_val))
confusion_matrix(y_val, clf.predict(X_val))
print(classification_report(y_val, clf.predict(X_val)))

import pandas as pd

importances = clf.feature_importances_
feature_importance = pd.DataFrame({
    "feature": X_val.columns,
    "importance": importances
}).sort_values(by="importance", ascending=False)

print(feature_importance)

"""Conclusion

- Logist regression and decision tree have been fitted on the imputed data
- The Logistic Regression model does a poor job predicting false positives for the 'Churn' class. The precision score of 43% means that when the model predicts a customer is going to churn, it is only correct 43% of the time.
- The feature importance shows that complain has the highest effect on churn followed by Preffered Payment Mode Cash on Delivery, Prefered Order Cat Mobile. Tenure and City Tier 1 have a negative effect on churn.
- Decicion Tree Classifier has a precision score 46% slightly better than the logistic regression.
- The feature importance of the decision tree classifier shows that Tenure is the most important feauter followed by Order Count, Complain, Order Amount Hike From last Year, Coupon Used, Number of Adress and Day since last order.

### Initial fit non imputed dataframe
"""

df_train_m.shape

# Start with your base dataframe copy
X_train_m = df_train_m.copy()

# Define the columns you want to remove
columns_to_drop = ['Churn']

# Drop the columns to create the final feature set (X_train)
X_train_m = X_train_m.drop(columns=columns_to_drop)
# 'errors="ignore"' ensures the code doesn't crash if one of the columns is already missing

# Check the shape of the new X_train to confirm the columns are gone
X_train_m.shape

# Start with your base dataframe copy
X_val_m = df_val_m.copy()

# Define the columns you want to remove
columns_to_drop = ['Churn']

# Drop the columns to create the final feature set (X_val)
X_val_m = X_val_m.drop(columns=columns_to_drop)
# 'errors="ignore"' ensures the code doesn't crash if one of the columns is already missing

# Check the shape of the new X_val to confirm the columns are gone
X_val_m.shape

# Start with your base dataframe copy
X_test_m = df_test_m.copy()

# Define the columns you want to remove
columns_to_drop = ['Churn']

# Drop the columns to create the final feature set (X_test)
X_test_m = X_test_m.drop(columns=columns_to_drop)
# 'errors="ignore"' ensures the code doesn't crash if one of the columns is already missing

# Check the shape of the new X_test to confirm the columns are gone
X_test_m.shape

if X_val_m.columns.equals(X_test_m.columns):
    print("Columns successfully align by name and order. You are ready to model.")
else:
    print("Column names or order still do not match! Review the previous set difference code.")

X_train_m.shape, y_train_m.shape
X_val_m.shape, y_val_m.shape
X_test_m.shape, y_test_m.shape

clf = DecisionTreeClassifier(criterion = 'entropy', class_weight = weights_dict, max_depth = 4, max_features = None,
                             min_samples_split = 25, min_samples_leaf = 15)

X_train_m.columns

clf.fit(X_train_m, y_train_m)

## Training metrics
roc_auc_score(y_train_m, clf.predict(X_train_m))
recall_score(y_train_m, clf.predict(X_train_m))
confusion_matrix(y_train_m, clf.predict(X_train_m))
print(classification_report(y_train_m, clf.predict(X_train_m)))

## Validation metrics
roc_auc_score(y_val_m, clf.predict(X_val_m))
recall_score(y_val_m, clf.predict(X_val_m))
confusion_matrix(y_val_m, clf.predict(X_val_m))
print(classification_report(y_val_m, clf.predict(X_val_m)))

#  X_train_m is the dataframe you used to TRAIN the model
#  clf is a tree-based model (Random Forest, Decision Tree, etc.)

# 1. Get the feature importances (already a 1D array)
coefficients = clf.feature_importances_

# 2. Get the feature names from your FEATURE dataframe
feature_names = X_train_m.columns # <-- Corrected from y_train_m.columns

# 3. Create a pandas DataFrame
feature_importance_df = pd.DataFrame({
    'Feature': feature_names,
    'Importance': coefficients, # Renamed to 'Importance' for clarity
    # Absolute value isn't needed for feature_importances_ as they are all positive
})

# 4. Sort the DataFrame by importance (highest magnitude first)
feature_importance_df = feature_importance_df.sort_values(
    by='Importance',
    ascending=False
)

# Display the result
(feature_importance_df)

plt.figure(figsize=(10, 6)) # Adjust figure size for better readability

sns.boxplot(
    x='Complain',  # Categories on the x-axis
    y='Tenure', # Numerical values on the y-axis
    hue='Churn',
    data=df_train,
    orient='v'  # 'v' makes the plot vertical, putting categories on the x-axis
)

plt.title('Tenure Distribution by Churn Status and Complaint Status')
plt.xlabel('Churn Status')
plt.ylabel('Tenure (e.g., Months)')
plt.legend(title='Complain Status')
plt.grid(True, axis='y', linestyle='--', alpha=0.6) # Grid on the y-axis for vertical plots
plt.show()

from sklearn import tree

# Assuming 'clf' is your trained DecisionTreeClassifier
# Assuming X_train_m is your feature dataframe used for training

fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(20, 10), dpi=300)

tree.plot_tree(
    clf,
    feature_names=X_train_m.columns, # Your column names
    class_names=['No Churn', 'Churn'], # Your class names
    filled=True,
    rounded=True,
    ax=axes
)

plt.show()
# You can also save it as an image file:
# fig.savefig('my_churn_tree_diagram.png')

"""Conclusion

- Decision tree classifier has been fitted on non imputed and with missing values dataframe.
- Decicion Tree Classifier has a precision score 0.45% which is slightly less than on the imputed dataframe.
- The feature importance of the decision tree classifier shows that Tenure is the most important feauter followed by Complain, NumberOfAddress, CouponUsed, SatisfactionScore, CashbackAmount, MaritalStatus_Single, WarehouseToHome.

### Spot-checking various ML algorithms

__Steps__ :

- Automate data preparation and model run through Pipelines

- Model Zoo : List of all models to compare/spot-check

- Evaluate using k-fold Cross validation framework

__Note__ : Restart the kernel and read the original dataset again followed by train-test split and then come directly to this section of the notebook

#### Automating data preparation and model run through Pipelines

### Model Zoo + k-fold Cross Validation

#### How are models selected ?

 - Why only tree models ? Why not SVM or ANNs?
"""

## Preparing data and a few common model parameters
#X = df_train.drop(columns = ['Exited'], axis = 1)
#y = y_train.ravel()

#weights_dict = {0 : 1.0, 1 : 3.93}
#_, num_samples = np.unique(y_train, return_counts = True)
#weight = (num_samples[0]/num_samples[1]).round(2)
#weight

#cols_to_scale = ['CreditScore', 'Age', 'Balance', 'EstimatedSalary', 'bal_per_product', 'bal_by_est_salary', 'tenure_age_ratio'
#                ,'age_surname_enc']

"""Read more about XGB parameters from here  : https://xgboost.readthedocs.io/en/latest/parameter.html

Tips to tune parameters for LightGBM : https://lightgbm.readthedocs.io/en/latest/Parameters-Tuning.html
"""

## Preparing a list of models to try out in the spot-checking process
def model_zoo(models = dict()):
    # Tree models
    for n_trees in [21, 1001]:
        models['rf_' + str(n_trees)] = RandomForestClassifier(n_estimators = n_trees, n_jobs = -1, criterion = 'entropy'
                                                              , class_weight = weights_dict, max_depth = 6, max_features = 0.6
                                                              , min_samples_split = 30, min_samples_leaf = 20)

        models['lgb_' + str(n_trees)] = LGBMClassifier(boosting_type='dart', num_leaves=31, max_depth= 6, learning_rate=0.1
                                                       , n_estimators=n_trees, class_weight=weights_dict, min_child_samples=20
                                                       , colsample_bytree=0.6, reg_alpha=0.3, reg_lambda=1.0, n_jobs=- 1
                                                       , importance_type = 'gain')

        models['xgb_' + str(n_trees)] = XGBClassifier(objective='binary:logistic', n_estimators = n_trees, max_depth = 6
                                                      , learning_rate = 0.03, n_jobs = -1, colsample_bytree = 0.6
                                                      , reg_alpha = 0.3, reg_lambda = 0.1, scale_pos_weight = imbalance_ratio)

        models['et_' + str(n_trees)] = ExtraTreesClassifier(n_estimators=n_trees, criterion = 'entropy', max_depth = 6
                                                            , max_features = 0.6, n_jobs = -1, class_weight = weights_dict
                                                            , min_samples_split = 30, min_samples_leaf = 20)

      # Naive-Bayes models
    #models['gauss_nb'] = GaussianNB()
    #models['multi_nb'] = MultinomialNB()
    #models['compl_nb'] = ComplementNB()
    #models['bern_nb'] = BernoulliNB()

    return models


def model_zoo_linear(models = dict()):
    # kNN models
    for n in [5,11, 23, 31, 53]:
        models['knn_' + str(n)] = KNeighborsClassifier(n_neighbors=n)

    # RL models
    for n in ["lbfgs", "liblinear", "newton-cg", "newton-cholesky", "sag", "saga"]:
        models['lr_' + str(n)] = LogisticRegression(penalty ='l2', dual=False, tol=0.0001, C=1.0, fit_intercept=True,
                                                    intercept_scaling=1, class_weight=weights_dict, solver=n ,max_iter=100, multi_class='auto',
                                                    verbose=0, warm_start=False, n_jobs=None, l1_ratio=None)

    return models

def evaluate_models(X, y, models, folds = 5, metric = 'recall'):
    results = dict()
    for name, model in models.items():
        # Evaluate model using cross_val_score directly on the model object
        # The 'pipeline' variable is removed as the 'model' itself is used
        scores = cross_val_score(model, X, y, cv = folds, scoring = metric, n_jobs = -1)

        # Store results of the evaluated model
        results[name] = scores
        mu, sigma = np.mean(scores), np.std(scores)
        # Printing individual model results
        print('Model {}: mean = {}, std_dev = {}'.format(name, mu, sigma))

    return results

X_train.sample(5)

X_train.columns

sc_en_X_train_II.columns

## Spot-checking in action
models_ = model_zoo_linear()
print('Recall metric')
results = evaluate_models(sc_en_X_train_II, y_train, models_, metric = 'recall')
print('F1-score metric')
results = evaluate_models(sc_en_X_train_II, y_train , models_, metric = 'f1')

## Spot-checking in action
models = model_zoo()
print('Recall metric')
results = evaluate_models(X_train, y_train , models, metric = 'recall') # USE X_train, y_train (no target inside)
print('F1-score metric')
results = evaluate_models(X_train, y_train , models, metric = 'f1')

sc_en_X_train_II_NoOutliers.columns

df_train_II_NoOutliers.drop('Churn', axis=1, inplace=True)

sc_en_X_train_II_NoOutliers.columns

## Spot-checking in action
models_ = model_zoo_linear()
print('Recall metric')
results = evaluate_models(sc_en_X_train_II_NoOutliers, y_train_NoOutliers, models_, metric = 'recall')
print('F1-score metric')
results = evaluate_models(sc_en_X_train_II_NoOutliers, y_train_NoOutliers , models_, metric = 'f1')

## Spot-checking in action
models = model_zoo()
print('Recall metric')
results = evaluate_models(df_train_II_NoOutliers, y_train_NoOutliers , models, metric = 'recall') # USE X_train, y_train (no target inside)
print('F1-score metric')
results = evaluate_models(df_train_II_NoOutliers, y_train_NoOutliers , models, metric = 'f1')

"""Spot-checking applied to missing values dataset"""

X_train_m.columns

## Spot-checking in action with MISSING VALUES DATA FRAME
models = model_zoo()
print('Recall metric')
results = evaluate_models(X_train_m, y_train_m , models, metric = 'recall')
print('F1-score metric')
results = evaluate_models(X_train_m, y_train_m , models, metric = 'f1')

"""Based on the relevant metric, a suitable model can be chosen for further hyperparameter tuning.

XG boost is chosen for further hyperparameter tuning because it has the best performance on recall metric regarding the dataset with missing values.

### Hyperparameter tuning

RandomSearchCV vs GridSearchCV

- Random Search is more suitable for large datasets, with a large number of parameter settings
- Grid Search results in a more precise hyperparameter tuning, thus resulting in better model performance. Intelligent tuning mechanism can also help reduce the time taken in GridSearch by a large factor

- Will optimize on F1 metric. We could easily reach 75% Recall from the default parameters as seen earlier
"""

# Replace 'path/to/your/filename.csv' with the path you copied
#file_path = '/content/drive/MyDrive/Churn Projectpro/df_train_m.csv'
#file_path = '/content/drive/MyDrive/Churn Projectpro/df_val_m.csv'
#file_path = '/content/drive/MyDrive/Churn Projectpro/df_test_m.csv'


# Training df:
#df_train_m_X = pd.read_csv(file_path)

# Val df:

#df_val_m_X = pd.read_csv(file_path)

# Test df:

#df_test_m_X = pd.read_csv(file_path)

# Drop columns:
#drop_features = ['Unnamed: 0', 'CustomerID', 'Churn',]

#df_train_m_X = df_train_m_X.drop(columns=drop_features)
#df_val_m_X = df_val_m_X.drop(columns=drop_features)
#df_test_m_X = df_test_m_X.drop(columns=drop_features)

X_train_m = df_train_m.copy()
X_val_m = df_val_m.copy()
X_test_m = df_test_m.copy()

# Drop columns:
drop_features = ['Churn', "Unnamed: 0","CustomerID"]

X_train_m = X_train_m.drop(columns=drop_features)
X_val_m = X_val_m.drop(columns=drop_features)
X_test_m = X_test_m.drop(columns=drop_features)

X_train_m.head(10).T

X_train_m.describe(include='all') # Describe all numerical columns

(X_train_m.sample(5))

from sklearn.base import BaseEstimator, TransformerMixin

class CategoricalEncoder(BaseEstimator, TransformerMixin):
    """
    Encodes categorical columns using LabelEncoding and OneHotEncoding.
    LabelEncoding is used for binary categorical columns.
    OneHotEncoding is used for columns with <= 10 distinct values.
    Numeric columns are passed through unchanged.
    """

    def __init__(self, cols=None, lcols=None, ohecols=None, reduce_df=False):
        """
        Parameters
        ----------
        cols : list of str
            Columns to encode. Default is to target/one-hot/label encode all categorical columns in the DataFrame.
        reduce_df : bool
            Whether to use reduced degrees of freedom for encoding
            (that is, add N-1 one-hot columns for a column with N
            categories). Default = False
        """

        if isinstance(cols, str):
            self.cols = [cols]
        else:
            self.cols = cols

        if isinstance(lcols, str):
            self.lcols = [lcols]
        else:
            self.lcols = lcols

        if isinstance(ohecols, str):
            self.ohecols = [ohecols]
        else:
            self.ohecols = ohecols

        self.reduce_df = reduce_df
        self.lmaps = None
        self.ohemaps = None
        self.numeric_cols = None


    def fit(self, X, y=None):
        """Fit label/one-hot encoder to X.
        """

        # Identify numeric columns to keep later and exclude them from self.cols handling
        self.numeric_cols = [c for c in X if str(X[c].dtype) != 'object' and (self.cols is None or c not in self.cols)]

        # Encode all remaining categorical cols by default if self.cols was None
        if self.cols is None:
            self.cols = [c for c in X if str(X[c].dtype)=='object']

        # Check columns are in X
        for col in self.cols:
            if col not in X:
                raise ValueError('Column \''+col+'\' not in X')

        # Separating out lcols and ohecols
        if self.lcols is None:
            self.lcols = [c for c in self.cols if X[c].nunique() <= 2]

        if self.ohecols is None:
            # The remaining columns in self.cols that weren't binary must be OHE
            self.ohecols = [c for c in self.cols if c not in self.lcols]

        ## Create Label Encoding mapping
        self.lmaps = dict()
        for col in self.lcols:
            self.lmaps[col] = dict(zip(X[col].values, X[col].astype('category').cat.codes.values))

        ## Create OneHot Encoding mapping
        self.ohemaps = dict()
        for col in self.ohecols:
            self.ohemaps[col] = []
            uniques = X[col].unique()
            for unique in uniques:
                self.ohemaps[col].append(unique)
            if self.reduce_df:
                del self.ohemaps[col][-1]

        return self


    def transform(self, X, y=None):
        """Perform label/one-hot encoding transformation.
        Returns a DataFrame with numeric columns passed through and
        categorical columns encoded.
        """

        # Initialize Xo with all columns that need to be in the final DF:
        # 1. The numeric columns we saved during fit()
        # 2. The categorical columns which we will transform/replace
        Xo = X[self.numeric_cols + self.cols].copy()

        ## Perform label encoding transformation (modifies in place)
        for col, lmap in self.lmaps.items():
            Xo[col] = Xo[col].map(lmap)
            Xo[col].fillna(-1, inplace=True)

        ## Perform one-hot encoding transformation (adds new columns, deletes old)
        for col, vals in self.ohemaps.items():
            for val in vals:
                new_col = col+'_'+str(val)
                Xo[new_col] = (Xo[col]==val).astype('uint8')
            del Xo[col] # Remove the original categorical column

        return Xo

from sklearn.base import BaseEstimator, TransformerMixin

class CategoricalEncoder(BaseEstimator, TransformerMixin):
    # ... (__init__ remains unchanged) ...
    def __init__(self, cols=None, lcols=None, ohecols=None, reduce_df=False):
        if isinstance(cols, str):
            self.cols = [cols]
        else:
            self.cols = cols

        if isinstance(lcols, str):
            self.lcols = [lcols]
        else:
            self.lcols = lcols

        if isinstance(ohecols, str):
            self.ohecols = [ohecols]
        else:
            self.ohecols = ohecols

        self.reduce_df = reduce_df
        self.lmaps = None
        self.ohemaps = None
        self.numeric_cols = None

    def fit(self, X, y=None):
        """Fit label/one-hot encoder to X."""

        # --- MODIFICATION START (FIT) ---
        # Work on a copy and force specified columns to be 'object' dtype for consistent handling
        X_processed = X.copy()
        if self.cols:
            for col in self.cols:
                if col in X_processed.columns:
                    X_processed[col] = X_processed[col].astype('object')
        # --- MODIFICATION END (FIT) ---

        # Identify numeric columns to keep later and exclude them from self.cols handling
        # This now correctly identifies truly numeric features not in 'self.cols'
        self.numeric_cols = [c for c in X_processed if str(X_processed[c].dtype) != 'object' and (self.cols is None or c not in self.cols)]

        # Encode all remaining categorical cols by default if self.cols was None
        if self.cols is None:
            self.cols = [c for c in X_processed if str(X_processed[c].dtype)=='object']

        # Check columns are in X
        for col in self.cols:
            if col not in X_processed:
                raise ValueError('Column \''+col+'\' not in X')

        # Separating out lcols and ohecols
        if self.lcols is None:
            self.lcols = [c for c in self.cols if X_processed[c].nunique() <= 2]

        if self.ohecols is None:
            # The remaining columns in self.cols that weren't binary must be OHE
            self.ohecols = [c for c in self.cols if c not in self.lcols]

        ## Create Label Encoding mapping
        self.lmaps = dict()
        for col in self.lcols:
            self.lmaps[col] = dict(zip(X_processed[col].values, X_processed[col].astype('category').cat.codes.values))

        ## Create OneHot Encoding mapping
        self.ohemaps = dict()
        for col in self.ohecols:
            self.ohemaps[col] = []
            uniques = X_processed[col].unique()
            for unique in uniques:
                self.ohemaps[col].append(unique)
            if self.reduce_df:
                del self.ohemaps[col][-1]

        return self

    def transform(self, X, y=None):
        """Perform label/one-hot encoding transformation."""

        # --- MODIFICATION START (TRANSFORM) ---
        # Ensure that columns targeted for encoding are cast to object type in the input X
        Xo = X.copy()
        if self.cols:
             for col in self.cols:
                if col in Xo.columns:
                    Xo[col] = Xo[col].astype('object')
        # --- MODIFICATION END (TRANSFORM) ---

        # Initialize Xo with all columns that need to be in the final DF:
        # 1. The numeric columns we saved during fit()
        # 2. The categorical columns which we will transform/replace
        Xo = Xo[self.numeric_cols + self.cols].copy()

        ## Perform label encoding transformation (modifies in place)
        for col, lmap in self.lmaps.items():
            Xo[col] = Xo[col].map(lmap)
            Xo[col].fillna(-1, inplace=True)

        ## Perform one-hot encoding transformation (adds new columns, deletes old)
        for col, vals in self.ohemaps.items():
            for val in vals:
                new_col = col+'_'+str(val)
                Xo[new_col] = (Xo[col]==val).astype('uint8')
            del Xo[col] # Remove the original categorical column

        return Xo

# Import the custom encoder class (assuming you have defined it above)
# from your_module import CategoricalEncoder

# 1. Instantiate the Encoder
# It will automatically find all 'object' dtype columns and apply the default logic
# (Label encoding for binary, OHE for <=10 unique values).
# You can optionally set reduce_df=True to save memory/columns.
encoder = CategoricalEncoder(reduce_df=False)

# 2. Fit and Transform the Training Data
# The .fit_transform() method learns the mappings (e.g., which values map to which OHE columns)
# and applies the transformation immediately.
X_train_m_processed = encoder.fit_transform(X_train_m, y_train_m)

# 3. Transform the Test Data
# The .transform() method uses the mappings already learned from the training data to transform the test set.
# You must NOT call .fit() or .fit_transform() on the test set.
#X_test_processed = encoder.transform(X_test)

# --- Verification Steps (Optional) ---
print("Original X_train_m shape:", X_train_m.shape)
print("Processed X_train_m shape:", X_train_m_processed.shape)
print("\nData Types in Processed Training Set:")
print(X_train_m_processed.info())

parameters = {
    'classifier__n_estimators':[10, 51, 100, 201, 350, 501],
    'classifier__max_depth': [3, 4, 6, 9],
    'classifier__learning_rate': [0.03, 0.05, 0.1, 0.5], # 0.5 and 1.0 are often too high
    'classifier__subsample': [0.7, 0.9, 1.0],
    'classifier__colsample_bytree': [0.3, 0.6, 0.8, 1.0],
    'classifier__min_child_weight': [1, 5, 10], # Added this common parameter
    'classifier__reg_alpha': [0, 0.3, 1, 5],
    'classifier__reg_lambda': [0.1, 1, 5, 10],
    'classifier__gamma': [0, 0.3, 1, 5],
    'classifier__scale_pos_weight': [2 ,3.0, 5.024324324324325, 7.0]}


    # include also min_child_weight, gamma, subsample

# 2. Instantiate the actual XGBoost model object with base parameters
# set use_label_encoder=False and set an objective/eval_metric for best practice.
xgb_base = xgb.XGBClassifier(
    objective='binary:logistic',
    eval_metric='logloss',
    use_label_encoder=False,
    random_state=42 # Good practice for reproducibility
)

# 3. Define the Pipeline with the ENCODER class and the MODEL OBJECT
pipeline = Pipeline(steps = [
    ('categorical_encoding', CategoricalEncoder(reduce_df=False)),
    ('classifier', xgb_base) # <-- Use the actual xgb_base object here


])

#pipeline.set_params(
#    categorical_encoding__cols=['CityTier'], # Specify which columns to process
#    categorical_encoding__ohecols=['CityTier']        # Specifically OHE CityTier
#)

search = RandomizedSearchCV(
    estimator=pipeline,           # Use the 'pipeline' object here
    param_distributions=parameters, # Use the parameters dictionary (which uses 'classifier__' prefixes)
    n_iter=20,                    # Number of iterations
    cv=5,                         # 5-fold cross validation
    scoring='f1',                 # Using F1 score is a good choice for imbalanced classification
    random_state=42,              # For reproducibility
    n_jobs=-1                     # Use all cores
)

X_train_m.columns

# 5. Fit the model and perform the search
print("Starting Randomized Search Fit...")
search.fit(X_train_m, y_train_m.ravel())
print("Search Complete.")

search.best_params_
search.best_score_

search.cv_results_

# Assuming you have X_val_m and y_val_m DataFrames/arrays available

best_pipeline = search.best_estimator_ # This is your finalized model

print("\n--- Evaluation on UNSEEN VALIDATION Data ---")

# Use the best pipeline to generate predictions on the VALIDATION data
# The pipeline automatically ensures the correct transformations are applied
y_val_m_pred = best_pipeline.predict(X_val_m)

# Calculate and print metrics for the validation set
val_roc_auc = roc_auc_score(y_val_m, y_val_m_pred)
val_recall = recall_score(y_val_m, y_val_m_pred)
val_conf_matrix = confusion_matrix(y_val_m, y_val_m_pred)
val_class_report = classification_report(y_val_m, y_val_m_pred)

print(f"ROC AUC Score (Validation): {val_roc_auc:.4f}")
print(f"Recall Score (Validation): {val_recall:.4f}")
print("Confusion Matrix (Validation):")
print(val_conf_matrix)
print("\nClassification Report (Validation):")
print(val_class_report)

"""#### Halving Search"""

parameters = {
    'classifier__n_estimators':[100, 300],
    'classifier__max_depth': [3, 4, 6],
    'classifier__learning_rate': [0.1], # 0.5 and 1.0 are often too high
    'classifier__subsample': [0.7, 0.9],
    'classifier__colsample_bytree': [0.6, 0.8],
    'classifier__min_child_weight': [1, 5, 10], # Added this common parameter
    'classifier__reg_alpha': [0, 1],
    'classifier__reg_lambda': [1, 5, 10],
    'classifier__gamma': [0, 0.3, 1],
    'classifier__scale_pos_weight': [2]}

# Instantiate GridSearchCV correctly
grid_search = GridSearchCV(
    estimator=pipeline,           # Use the 'pipeline' object here
    param_grid=parameters,        # Use 'param_grid' argument for GridSearchCV
    # n_iter and random_state are removed
    cv=5,                         # 5-fold cross validation
    scoring='f1',                 # Using F1 score
    n_jobs=-1                     # Use all cores
)

from sklearn.experimental import enable_halving_search_cv  # Required for experimental feature
from sklearn.model_selection import HalvingGridSearchCV

halving_search = HalvingGridSearchCV(
    estimator=pipeline,
    param_grid=parameters,
    factor=3,           # Number of candidates eliminated in each round
    resource='n_samples',
    scoring='f1',
    cv=5,
    n_jobs=-1,
    verbose=1,
    min_resources=300,
    aggressive_elimination=True,
)

# Fit the model and perform the exhaustive search
print("Starting halving_search Search Fit")
halving_search.fit(X_train_m, y_train_m.ravel())
print("halving_search Search Complete.")

# View the best results
print(f"Best F1 Score (CV average): {halving_search.best_score_}")
print(f"Best Parameters: {halving_search.best_params_}")

halving_search.cv_results_

# Assuming you have X_val_m and y_val_m DataFrames/arrays available

best_pipeline = search.best_estimator_ # This is your finalized model

print("\n--- Evaluation on UNSEEN VALIDATION Data ---")

# Use the best pipeline to generate predictions on the VALIDATION data
# The pipeline automatically ensures the correct transformations are applied
y_val_m_pred = best_pipeline.predict(X_val_m)

# Calculate and print metrics for the validation set
val_roc_auc = roc_auc_score(y_val_m, y_val_m_pred)
val_recall = recall_score(y_val_m, y_val_m_pred)
val_conf_matrix = confusion_matrix(y_val_m, y_val_m_pred)
val_class_report = classification_report(y_val_m, y_val_m_pred)

print(f"ROC AUC Score (Validation): {val_roc_auc:.4f}")
print(f"Recall Score (Validation): {val_recall:.4f}")
print("Confusion Matrix (Validation):")
print(val_conf_matrix)
print("\nClassification Report (Validation):")
print(val_class_report)

# Assuming 'grid_search' is the name of your fitted GridSearchCV object
best_pipeline = halving_search.best_estimator_

# 1. Extract the trained XGBoost classifier object from the pipeline
best_xgb_model = best_pipeline.named_steps['classifier']

# 2. Get the feature importances array
feature_importances = best_xgb_model.feature_importances_

# 3. CRITICAL STEP: Get the column names *after* the encoder runs on the training data.
# We apply the transform method of the ENCODER step on the original X_train data
# to get the exact DataFrame structure the model was trained on.

encoder_step = best_pipeline.named_steps['categorical_encoding']
X_train_transformed_columns = encoder_step.transform(X_train_m).columns

feature_names = X_train_transformed_columns # Use these as your names

# 4. Create the DataFrame using the correct lengths
importance_df = pd.DataFrame({
    'Feature': feature_names,
    'Importance': feature_importances
})

# Sort and display as before
importance_df = importance_df.sort_values(by='Importance', ascending=False)
print("Top 10 Feature Importances (Gain):")
print(importance_df.head(40))

# Plotting the results
plt.figure(figsize=(10, 6))
sns.barplot(x='Importance', y='Feature', data=importance_df.head(50))
plt.title('Top 15 Feature Importances')
plt.tight_layout()
plt.show()

"""### Can we do better? - Ensembles"""

parameters = {
    'classifier__n_estimators':[10, 51, 100, 201, 350, 501],
    'classifier__max_depth': [3, 4, 6, 9],
    'classifier__learning_rate': [0.03, 0.05, 0.1, 0.5], # 0.5 and 1.0 are often too high
    'classifier__subsample': [0.7, 0.9, 1.0],
    'classifier__colsample_bytree': [0.3, 0.6, 0.8, 1.0],
    'classifier__min_child_weight': [1, 5, 10], # Added this common parameter
    'classifier__reg_alpha': [0, 0.3, 1, 5],
    'classifier__reg_lambda': [0.1, 1, 5, 10],
    'classifier__gamma': [0, 0.3, 1, 5],
    'classifier__scale_pos_weight': [2 ,3.0, 5.024324324324325, 7.0]}

#{'classifier__subsample': 0.9,
#'classifier__scale_pos_weight': 2,
#'classifier__reg_lambda': 10,
#'classifier__reg_alpha': 0.3,
#'classifier__n_estimators': 501,
# 'classifier__min_child_weight': 1,
# 'classifier__max_depth': 9,
# 'classifier__learning_rate': 0.5,
# 'classifier__gamma': 0.3,
# 'classifier__colsample_bytree': 0.6}
#np.float64(0.8127831036926777)

## Three versions of the final model with best params for F1-score metric

# Equal weights to both target classes (no class imbalance correction)

xgb1 = xgb.XGBClassifier(
    objective='binary:logistic',
    scale_pos_weight=1,  # Different than grid Search
    colsample_bytree=0.6,
    learning_rate=0.5,
    max_depth=9,
    n_estimators=501,
    reg_alpha=0.3,
    reg_lambda=10,
    subsample=0.9,
    min_child_weight=1,
    gamma=0.3,
    use_label_encoder=False,
    eval_metric='logloss',
    random_state=42
)

# Addressing class imbalance completely by weighting the undersampled class by the class imbalance ratio

xgb2 = xgb.XGBClassifier(
    objective='binary:logistic',
    scale_pos_weight=10,  # Different than grid Search
    colsample_bytree=0.6,
    learning_rate=0.5,
    max_depth=9,
    n_estimators=501,
    reg_alpha=0.3,
    reg_lambda=10,
    subsample=0.9,
    min_child_weight=1,
    gamma=0.3,
    use_label_encoder=False,
    eval_metric='logloss',
    random_state=42
)



# Best class_weight parameter settings (partial class imbalance correction)

xgb3 = xgb.XGBClassifier(
    objective='binary:logistic',
    scale_pos_weight=2,  # Same as best params grid Search
    colsample_bytree=0.6,
    learning_rate=0.5,
    max_depth=9,
    n_estimators=501,
    reg_alpha=0.3,
    reg_lambda=10,
    subsample=0.9,
    min_child_weight=1,
    gamma=0.3,
    use_label_encoder=False,
    eval_metric='logloss',
    random_state=42
)

# Assuming CategoricalEncoder class and xgb1, xgb2, xgb3 objects are defined

# --- Model 1 Definition and Parameter Setting ---
model_1 = Pipeline(steps = [
    ('categorical_encoding', CategoricalEncoder(reduce_df=False)),
    ('classifier', xgb1)
])

# Set parameters for model_1 (must be a separate line after definition)
#model_1.set_params(
#    categorical_encoding__cols=['CityTier'],
#    categorical_encoding__ohecols=['CityTier']
#)


# --- Model 2 Definition and Parameter Setting ---
model_2 = Pipeline(steps = [
    ('categorical_encoding', CategoricalEncoder(reduce_df=False)),
    ('classifier', xgb2)
])

# Set parameters for model_2 (must use the model_2 variable, not model_1)
#model_2.set_params(
#    categorical_encoding__cols=['CityTier'],
#    categorical_encoding__ohecols=['CityTier']
#)


# --- Model 3 Definition and Parameter Setting ---
model_3 = Pipeline(steps = [
    ('categorical_encoding', CategoricalEncoder(reduce_df=False)),
    ('classifier', xgb3)
])

# Set parameters for model_3 (must use the model_3 variable, not model_1)
#model_3.set_params(
#    categorical_encoding__cols=['CityTier'],
#    categorical_encoding__ohecols=['CityTier']
#)

model_1.fit(X_train_m, y_train_m.ravel())
model_2.fit(X_train_m, y_train_m.ravel())
model_3.fit(X_train_m, y_train_m.ravel())

## Getting prediction probabilities from each of these models
m1_pred_probs_trn = model_1.predict_proba(X_train_m)
m2_pred_probs_trn = model_2.predict_proba(X_train_m)
m3_pred_probs_trn = model_3.predict_proba(X_train_m)

## Checking correlations between the predictions of the 3 models
df_t = pd.DataFrame({'m1_pred': m1_pred_probs_trn[:,1], 'm2_pred': m2_pred_probs_trn[:,1], 'm3_pred': m3_pred_probs_trn[:,1]})
df_t.shape
df_t.corr()

"""Although models m1 and m2 are highly correlated (0.9), they are still less closely associated than m2 and m3.
Thus, we'll try to form an ensemble of m1 and m2 (model averaging/stacking) and see if that improves the model accuracy
"""

## Getting prediction probabilities from each of these models
m1_pred_probs_val = model_1.predict_proba(X_val_m)
m2_pred_probs_val = model_2.predict_proba(X_val_m)
m3_pred_probs_val = model_3.predict_proba(X_val_m)

threshold = 0.5

## Best model (Model 3) predictions
m3_preds = np.where(m3_pred_probs_val[:,1] >= threshold, 1, 0)

## Model averaging predictions (Weighted average)
m1_m2_preds = np.where(((0.1*m1_pred_probs_val[:,1]) + (0.9*m2_pred_probs_val[:,1])) >= threshold, 1, 0)

## Model 3 (Best model, tuned by GridSearch) performance on tey_val_midation set
roc_auc_score(y_val_m, m3_preds)
recall_score(y_val_m, m3_preds)
confusion_matrix(y_val_m, m3_preds)
print(classification_report(y_val_m, m3_preds))

## Ensemble model prediction on validation set
roc_auc_score(y_val_m, m1_m2_preds)
recall_score(y_val_m, m1_m2_preds)
confusion_matrix(y_val_m, m1_m2_preds)
print(classification_report(y_val_m, m1_m2_preds))

"""### Hyperparameter tuning dropped features

#### RandomizedSearch
"""

# Dropping Tenure first after 'Complain', NumberOfAddress
#drop_columns = ["Complain", "Tenure", "CityTier" ,"NumberOfAddress", "SatisfactionScore", "CouponUsed",
#                "DaySinceLastOrder", "NumberOfDeviceRegistered", "CashbackAmount", "OrderCount",
#                "WarehouseToHome"]
#X_train_m_drop = X_train_m.drop(columns=drop_columns)
#X_val_m_drop = X_val_m.drop(columns=drop_columns)
#3X_test_m_drop = X_test_m.drop(columns=drop_columns)

# Dropping Tenure first after 'Complain', NumberOfAddress
drop_columns = ["Complain", "Tenure"]
X_train_m_drop = X_train_m.drop(columns=drop_columns)
X_val_m_drop = X_val_m.drop(columns=drop_columns)
X_test_m_drop = X_test_m.drop(columns=drop_columns)

X_train_m_drop.columns

# 2. Instantiate the actual XGBoost model object with base parameters
# set use_label_encoder=False and set an objective/eval_metric for best practice.
xgb_base = xgb.XGBClassifier(
    objective='binary:logistic',
    eval_metric='logloss',
    use_label_encoder=False,
    random_state=42 # Good practice for reproducibility
)

# 3. Define the Pipeline with the ENCODER class and the MODEL OBJECT
pipeline = Pipeline(steps = [
    ('categorical_encoding', CategoricalEncoder(reduce_df=False)),
    ('classifier', xgb_base) # <-- Use the actual xgb_base object here


])

#pipeline.set_params(
#    categorical_encoding__cols=['CityTier'], # Specify which columns to process
#    categorical_encoding__ohecols=['CityTier']        # Specifically OHE CityTier
#)

search = RandomizedSearchCV(
    estimator=pipeline,           # Use the 'pipeline' object here
    param_distributions=parameters, # Use the parameters dictionary (which uses 'classifier__' prefixes)
    n_iter=20,                    # Number of iterations
    cv=5,                         # 5-fold cross validation
    scoring='f1',                 # Using F1 score is a good choice for imbalanced classification
    random_state=42,              # For reproducibility
    n_jobs=-1                     # Use all cores
)

# 5. Fit the model and perform the search
print("Starting Randomized Search Fit...")
search.fit(X_train_m_drop, y_train_m.ravel())
print("Search Complete.")

search.best_params_
search.best_score_

search.cv_results_

# X_val_m and y_val_m DataFrames/arrays available

best_pipeline = search.best_estimator_ # This is your finalized model

print("\n--- Evaluation on UNSEEN VALIDATION Data ---")

# Use the best pipeline to generate predictions on the VALIDATION data
# The pipeline automatically ensures the correct transformations are applied
y_val_m_pred = best_pipeline.predict(X_val_m_drop)

# Calculate and print metrics for the validation set
val_roc_auc = roc_auc_score(y_val_m, y_val_m_pred)
val_recall = recall_score(y_val_m, y_val_m_pred)
val_conf_matrix = confusion_matrix(y_val_m, y_val_m_pred)
val_class_report = classification_report(y_val_m, y_val_m_pred)

print(f"ROC AUC Score (Validation): {val_roc_auc:.4f}")
print(f"Recall Score (Validation): {val_recall:.4f}")
print("Confusion Matrix (Validation):")
print(val_conf_matrix)
print("\nClassification Report (Validation):")
print(val_class_report)

"""#### Halving Search"""

parameters = {
    'classifier__n_estimators':[100, 300],
    'classifier__max_depth': [3, 4, 6],
    'classifier__learning_rate': [0.1], # 0.5 and 1.0 are often too high
    'classifier__subsample': [0.7, 0.9],
    'classifier__colsample_bytree': [0.6, 0.8],
    'classifier__min_child_weight': [1, 5, 10], # Added this common parameter
    'classifier__reg_alpha': [0, 1],
    'classifier__reg_lambda': [1, 5, 10],
    'classifier__gamma': [0, 0.3, 1],
    'classifier__scale_pos_weight': [2]}

from sklearn.experimental import enable_halving_search_cv  # Required for experimental feature
from sklearn.model_selection import HalvingGridSearchCV

halving_search = HalvingGridSearchCV(
    estimator=pipeline,
    param_grid=parameters,
    factor=3,           # Number of candidates eliminated in each round
    resource='n_samples',
    scoring='f1',
    cv=5,
    n_jobs=-1,
    verbose=1,
    min_resources=300,
    aggressive_elimination=True,
)

# Fit the model and perform the exhaustive search
print("Starting halving_search Fit")
halving_search.fit(X_train_m_drop, y_train_m.ravel())
print("halving_search Search Complete.")

# View the best results
print(f"Best F1 Score (CV average): {halving_search.best_score_}")
print(f"Best Parameters: {halving_search.best_params_}")

# Assuming you have X_val_m and y_val_m DataFrames/arrays available

best_pipeline = search.best_estimator_ # This is your finalized model

print("\n--- Evaluation on UNSEEN VALIDATION Data ---")

# Use the best pipeline to generate predictions on the VALIDATION data
# The pipeline automatically ensures the correct transformations are applied
y_val_m_pred = best_pipeline.predict(X_val_m)

# Calculate and print metrics for the validation set
val_roc_auc = roc_auc_score(y_val_m, y_val_m_pred)
val_recall = recall_score(y_val_m, y_val_m_pred)
val_conf_matrix = confusion_matrix(y_val_m, y_val_m_pred)
val_class_report = classification_report(y_val_m, y_val_m_pred)

print(f"ROC AUC Score (Validation): {val_roc_auc:.4f}")
print(f"Recall Score (Validation): {val_recall:.4f}")
print("Confusion Matrix (Validation):")
print(val_conf_matrix)
print("\nClassification Report (Validation):")
print(val_class_report)

# Assuming 'grid_search' is the name of your fitted GridSearchCV object
best_pipeline = halving_search.best_estimator_

# 1. Extract the trained XGBoost classifier object from the pipeline
best_xgb_model = best_pipeline.named_steps['classifier']

# 2. Get the feature importances array
feature_importances = best_xgb_model.feature_importances_

# 3. CRITICAL STEP: Get the column names *after* the encoder runs on the training data.
# We apply the transform method of the ENCODER step on the original X_train data
# to get the exact DataFrame structure the model was trained on.

encoder_step = best_pipeline.named_steps['categorical_encoding']
X_train_transformed_columns = encoder_step.transform(X_train_m).columns

feature_names = X_train_transformed_columns # Use these as your names

# 4. Create the DataFrame using the correct lengths
importance_df = pd.DataFrame({
    'Feature': feature_names,
    'Importance': feature_importances
})

# Sort and display as before
importance_df = importance_df.sort_values(by='Importance', ascending=False)
print("Top 10 Feature Importances (Gain):")
print(importance_df.head(40))

# Plotting the results
plt.figure(figsize=(10, 6))
sns.barplot(x='Importance', y='Feature', data=importance_df.head(50))
plt.title('Top 15 Feature Importances')
plt.tight_layout()
plt.show()

"""### Can we do better 2.0? - Ensembles"""

#Best Parameters: {'classifier__colsample_bytree': 0.8, 'classifier__learning_rate': 0.1, 'classifier__max_depth': 6, 'classifier__n_estimators': 201, 'classifier__reg_alpha': 1, 'classifier__reg_lambda': 0.1, 'classifier__scale_pos_weight': 5.3}

#{'classifier__subsample': 0.9,
# 'classifier__scale_pos_weight': 2,
# 'classifier__reg_lambda': 10,
# 'classifier__reg_alpha': 0.3,
# 'classifier__n_estimators': 501,
# 'classifier__min_child_weight': 1,
# 'classifier__max_depth': 9,
# 'classifier__learning_rate': 0.5,
# 'classifier__gamma': 0.3,
# 'classifier__colsample_bytree': 0.6}
#np.float64(0.714304173713135)

## Three versions of the final model with best params for F1-score metric

# Equal weights to both target classes (no class imbalance correction)

xgb1 = xgb.XGBClassifier(
    objective='binary:logistic',
    scale_pos_weight=1, # Different than Random Search
    colsample_bytree=0.6,
    learning_rate=0.5,
    max_depth=9,
    n_estimators=501,
    reg_alpha=0.3,
    reg_lambda=10,
    subsample = 0.9,
    gamma = 0.3,
    use_label_encoder=False,
    eval_metric='logloss',
    random_state=42
)

# Addressing class imbalance completely by weighting the undersampled class by the class imbalance ratio

xgb2 = xgb.XGBClassifier(
    objective='binary:logistic',
    scale_pos_weight=10, # Different than Random Search
    colsample_bytree=0.6,
    learning_rate=0.5,
    max_depth=9,
    n_estimators=501,
    reg_alpha=0.3,
    reg_lambda=10,
    subsample = 0.9,
    gamma = 0.3,
    use_label_encoder=False,
    eval_metric='logloss',
    random_state=42
)


# Best class_weight parameter settings (partial class imbalance correction)

xgb3 = xgb.XGBClassifier(
    objective='binary:logistic',
    scale_pos_weight=2, # Same as  Random Search best params
    colsample_bytree=0.6,
    learning_rate=0.5,
    max_depth=9,
    n_estimators=501,
    reg_alpha=0.3,
    reg_lambda=10,
    subsample = 0.9,
    gamma = 0.3,
    use_label_encoder=False,
    eval_metric='logloss',
    random_state=42
)

# Assuming CategoricalEncoder class and xgb1, xgb2, xgb3 objects are defined

# --- Model 1 Definition and Parameter Setting ---
model_1 = Pipeline(steps = [
    ('categorical_encoding', CategoricalEncoder(reduce_df=False)),
    ('classifier', xgb1)
])

# Set parameters for model_1 (must be a separate line after definition)
#model_1.set_params(
#   categorical_encoding__cols=['CityTier'],
#   categorical_encoding__ohecols=['CityTier']
#)


# --- Model 2 Definition and Parameter Setting ---
model_2 = Pipeline(steps = [
    ('categorical_encoding', CategoricalEncoder(reduce_df=False)),
    ('classifier', xgb2)
])

# Set parameters for model_2 (must use the model_2 variable, not model_1)
#model_2.set_params(
#   categorical_encoding__cols=['CityTier'],
#   categorical_encoding__ohecols=['CityTier']
#)


# --- Model 3 Definition and Parameter Setting ---
model_3 = Pipeline(steps = [
    ('categorical_encoding', CategoricalEncoder(reduce_df=False)),
    ('classifier', xgb3)# ('classifier', xgb3)
])

# Set parameters for model_3 (must use the model_3 variable, not model_1)
#model_3.set_params(
#   categorical_encoding__cols=['CityTier'],    categorical_encoding__ohecols=['CityTier']
#)

X_train_m_drop.columns

model_1.fit(X_train_m_drop, y_train_m.ravel())
model_2.fit(X_train_m_drop, y_train_m.ravel())
model_3.fit(X_train_m_drop, y_train_m.ravel())

## Getting prediction probabilities from each of these models
m1_pred_probs_trn = model_1.predict_proba(X_train_m_drop)
m2_pred_probs_trn = model_2.predict_proba(X_train_m_drop)
m3_pred_probs_trn = model_3.predict_proba(X_train_m_drop)

## Checking correlations between the predictions of the 3 models
df_t = pd.DataFrame({'m1_pred': m1_pred_probs_trn[:,1], 'm2_pred': m2_pred_probs_trn[:,1], 'm3_pred': m3_pred_probs_trn[:,1]})
df_t.shape
df_t.corr()

## Getting prediction probabilities from each of these models
m1_pred_probs_val = model_1.predict_proba(X_val_m_drop)
m2_pred_probs_val = model_2.predict_proba(X_val_m_drop)
m3_pred_probs_val = model_3.predict_proba(X_val_m_drop)

threshold = 0.5

## Best model (Model 3) predictions
m3_preds = np.where(m3_pred_probs_val[:,1] >= threshold, 1, 0)

## Model averaging predictions (Weighted average)
m1_m2_preds = np.where(((0.1*m1_pred_probs_val[:,1]) + (0.9*m2_pred_probs_val[:,1])) >= threshold, 1, 0)

## Model 3 (Best model, tuned by GridSearch) performance on tey_val_midation set
roc_auc_score(y_val_m, m3_preds)
recall_score(y_val_m, m3_preds)
confusion_matrix(y_val_m, m3_preds)
print(classification_report(y_val_m, m3_preds))

## Ensemble model prediction on validation set
roc_auc_score(y_val_m, m1_m2_preds)
recall_score(y_val_m, m1_m2_preds)
confusion_matrix(y_val_m, m1_m2_preds)
print(classification_report(y_val_m, m1_m2_preds))

# Validation metrics
y_val_pred = model_3.predict(X_val_m_drop)
y_val_prob = model_3.predict_proba(X_val_m_drop)[:, 1]

roc_auc_val = roc_auc_score(y_val_m, y_val_prob)
recall_val = recall_score(y_val_m, y_val_pred)
cm_val = confusion_matrix(y_val_m, y_val_pred)
report_val = classification_report(y_val_m, y_val_pred)

print("ROC AUC (val):", roc_auc_val)
print("Recall (val, class 1):", recall_val)
print("Confusion Matrix (val):\n", cm_val)
print("Classification Report (val):\n", report_val)

"""Conclusion

- Excluding tenure and complain feature foreces the model to utilzes other featufres

- The accuracy metrics drops but still has a good enough recall and f1-score

### Partial Dependence Plots
"""

df_train_m_X.columns

df_train_m.columns

features_to_plot = [
    'PreferedOrderCat', # Categorical feature
    'MaritalStatus', #Categorical feature
    "PreferredPaymentMode", #Categorical feature
    "MaritalStatus_Single", # Numeric features
    "MaritalStatus_Married", # Categorical feature
    "PreferredPaymentMode_E wallet" # Categorical feature
]

features=features_to_plot,

categorical_features_indices_or_names = ['PreferedOrderCat', 'MaritalStatus', 'PreferredPaymentMode',
                                         "PreferredLoginDevice", "Gender"]
categorical_features=categorical_features_indices_or_names,

df_train_m_X.columns

model_3_updated = xgb.XGBClassifier(
    objective='binary:logistic',
    scale_pos_weight=7,
    colsample_bytree=0.8,
    learning_rate=0.1,
    max_depth=6,
    n_estimators=201,
    reg_alpha=0,
    reg_lambda=0.1,
    use_label_encoder=False,
    eval_metric='logloss',
    random_state=42
)

xgb3 = xgb.XGBClassifier(
    objective='binary:logistic',
    scale_pos_weight=2, # Same as  Random Search best params
    colsample_bytree=0.6,
    learning_rate=0.5,
    max_depth=9,
    n_estimators=501,
    reg_alpha=0.3,
    reg_lambda=10,
    subsample = 0.9,
    gamma = 0.3,
    use_label_encoder=False,
    eval_metric='logloss',
    random_state=42
)

# Create a list of the 3 columns you want to remove
cols_to_drop = ['Unnamed: 0', 'CustomerID', 'Churn']

# Returns a new DataFrame (Original stays unchanged)
df_train_m_X = df_train_m_X.drop(columns=cols_to_drop)

df_train_m_X.columns

import numpy as np
import xgboost as xgb
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder
from sklearn.pipeline import Pipeline

gender_mapping = {'Male': 1, 'Female': 0}
df_train_m_X['Gender_cat'] = df_train_m_X['Gender'].map(gender_mapping)

# 1. Define feature lists
# These will keep their NaNs for XGBoost's native handling
numerical_with_nans = [
    "WarehouseToHome", "HourSpendOnApp", "NumberOfDeviceRegistered",
    "NumberOfAddress", "OrderAmountHikeFromlastYear", "CouponUsed",
    "OrderCount", "DaySinceLastOrder", "CashbackAmount",
    "SatisfactionScore", "Gender_cat", "Tenure", "Complain"
]

# Categorical features
categorical_features = [
    "PreferredLoginDevice", "PreferredPaymentMode",
    "PreferedOrderCat", "MaritalStatus", "CityTier"
]



# 2. Drop the original 'Gender' string column so the model doesn't see it
df_train_m_X = df_train_m_X.drop(columns=['Gender'])

# 2. Refined Transformer
preprocessor = ColumnTransformer(
    transformers=[
        # Note: handle_unknown='ignore' ensures that if an unknown value (or NaN if pre-processed)
        # is found, the resulting row is all zeros. XGBoost interprets 'all zeros'
        # for a OHE group as the missing state.
        ('cat', OneHotEncoder(handle_unknown='ignore', sparse_output=False), categorical_features),
        ('num', 'passthrough', numerical_with_nans)
    ]
)

# 3. XGBoost Model (using your 2026 optimized parameters)
xgb3 = xgb.XGBClassifier(
    objective='binary:logistic',
    scale_pos_weight=2,
    colsample_bytree=0.6,
    learning_rate=0.5,
    max_depth=9,
    n_estimators=501,
    reg_alpha=0.3,
    reg_lambda=10,
    subsample=0.9,
    gamma=0.3,
    eval_metric='logloss',
    random_state=42
)

# 4. Final Pipeline
model_3_updated = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('classifier', xgb3)
])

# 5. CRITICAL: Handle SatisfactionScore outlier and Categorical NaNs before fit
# Pipelines don't handle outliers or Categorical NaNs automatically without an imputer
df_train_m_X['SatisfactionScore'] = df_train_m_X['SatisfactionScore'].replace(589314, np.nan)
df_train_m_X[categorical_features] = df_train_m_X[categorical_features].fillna('Missing')

model_3_updated.fit(df_train_m_X, y_train_m.ravel())

# Fit the model: This MUST work now
model_3_updated.fit(df_train_m_X, y_train_m.ravel())

import numpy as np
import xgboost as xgb
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder
from sklearn.pipeline import Pipeline

# --- 1. DATA CLEANING (Only existing columns) ---
# Fix the SatisfactionScore outlier
df_train_m_X['SatisfactionScore'] = df_train_m_X['SatisfactionScore'].replace(589314, np.nan)

# Fill Categorical NaNs so OneHotEncoder doesn't crash
categorical_features = ["PreferredLoginDevice", "PreferredPaymentMode", "PreferedOrderCat", "MaritalStatus", "CityTier"]
df_train_m_X[categorical_features] = df_train_m_X[categorical_features].fillna('Missing')

# --- 2. DEFINE FEATURE LISTS ---
# Using Gender_cat as it is already in your Index
numerical_with_nans = [
    "WarehouseToHome", "HourSpendOnApp", "NumberOfDeviceRegistered",
    "NumberOfAddress", "OrderAmountHikeFromlastYear", "CouponUsed",
    "OrderCount", "DaySinceLastOrder", "CashbackAmount",
    "SatisfactionScore", "Gender_cat", "Tenure", "Complain"
]

# --- 3. REFINED TRANSFORMER ---
preprocessor = ColumnTransformer(
    transformers=[
        ('cat', OneHotEncoder(handle_unknown='ignore', sparse_output=False), categorical_features),
        ('num', 'passthrough', numerical_with_nans)
    ],
    remainder='drop'
)

# --- 4. FINAL PIPELINE ---
model_3_updated = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('classifier', xgb3)
])

# Fit the model
model_3_updated.fit(df_train_m_X, y_train_m.ravel())

df_train_m_X.columns

from sklearn.inspection import PartialDependenceDisplay
import matplotlib.pyplot as plt

# Use the same dataframe you used for model_3_updated.fit()
PartialDependenceDisplay.from_estimator(
    estimator=model_3_updated,
    X=df_train_m_X,  # Use the cleaned dataframe with 'Gender_cat'
    features=[
        "CashbackAmount",
        "OrderCount",
        "DaySinceLastOrder",
        "Tenure"
    ],
    grid_resolution=50
)

plt.show()

n_features = 11

fig, ax = plt.subplots(
    nrows=n_features,
    ncols=1,
    figsize=(5, 3.2 * n_features)  # ⬅️ width, height
)

PartialDependenceDisplay.from_estimator(
    model_3_updated,
    df_train_m_X,
    features=[
        "CashbackAmount",
        "OrderCount",
        "DaySinceLastOrder",
        "CityTier",
        "WarehouseToHome",
        "HourSpendOnApp",
        "NumberOfDeviceRegistered",
        "NumberOfAddress",
        "OrderAmountHikeFromlastYear",
        "CouponUsed",
        "Tenure"
    ],
    grid_resolution=50,
    ax=ax
)

plt.subplots_adjust(hspace=1)
plt.show()

disp = PartialDependenceDisplay.from_estimator(
    model_3_updated,
    df_train_m_X,
    features=[
        "CashbackAmount",
        "OrderCount",
        "DaySinceLastOrder",
        "CityTier",
        "WarehouseToHome",
        "HourSpendOnApp",
        "NumberOfDeviceRegistered",
        "NumberOfAddress",
        "OrderAmountHikeFromlastYear",
        "CouponUsed",
        "Tenure"
    ],
    grid_resolution=50
)

disp.figure_.set_size_inches(10, 32)
plt.subplots_adjust(hspace=1)
plt.show()

# The same code without tenure
disp = PartialDependenceDisplay.from_estimator(
    model_3_updated,
    df_train_m_X,
    features=[
        "CashbackAmount",
        "OrderCount",
        "DaySinceLastOrder",
        "CityTier",
        "WarehouseToHome",
        "HourSpendOnApp",
        "NumberOfDeviceRegistered",
        "NumberOfAddress",
        "OrderAmountHikeFromlastYear",
        "CouponUsed"
    ],
    grid_resolution=50
)

disp.figure_.set_size_inches(10, 32)
plt.subplots_adjust(hspace=1)
plt.show()

from sklearn.impute import SimpleImputer
from sklearn.inspection import PartialDependenceDisplay
import matplotlib.pyplot as plt
from itertools import combinations
import numpy as np
from mpl_toolkits.axes_grid1 import make_axes_locatable


# 1. Extract pieces from your existing pipeline
preprocessor = model_3_updated.named_steps['preprocessor']
xgb_model = model_3_updated.named_steps['classifier']

# 2. Define the features
numeric_features = [
    "CashbackAmount", "OrderCount", "DaySinceLastOrder", "WarehouseToHome",
    "HourSpendOnApp", "NumberOfDeviceRegistered", "NumberOfAddress",
    "OrderAmountHikeFromlastYear", "CouponUsed", "Tenure"
]

# 3. Transform and Impute (for plotting stability)
X_transformed = preprocessor.transform(df_train_m_X)
all_feature_names = preprocessor.get_feature_names_out()

imputer = SimpleImputer(strategy='median')
X_imputed = imputer.fit_transform(X_transformed)

# 4. Create the map
feature_map = {name: f"num__{name}" for name in numeric_features}

# 5. Generate combinations using the correct list



from mpl_toolkits.axes_grid1 import make_axes_locatable

for feat1, feat2 in combinations(numeric_features, 2):
    display_features = [(feature_map[feat1], feature_map[feat2])]

    # Increase figsize width even more to 12
    fig, ax = plt.subplots(figsize=(12, 6))

    display = PartialDependenceDisplay.from_estimator(
        xgb_model,
        X_imputed,
        features=display_features,
        feature_names=all_feature_names,
        grid_resolution=30,
        ax=ax
    )

    # 1. Push Legend further right (increase 1.2 to 1.4 or 1.5)
    legend = ax.get_legend()
    if legend:
        # (x, y) coordinates: 1.4 is far to the right of the axis
        legend.set_bbox_to_anchor((1.45, 1.0))

    # 2. Push Colorbar further right (increase pad from 0.5 to 0.8)
    if display.contours_ is not None:
        divider = make_axes_locatable(ax)
        # Size is the thickness of the bar, pad is the distance from plot
        cax = divider.append_axes("right", size="5%", pad=0.8)
        plt.colorbar(display.contours_[0][0], cax=cax)

    ax.set_title(f"2D PDP: {feat1} vs {feat2}", fontsize=12, pad=25)

    # rect=[left, bottom, right, top] - decrease the 'right' boundary to 0.8
    # This squeezes the plot to the left, giving the legend/colorbar more room
    plt.tight_layout(rect=[0, 0, 0.8, 1])
    plt.show()

# In order to visualize Gencer and Complain install sickit-learn 1.5.2

#pip install scikit-learn==1.5.2

#import sklearn
#print(f"Active version: {sklearn.__version__}")

X_pdp = df_train_m_X.copy()


categorical_pdp_features = [
    "PreferredLoginDevice",
    "PreferredPaymentMode",
    "Gender_cat",
    "MaritalStatus",
    "PreferedOrderCat",
    "Complain"
]

for col in categorical_pdp_features:
    X_pdp[col] = X_pdp[col].astype(object).fillna("Missing")

PartialDependenceDisplay.from_estimator(
    model_3_updated,
    X_pdp,
    features=categorical_pdp_features,
    categorical_features=categorical_pdp_features,


)

plt.show()

import matplotlib.pyplot as plt
from sklearn.inspection import PartialDependenceDisplay

n_features = len(categorical_pdp_features)

fig, ax = plt.subplots(
    nrows=n_features,
    ncols=1,
    figsize=(5, 4 * n_features)  #  increase height
)

PartialDependenceDisplay.from_estimator(
    model_3_updated,
    X_pdp,
    features=categorical_pdp_features,
    categorical_features=categorical_pdp_features,
    ax=ax
)

fig.subplots_adjust(hspace=1)  #  vertical spacing
plt.show()

"""Conclusion

- As the partial dependence plots show, Complain and Tenure have huge impacts on churn.

Interpretaions:
- On average, a customer who change worm not complaining to complain the likelyhood of churn increases by about 20 percentage points, holding all other values constant. Given the low correlation (0.20) between this and other features, this represents a reliable isolated effect of customer complaints on churn risk.
- Based on the partial dependence analysis, a customer with 20 months of tenure has a predicted probability of churn that is 47 percentage points lower than a customer with less than one month of tenure, assuming all other customer characteristics remain equal.
- While tenure and complaint are powerful predictors, the partial dependence plots for the other numeric features show non-linear relationships with churn.This non-linearity warrants a deeper, dedicated analysis to understand their specific contributions to churn risk.

### Ablation Study
"""

# --- Model 3 Definition and Parameter Setting ---
model_3 = Pipeline(steps = [
    ('categorical_encoding', CategoricalEncoder(reduce_df=False)),
    ('classifier', xgb3)# ('classifier', xgb3)
])

from sklearn.metrics import f1_score
import copy


def run_ablation(pipeline, X_train_m_drop, y_train, X_val_m_drop, y_val, features_to_drop):
    X_train_ab = X_train_m_drop.drop(columns=features_to_drop)
    X_val_ab = X_val_m_drop.drop(columns=features_to_drop)

    pipe = copy.deepcopy(pipeline)
    pipe.fit(X_train_ab, y_train)

    preds = pipe.predict(X_val_ab)
    return f1_score(y_val, preds, average="binary", pos_label=1)

X_train_m_drop['CityTier'] = X_train_m_drop['CityTier'].map({1: 'Tier1', 2: 'Tier2', 3: 'Tier3'})
X_val_m_drop['CityTier'] = X_val_m_drop['CityTier'].map({1: 'Tier1', 2: 'Tier2', 3: 'Tier3'})

print(X_val_m_drop['CityTier'].unique())

X_train_m_drop.columns

#X_train_m_drop = X_train_m_drop.drop(columns=['Unnamed: 0', 'CustomerID'])

print(X_train_m_drop.columns.tolist())

baseline_f1 = run_ablation(
    model_3,
    X_train_m_drop,
    y_train_m,
    X_val_m_drop,
    y_val_m,
    features_to_drop=[]
)

print("Baseline F1:", baseline_f1)

from sklearn.metrics import f1_score
import copy
import pandas as pd

#With deepcopy: Ensure the model is "born again" as a completely untrained, blank slate for every
#feature you test. This ensures that the drop in F1-score is caused only by the missing feature,
#not by leftover patterns from a previous run, basically a Fresh Start for every turn.

def run_ablation(pipeline, X_train, y_train, X_val, y_val, features_to_drop):
    X_train_ab = X_train.drop(columns=features_to_drop)
    X_val_ab = X_val.drop(columns=features_to_drop)

    pipe = copy.deepcopy(pipeline)
    pipe.fit(X_train_ab, y_train)

    preds = pipe.predict(X_val_ab)
    # F1 score for class 1
    return f1_score(y_val, preds, average="binary", pos_label=1)

# ---- Baseline
baseline_f1 = run_ablation(
    model_3,
    X_train_m_drop,
    y_train_m,
    X_val_m_drop,
    y_val_m,
    features_to_drop=[]
)

print("Baseline F1 (class 1):", baseline_f1)

results = []

for feature in X_train_m_drop.columns:


    f1_class1 = run_ablation(
        model_3,
        X_train_m_drop,
        y_train_m,
        X_val_m_drop,
        y_val_m,
        features_to_drop=[feature]
    )

    results.append({
        "dropped_features": feature,
        "f1_class1": f1_class1,
        "delta_from_baseline": f1_class1 - baseline_f1
    })

ablation_df = pd.DataFrame(results).sort_values("f1_class1")
print(ablation_df)

# CashbackAmount, the F1-score drops by 2.5%
# Gender, your F1-score increases by nearly 5%
# recursive feature elimination could be applied

# backward Elimination or recursive elimination (RFE see sickit learn not applied)

# Dropping Features that decrease F1 score of class 1

# Fixed: Added the missing comma after "PreferredPaymentMode"
cols_to_drop = [
    'DaySinceLastOrder', "NumberOfAddress", "CityTier",
    "PreferedOrderCat", "MaritalStatus", "HourSpendOnApp",
    "CouponUsed", "PreferredLoginDevice", "PreferredPaymentMode", # Comma added here
    "NumberOfDeviceRegistered", "Gender"
]

X_train_m_drop_2 = X_train_m_drop.drop(columns=cols_to_drop).copy()
X_val_m_drop_2 = X_val_m_drop.drop(columns=cols_to_drop).copy()
X_test_m_drop_2 = X_test_m_drop.drop(columns=cols_to_drop).copy()

X_train_m_drop_2.columns

model_3.fit(X_train_m_drop_2, y_train_m.ravel())

# Validation metrics
y_val_pred = model_3.predict(X_val_m_drop_2)
y_val_prob = model_3.predict_proba(X_val_m_drop_2)[:, 1]

roc_auc_val = roc_auc_score(y_val_m, y_val_prob)
recall_val = recall_score(y_val_m, y_val_pred)
cm_val = confusion_matrix(y_val_m, y_val_pred)
report_val = classification_report(y_val_m, y_val_pred)

print("ROC AUC (val):", roc_auc_val)
print("Recall (val, class 1):", recall_val)
print("Confusion Matrix (val):\n", cm_val)
print("Classification Report (val):\n", report_val)

#np.float64(0.8442941764202113)
#0.7168141592920354
#array([[482,  14],
#       [ 32,  81]])
#              precision    recall  f1-score   support

#           0       0.94      0.97      0.95       496
#           1       0.85      0.72      0.78       113

#    accuracy                           0.92       609
#   macro avg       0.90      0.84      0.87       609
#weighted avg       0.92      0.92      0.92       609

import matplotlib.pyplot as plt
from sklearn.model_selection import LearningCurveDisplay

# Use from_estimator to automatically handle cross-validation
LearningCurveDisplay.from_estimator(
    model_3,
    X_train_m_drop_2,
    y_train_m.ravel(),
    cv=5,
    score_name="Accuracy"
)

plt.title("Learning Curve (Check for Overfitting)")
plt.show()

import pandas as pd
import numpy as np
from sklearn.model_selection import PredefinedSplit, LearningCurveDisplay

# 1. Combine sets using pd.concat instead of np.vstack to keep column names
X_combined = pd.concat([X_train_m_drop_2, X_val_m_drop_2])
y_combined = pd.concat([pd.Series(y_train_m.ravel()), pd.Series(y_val_m.ravel())])

# 2. Create the mask (exactly as before)
split_indices = np.concatenate([
    -1 * np.ones(len(X_train_m_drop_2)),
     0 * np.ones(len(X_val_m_drop_2))
])
pds = PredefinedSplit(test_fold=split_indices)

# 3. Plotting should now work
LearningCurveDisplay.from_estimator(
    model_3,
    X_combined,
    y_combined,
    cv=pds,
    score_name="F1"
)

parameters = {
    'classifier__n_estimators':[10, 51, 100, 201, 350, 501],
    'classifier__max_depth': [3, 4, 6, 9],
    'classifier__learning_rate': [0.03, 0.05, 0.1, 0.5], # 0.5 and 1.0 are often too high
    'classifier__subsample': [0.7, 0.9, 1.0],
    'classifier__colsample_bytree': [0.3, 0.6, 0.8, 1.0],
    'classifier__min_child_weight': [1, 5, 10], # Added this common parameter
    'classifier__reg_alpha': [0, 0.3, 1, 5],
    'classifier__reg_lambda': [0.1, 1, 5, 10],
    'classifier__gamma': [0, 0.3, 1, 5],
    'classifier__scale_pos_weight': [2 ,3.0, 5.024324324324325, 7.0]}


    # include also min_child_weight, gamma, subsample

# 2. Instantiate the actual XGBoost model object with base parameters
# set use_label_encoder=False and set an objective/eval_metric for best practice.
xgb_base = xgb.XGBClassifier(
    objective='binary:logistic',
    eval_metric='logloss',
    use_label_encoder=False,
    random_state=42 # Good practice for reproducibility
)

# 3. Define the Pipeline with the ENCODER class and the MODEL OBJECT
pipeline = Pipeline(steps = [
    ('categorical_encoding', CategoricalEncoder(reduce_df=False)),
    ('classifier', xgb_base) # <-- Use the actual xgb_base object here


])

#pipeline.set_params(
#    categorical_encoding__cols=['CityTier'], # Specify which columns to process
#    categorical_encoding__ohecols=['CityTier']        # Specifically OHE CityTier
#)

search = RandomizedSearchCV(
    estimator=pipeline,           # Use the 'pipeline' object here
    param_distributions=parameters, # Use the parameters dictionary (which uses 'classifier__' prefixes)
    n_iter=20,                    # Number of iterations
    cv=5,                         # 5-fold cross validation
    scoring='f1',                 # Using F1 score is a good choice for imbalanced classification
    random_state=42,              # For reproducibility
    n_jobs=-1                     # Use all cores
)

# 5. Fit the model and perform the search
print("Starting Randomized Search Fit...")
search.fit(X_train_m_drop_2, y_train_m.ravel())
print("Search Complete.")

search.best_params_
search.best_score_

{'classifier__subsample': 0.9,
 'classifier__scale_pos_weight': 2,
 'classifier__reg_lambda': 10,
 'classifier__reg_alpha': 0.3,
 'classifier__n_estimators': 501,
 'classifier__min_child_weight': 1,
 'classifier__max_depth': 9,
 'classifier__learning_rate': 0.5,
 'classifier__gamma': 0.3,
 'classifier__colsample_bytree': 0.6}
np.float64(0.714304173713135)

# Assuming you have X_val_m and y_val_m DataFrames/arrays available

best_pipeline = search.best_estimator_ # This is your finalized model

print("\n--- Evaluation on UNSEEN VALIDATION Data ---")

# Use the best pipeline to generate predictions on the VALIDATION data
# The pipeline automatically ensures the correct transformations are applied
y_val_m_pred = best_pipeline.predict(X_val_m)

# Calculate and print metrics for the validation set
val_roc_auc = roc_auc_score(y_val_m, y_val_m_pred)
val_recall = recall_score(y_val_m, y_val_m_pred)
val_conf_matrix = confusion_matrix(y_val_m, y_val_m_pred)
val_class_report = classification_report(y_val_m, y_val_m_pred)

print(f"ROC AUC Score (Validation): {val_roc_auc:.4f}")
print(f"Recall Score (Validation): {val_recall:.4f}")
print("Confusion Matrix (Validation):")
print(val_conf_matrix)
print("\nClassification Report (Validation):")
print(val_class_report)

#ROC AUC (val): 0.926687838995147
#Recall (val, class 1): 0.8141592920353983
#Confusion Matrix (val):
# [[463  33]
# [ 21  92]]
#Classification Report (val):
#               precision    recall  f1-score   support

#           0       0.96      0.93      0.94       496
#           1       0.74      0.81      0.77       113

#    accuracy                           0.91       609
#   macro avg       0.85      0.87      0.86       609
#weighted avg       0.92      0.91      0.91       609

"""The Leave-One-Feature-Out (LOFO) Ablation Study suggests that with only five features, we can achieve superior recall (increasing from 0.72 to 0.81). Beyond the performance gains, this streamlined model is significantly more cost-effective. It reduces the technical overhead of maintaining redundant data streams, lowers computational costs for real-time inference, and focuses our retention strategy on the most impactful behavioral triggers.
To sum up the two features Complain and Tenure have been excluded as well as the further features leaving the model with only five features.

### Error analysis
"""

X_train_m_drop_2.shape
X_val_m_drop_2.shape

## Final model with best params from ablation study above


xgb3 = xgb.XGBClassifier(
    objective='binary:logistic',
    scale_pos_weight=2,
    colsample_bytree=0.6,
    learning_rate=0.5,
    max_depth=9,
    n_estimators=501,
    reg_alpha=0.3,
    reg_lambda=10,
    use_label_encoder=False,
    eval_metric='logloss',
    min_child_weight = 1,
    gamma = 0.3,
    random_state=42
)

print(model_3.get_params())

# --- Model 3 Definition and Parameter Setting ---
model_3 = Pipeline(steps = [
    ('categorical_encoding', CategoricalEncoder(reduce_df=False)),
    ('classifier', xgb3)# ('classifier', xgb3)
])

model_3.fit(X_train_m_drop_2, y_train_m.ravel())

## Making predictions on a copy of validation set
df_ea = X_val_m_drop_2.copy()
df_ea['y_pred'] = model_3.predict(X_val_m_drop_2)
df_ea['y_pred_prob'] = model_3.predict_proba(X_val_m_drop_2)[:,1]

X_val_m_drop_2.columns

df_ea.shape
df_ea.sample(20)

sns.violinplot(
    x=y_val_m.ravel(),
    y=df_ea['y_pred_prob']
)

df_plot = df_ea.copy()
df_plot['y_true'] = y_val_m.ravel()

sns.histplot(
    data=df_plot,
    x='y_pred_prob',
    hue='y_true',
    bins=55,
    kde=False,
    stat='density',
    common_norm=False,
    element='bars'
)

plt.xlabel("Predicted Probability")
plt.title("Histogram of Predicted Probabilities")
plt.show()

import matplotlib.pyplot as plt
import seaborn as sns

df_plot = df_ea.copy()
df_plot['y_true'] = y_val_m.ravel()

# 1. Use high-visibility hex codes
# Neon Green for Class 0 (Negative), Neon Pink/Magenta for Class 1 (Positive)
bright_colors = ['#00FF00', '#FF00FF']

ax = sns.histplot(
    data=df_plot,
    x='y_pred_prob',
    hue='y_true',
    bins=55,
    kde=False,
    stat='density',
    common_norm=False,
    element='step',      # 'step' or 'poly' often makes overlapping areas clearer than bars
    fill=True,           # Ensure the area is filled
    alpha=0.4,           # Lower alpha makes the overlap (the "errors") look like a third color
    palette=bright_colors,
    linewidth=2          # Thicker lines make the distribution boundaries pop
)

leg = ax.get_legend()
# Note: If your legend labels are '0' and '1', these represent 'Actual Negative' and 'Actual Positive'
new_labels = ['Actual Negative (Check for FP)', 'Actual Positive (Check for FN)']

for t, l in zip(leg.texts, new_labels):
    t.set_text(l)

plt.xlabel("Predicted Probability")
plt.title("Distribution of Probabilities (Brightened for Error Analysis)")
plt.grid(axis='y', alpha=0.3) # Adding a subtle grid helps see density differences
plt.show()

"""The true negative side, predicted probilities less than 0.5 (light green area), has log  normal distribution, is very confident but the density decreases and has several data points that are more than 0.5 leading to false negative predictions. On the other hand, true postivie side (light pruple area), is less confident and has a shape that is flatt and more scattered leadeing to predicted probilities that are false positive.

#### Revisiting bivariate plots of important features

The difference in distribution of these features across the two classes help us to test a few hypotheses
"""

# 1. Create the analysis copy
df_ea = X_val_m_drop_2.copy()

# 2. Add the predictions and labels to THIS copy (df_ea)
df_ea['y_pred'] = model_3.predict(X_val_m_drop_2)
df_ea['y_pred_prob'] = model_3.predict_proba(X_val_m_drop_2)[:, 1]
df_ea['y_true'] = y_val_m.ravel()

df_ea.sample(10)

# Individuas predicted to churn but did NOT churn

false_positives = df_ea[
    (df_plot['y_pred'] == 1) &
    (df_plot['y_true'] == 0)
].copy()   #

false_positives['Churn'] = 0

# Individuas predicted not to churn but DID churn
false_negatives = df_ea[
    (df_plot['y_pred'] == 0) &
    (df_plot['y_true'] == 1)
].copy()
false_negatives['Churn'] = 1

false_positives.sample(10)

df_errors = pd.concat(
    [false_positives, false_negatives],
    axis=0,
    ignore_index=True
)

df_errors.sample(20)

"""""
# Define numeric features
num_list = [
    'WarehouseToHome', "NumberOfAddress", "OrderAmountHikeFromlastYear", "OrderCount", "CashbackAmount"
]

# Add error type column
df_errors['error_type'] = np.select(
    [
        (df_errors['y_true'] == 0) & (df_errors['y_pred'] == 0),
        (df_errors['y_true'] == 1) & (df_errors['y_pred'] == 1),
        (df_errors['y_true'] == 0) & (df_errors['y_pred'] == 1),
        (df_errors['y_true'] == 1) & (df_errors['y_pred'] == 0),
    ],
    [
        'True Negative',
        'True Positive',
        'False Positive',
        'False Negative'
    ],
    default='Unknown'  # <- explicitly make it a string
)
"""""


def numeric_features_visuals_by_error_type(df, features):
    sns.set(style="whitegrid")

    for feature in features:
        plt.figure(figsize=(7, 4))
        sns.boxplot(
            data=df,
            x='Churn', # Matches the column created earlier
            y=feature,
            notch=True
        )
        plt.title(f"{feature} by Prediction Error Type")
        plt.xlabel("Prediction Outcome")
        plt.tight_layout()
        plt.show()

num_list = [
    "WarehouseToHome",	"OrderAmountHikeFromlastYear",	"OrderCount",	"CashbackAmount"
]
numeric_features_visuals_by_error_type(df_errors, num_list)

df_errors.columns

cat_feats = [
    'SatisfactionScore'
]

def categorical_features_visuals_by_error_type(df):
    sns.set(style="whitegrid")

    for feature in cat_feats:
        plt.figure(figsize=(8, 4))
        sns.countplot(
            data=df,
            x=feature,
            hue='Churn'
        )
        plt.title(f"{feature} by Prediction Error Type")
        plt.xlabel(feature)
        plt.ylabel("Count")
        plt.xticks(rotation=45)
        plt.legend(title="Prediction Outcome")
        plt.tight_layout()
        plt.show()
categorical_features_visuals_by_error_type(df_errors)

plt.figure(figsize=(7, 4))
sns.boxplot(
    data=df_errors,
    x='SatisfactionScore',
    y='CashbackAmount'
)
plt.title("Cashback Amount by Satisfaction Score")
plt.xlabel("Satisfaction Score")
plt.ylabel("Cashback Amount")
plt.tight_layout()
plt.show()

df_errors.sample(5)

error_subset = df_errors[
    df_errors['Churn'].isin([0, 1])
]

plt.figure(figsize=(9, 4))
sns.boxplot(
    data=error_subset,
    x='SatisfactionScore',
    y='CashbackAmount',
    hue='Churn'
)
plt.title("Cashback Amount by Satisfaction Score (FP vs FN)")
plt.xlabel("Satisfaction Score")
plt.ylabel("Cashback Amount")
plt.legend(title="Error Type")
plt.tight_layout()
plt.show()

# The correct way to filter the DataFrame
SatisfactionScore_3 = df_errors[df_errors["SatisfactionScore"] == 3].copy()

import matplotlib.pyplot as plt
import seaborn as sns

num_list = [
    "WarehouseToHome", "OrderAmountHikeFromlastYear", "OrderCount", "CashbackAmount"
]

# 1. Added () and :
def plot_satisfaction_errors():
    for i in num_list:
        plt.figure(figsize=(7, 4))
        sns.boxplot(
            data=SatisfactionScore_3,
            x='Churn', # This shows actual churners vs non-churners in the '3' segment
            y=i
        )
        plt.title(f"Distribution of {i} for Satisfaction Score 3")
        plt.show() # 2. Added plt.show() to render each plot separately

# 3. Call the function to actually generate the plots
plot_satisfaction_errors()

SatisfactionScore_3

df_errors.sample(10)

sns.boxplot(x = 'y_true', y = 'WarehouseToHome', data = df_errors)

# 1. Create the analysis copy
df_ea = X_val_m_drop_2.copy()

# 2. Add the predictions and labels to THIS copy (df_ea)
df_ea['y_pred'] = model_3.predict(X_val_m_drop_2)
df_ea['y_pred_prob'] = model_3.predict_proba(X_val_m_drop_2)[:, 1]
df_ea['y_true'] = y_val_m.ravel()

## Are we able to correctly identify pockets of high-churn customer regions in feature space?
df_errors.y_true.value_counts(normalize=True).sort_index()
df_errors[(df_errors.WarehouseToHome > 9) & (df_errors.WarehouseToHome < 24)].y_true.value_counts(normalize=True).sort_index()
df_errors[(df_errors.WarehouseToHome > 9) & (df_errors.WarehouseToHome < 24)].y_pred.value_counts(normalize=True).sort_index()

num_feats = [ 'WarehouseToHome', "OrderAmountHikeFromlastYear", "OrderCount", "CashbackAmount"]

# 1. Create the analysis copy
df_ea = X_val_m_drop_2.copy()

# 2. Add the predictions and labels to THIS copy (df_ea)
df_ea['y_pred'] = model_3.predict(X_val_m_drop_2)
df_ea['y_pred_prob'] = model_3.predict_proba(X_val_m_drop_2)[:, 1]
df_ea['y_true'] = y_val_m.ravel()

# 3. Run the correlation on df_ea, which now contains all columns
correlation_matrix = df_ea[num_list + ['y_pred', 'y_true']].corr()

# 4. Display the results
print(correlation_matrix[['y_pred', 'y_true']])

"""#### Extracting the subset of incorrect predictions

All incorrect predictions are extracted and categorized into false positives (low precision) and false negatives (low recall)
"""

df_ea.columns

low_recall = df_ea[(df_ea.y_true == 1) & (df_ea.y_pred == 0)]
low_prec = df_ea[(df_ea.y_true == 0) & (df_ea.y_pred == 1)]
low_recall.shape
low_prec.shape
low_recall.head()
low_prec.head()

## Prediction probabilty distribution of errors causing low recall
sns.distplot(low_recall.y_pred_prob, hist=False)

## Prediction probabilty distribution of errors causing low precision
sns.distplot(low_prec.y_pred_prob, hist=False)

"""#### Tweaking the threshold of classifier"""

threshold = 0.50

X_val_m_drop_2.columns

#X_val_m_drop_2 = X_val_m_drop_2.drop(columns=['Churn', "y_true"]).copy()

## Predict on validation set with adjustable decision threshold
probs = model_3.predict_proba(X_val_m_drop_2)[:,1]
val_preds = np.where(probs > threshold, 1, 0)

## Default params : 0.5 threshold
confusion_matrix(y_val_m, val_preds)
print(classification_report(y_val_m, val_preds))

## Tweaking threshold between 0.4 and 0.6
confusion_matrix(y_val_m, val_preds)
print(classification_report(y_val_m, val_preds))

gap_samples = probs[(probs > 0.50) & (probs <= 0.55)]
print(f"Samples between 0.50 and 0.55: {len(gap_samples)}")

"""#### Checking whether there's too much dependence on certain features

We'll compare a few important features : NumOfProducts, IsActiveMember, Age, Balance
"""

df_ea.columns

df_ea.SatisfactionScore.value_counts(normalize=True).sort_index()
low_recall.SatisfactionScore.value_counts(normalize=True).sort_index()
low_prec.SatisfactionScore.value_counts(normalize=True).sort_index()

"""### Train final, best model ; Save model and its parameters"""

# Define the constraints using the exact column names after One-Hot Encoding

xgb3 = xgb.XGBClassifier(
    objective='binary:logistic',
    scale_pos_weight=2.5, # Same as  Random Search best params best  weight so far 2.5
    colsample_bytree=0.6,
    learning_rate=0.5,
    max_depth=9,
    n_estimators=501,
    reg_alpha=0.3,
    reg_lambda=10,
    subsample = 0.9,
    gamma = 0.3,
    use_label_encoder=False,
    eval_metric='logloss',
    random_state=42
)

# Your pipeline
model_3 = Pipeline(steps = [
    ('categorical_encoding', CategoricalEncoder(reduce_df=False)),
    ('classifier', xgb3)
])

# To train with your custom weights for CityTier and PaymentMode:
# model_3.fit(X_train_m_drop_2, y_train_m.ravel(), classifier__sample_weight=weights_array)

model_3.fit(X_train_m_drop_2, y_train_m.ravel())

# Predict target probabilities
val_probs = model_3.predict_proba(X_val_m_drop_2)[:,1]

# Predict target values on val data
val_preds = np.where(val_probs > 0.65, 1, 0) # The probability threshold can be tweaked

# Explicitly name the x and y arguments
sns.boxplot(x=y_val_m.ravel(), y=val_probs)

# Adding labels for clarity
plt.xlabel("True Label (0 or 1)")
plt.ylabel("Predicted Probability")
plt.title("Distribution of Probabilities by Actual Class")
plt.axhline(0.65, color='red', linestyle='--') # Optional: visually show your threshold
plt.show()

# Validation metrics
y_val_pred = model_3.predict(X_val_m_drop_2)
y_val_prob = model_3.predict_proba(X_val_m_drop_2)[:, 1]

roc_auc_val = roc_auc_score(y_val_m, y_val_prob)
recall_val = recall_score(y_val_m, y_val_pred)
cm_val = confusion_matrix(y_val_m, y_val_pred)
report_val = classification_report(y_val_m, y_val_pred)

print("ROC AUC (val):", roc_auc_val)
print("Recall (val, class 1):", recall_val)
print("Confusion Matrix (val):\n", cm_val)
print("Classification Report (val):\n", report_val)

"""### Resampling: Testing the robustness of the model"""

X_test_m_drop_2.shape

X_test_m_drop_2.columns

XG_existing_thresholds = [0.0027162, 0.0028023, 0.00299955, 0.00309932, 0.00310024, 0.00329922, 0.00339915, 0.00349915,
                      0.0035, 0.0036, 0.0037, 0.0038, 0.0039, 0.0040, 0.0041, 0.0041, 0.0042, 0.0043, 0.0044, 0.0045,
                      0.0046, 0.0047, 0.0048, 0.0049, 0.0050, 0.0051, 0.0052, 0.0053, 0.0054, 0.0055, 0.0056, 0.0057,
                      0.0058, 0.0059, 0.0060, 0.0061, 0.0062, 0.0063, 0.0064, 0.0065, 0.0066, 0.0067, 0.0068, 0.0069,
                      0.0070, 0.0071, 0.0072, 0.0073, 0.0074, 0.0075, 0.0076, 0.0077, 0.0078, 0.0079, 0.0080, 0.0081,
                      0.0082, 0.0083, 0.0084, 0.0085, 0.0086, 0.0087, 0.0088, 0.0089, 0.0090, 0.0091, 0.0092, 0.0093,
                      0.0094, 0.0095, 0.0096, 0.0097, 0.0098, 0.0099, 0.010, 0.011, 0.012, 0.013, 0.014, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15,
                      0.16, 0.17, 0.18, 0.19, 0.2, 0.3, 0.35 ,0.4, 0.5, 0.6, 0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67,
                      0.68, 0.69 ,0.7, 0.8, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.961, 0.962, 0.963, 0.964, 0.965, 0.966,
                      0.967, 0.968, 0.969, 0.97, 0.98, 0.99, 0.991, 0.992, 0.993, 0.994, 0.995, 0.996, 0.997, 0.998, 0.999,
                      0.9991, 0.9992, 0.9993, 0.9994, 0.99946, 0.99947, 0.99948, 0.9995, 0.9996]

# Generating new threshold values from 0.014 to 0.1 in increments of 0.001
XG_new_thresholds = np.arange(0.014, 0.101, 0.001).tolist()

# Combine existing and new thresholds
XG_combined_thresholds = XG_existing_thresholds + XG_new_thresholds
XG_combined_thresholds.sort()

print(len(XG_combined_thresholds))

def XG_Precision_Recall_MyFunction(df, thresholds, min_positives=1):
    results = {}

    # Ensure inputs are handled correctly
    true_labels = df['True_Label'].values
    probs = df['Predicted_Probability'].values

    for threshold in thresholds:
        # Vectorized prediction based on threshold
        preds = (probs >= threshold).astype(int)

        TP = ((true_labels == 1) & (preds == 1)).sum()
        FP = ((true_labels == 0) & (preds == 1)).sum()
        FN = ((true_labels == 1) & (preds == 0)).sum()

        num_predicted_positives = TP + FP

        if num_predicted_positives >= min_positives:
            precision = TP / num_predicted_positives
            recall = TP / (TP + FN) if (TP + FN) != 0 else 0
        else:
            # Handling edge cases where no positives are predicted
            precision = 1.0
            recall = 0.0

        results[threshold] = {
            'TP': TP, 'FP': FP, 'FN': FN,
            'Precision': precision, 'Recall': recall
        }

    return results

xgb3 = xgb.XGBClassifier(
    objective='binary:logistic',
    scale_pos_weight=2,
    colsample_bytree=0.6,
    learning_rate=0.5,
    max_depth=9,
    n_estimators=501,
    reg_alpha=0.3,
    reg_lambda=10,
    use_label_encoder=False,
    eval_metric='logloss',
    min_child_weight = 1,
    gamma = 0.3,
    random_state=42
)

# --- 1. Model & Pipeline Definition ---
#xgb3 = xgb.XGBClassifier(
#    objective='binary:logistic', scale_pos_weight=2.5, colsample_bytree=0.6,
 #   learning_rate=0.5, max_depth=9, n_estimators=501, reg_alpha=0.3,
 #   reg_lambda=25, subsample=0.9, gamma=0.3, eval_metric='logloss', random_state=42
#)

# --- 1. Model & Pipeline Definition ---
base_pipeline = Pipeline(steps=[
    ('categorical_encoding', CategoricalEncoder(reduce_df=False)),
    ('classifier', xgb3)
])

# --- 2. Data Preparation ---
Features_to_Use = ['WarehouseToHome', 'SatisfactionScore', 'OrderAmountHikeFromlastYear', 'OrderCount', 'CashbackAmount']
X_train_slim = X_train_m_drop_2[Features_to_Use].copy()
X_test_slim = X_test_m_drop_2[Features_to_Use].copy()
y_test_series = pd.Series(y_test_m.ravel(), index=X_test_slim.index)

# --- 3. Training ---
# Fit ONCE on the slim training set
base_pipeline.fit(X_train_slim, y_train_m)


# Ensure validation data uses the SAME features as training
X_val_slim = X_val_m_drop_2[Features_to_Use].copy()



# --- 4. Bootstrap Robustness Loop ---
rounds = 1000
true_labels_pred_list = []
custom_threshold = 0.50

print(f"Executing {rounds} bootstrap rounds...")

for i in range(rounds):
    X_resample = X_test_slim.sample(n=len(X_test_slim), replace=True, random_state=42 + i)
    # Use the .loc accessor on the pandas Series y_test_series
    y_resample = y_test_series.loc[X_resample.index]
    # Using the existing calibrated_final here
    y_prob = base_pipeline.predict_proba(X_resample)[:, 1]


    true_labels_pred_list.append(pd.DataFrame({
        'True_Label': y_resample.values,
        'Predicted_Probability': y_prob
    }))

unique_starts = len(set([df['Predicted_Probability'].iloc[0] for df in true_labels_pred_list]))
print(f"Verification: {unique_starts}/{rounds} unique rounds generated.")

# --- 4. Precision-Recall Calculations (Per Round) ---
def create_prc_df(result_dict):
    return pd.DataFrame({
        'Threshold': list(result_dict.keys()),
        'Precision': [v['Precision'] for v in result_dict.values()],
        'Recall': [v['Recall'] for v in result_dict.values()]
    })

# Generate PR DataFrames for all 50 rounds
XG_PRC_all_rounds = [create_prc_df(XG_Precision_Recall_MyFunction(df, XG_combined_thresholds))
                     for df in true_labels_pred_list]

def sort_and_calculate_auc(recall, precision):
    """
    Sorts recall values and calculates the Area Under the Curve (AUC)
    using the trapezoidal rule.
    """
    import numpy as np
    # Recall must be monotonic (sorted) for the auc function to work correctly
    sorted_indices = np.argsort(recall)
    sorted_recall = recall[sorted_indices]
    sorted_precision = precision[sorted_indices]

    return auc(sorted_recall, sorted_precision)

# --- 5. NEW: Calculate Median Precision/Recall PER ROUND ---
round_medians = []
for i, df in enumerate(XG_PRC_all_rounds):
    m_p = df['Precision'].median()
    m_r = df['Recall'].median()
    round_medians.append({'Round': i+1, 'Median_Precision': m_p, 'Median_Recall': m_r})

df_round_medians = pd.DataFrame(round_medians)
print("\n--- Summary of Medians per Round ---")
print(df_round_medians.head())

# --- 6. Final Stability Plotting ---
plt.figure(figsize=(15, 8))
cmap = plt.get_cmap('tab20')

all_p = np.array([df['Precision'].values for df in XG_PRC_all_rounds])
all_r = np.array([df['Recall'].values for df in XG_PRC_all_rounds])

# Plot individual curves
for i in range(rounds):
    plt.plot(all_r[i], all_p[i], color=cmap(i % 20), alpha=0.1, lw=0.5)

# Calculate Global Median (The black dashed line)
median_p = np.median(all_p, axis=0)
median_r = np.median(all_r, axis=0)
XG_AP_M = sort_and_calculate_auc(median_r, median_p)

plt.plot(median_r, median_p, color='black', linestyle='--', lw=3, label=f'Global Median PRC (AUC={XG_AP_M:.2f})')

# Optional: Scatter plot the round medians to see the "center" of the fluctuations
plt.scatter(df_round_medians['Median_Recall'], df_round_medians['Median_Precision'],
            color='red', s=10, alpha=0.5, label='Individual Round Medians', zorder=5)

plt.title(f'XGBoost Stability Analysis (2025): {rounds} Bootstrap PR Curves')
plt.xlabel('Recall')
plt.ylabel('Precision')
plt.grid(True, linestyle=':', alpha=0.6)
plt.legend(loc='upper left', bbox_to_anchor=(1, 1))
plt.tight_layout()
plt.show()

"""Experimentation"""

# --- 1. Model & Pipeline Definition ---
#xgb3 = xgb.XGBClassifier(
#    objective='binary:logistic', scale_pos_weight=2.5, colsample_bytree=0.6,
 #   learning_rate=0.5, max_depth=9, n_estimators=501, reg_alpha=0.3,
 #   reg_lambda=25, subsample=0.9, gamma=0.3, eval_metric='logloss', random_state=42
#)

# --- 1. Model & Pipeline Definition ---
base_pipeline = Pipeline(steps=[
    ('categorical_encoding', CategoricalEncoder(reduce_df=False)),
    ('classifier', xgb3)
])

# --- 2. Data Preparation ---
Features_to_Use = ['WarehouseToHome', 'SatisfactionScore', 'OrderAmountHikeFromlastYear', 'OrderCount', 'CashbackAmount']
X_train_slim = X_train_m_drop_2[Features_to_Use].copy()
X_test_slim = X_test_m_drop_2[Features_to_Use].copy()
y_test_series = pd.Series(y_test_m.ravel(), index=X_test_slim.index)

# --- 3. Training ---
# Fit ONCE on the slim training set
base_pipeline.fit(X_train_slim, y_train_m)


# Ensure validation data uses the SAME features as training
X_val_slim = X_val_m_drop_2[Features_to_Use].copy()



# --- 4. Bootstrap Robustness Loop ---
rounds = 1000
true_labels_pred_list = []
custom_threshold = 0.50

print(f"Executing {rounds} bootstrap rounds...")

for i in range(rounds):
    X_resample = X_test_slim.sample(n=len(X_test_slim), replace=True, random_state=42 + i)
    # Use the .loc accessor on the pandas Series y_test_series
    y_resample = y_test_series.loc[X_resample.index]
    # Using the existing calibrated_final here
    y_prob = base_pipeline.predict_proba(X_resample)[:, 1]


    true_labels_pred_list.append(pd.DataFrame({
        'True_Label': y_resample.values,
        'Predicted_Probability': y_prob
    }))

unique_starts = len(set([df['Predicted_Probability'].iloc[0] for df in true_labels_pred_list]))
print(f"Verification: {unique_starts}/{rounds} unique rounds generated.")

# --- 4. Precision-Recall Calculations (Per Round) ---
def create_prc_df(result_dict):
    return pd.DataFrame({
        'Threshold': list(result_dict.keys()),
        'Precision': [v['Precision'] for v in result_dict.values()],
        'Recall': [v['Recall'] for v in result_dict.values()]
    })

# Generate PR DataFrames for all 50 rounds
XG_PRC_all_rounds = [create_prc_df(XG_Precision_Recall_MyFunction(df, XG_combined_thresholds))
                     for df in true_labels_pred_list]

def sort_and_calculate_auc(recall, precision):
    """
    Sorts recall values and calculates the Area Under the Curve (AUC)
    using the trapezoidal rule.
    """
    import numpy as np
    # Recall must be monotonic (sorted) for the auc function to work correctly
    sorted_indices = np.argsort(recall)
    sorted_recall = recall[sorted_indices]
    sorted_precision = precision[sorted_indices]

    return auc(sorted_recall, sorted_precision)

# --- 5. NEW: Calculate Median Precision/Recall PER ROUND ---
round_medians = []
for i, df in enumerate(XG_PRC_all_rounds):
    m_p = df['Precision'].median()
    m_r = df['Recall'].median()
    round_medians.append({'Round': i+1, 'Median_Precision': m_p, 'Median_Recall': m_r})

df_round_medians = pd.DataFrame(round_medians)
print("\n--- Summary of Medians per Round ---")
print(df_round_medians.head())

# --- 6. Final Stability Plotting ---
plt.figure(figsize=(15, 8))
cmap = plt.get_cmap('tab20')

all_p = np.array([df['Precision'].values for df in XG_PRC_all_rounds])
all_r = np.array([df['Recall'].values for df in XG_PRC_all_rounds])

# Plot individual curves
for i in range(rounds):
    plt.plot(all_r[i], all_p[i], color=cmap(i % 20), alpha=0.1, lw=0.5)

# Calculate Global Median (The black dashed line)
median_p = np.median(all_p, axis=0)
median_r = np.median(all_r, axis=0)
XG_AP_M = sort_and_calculate_auc(median_r, median_p)

plt.plot(median_r, median_p, color='black', linestyle='--', lw=3, label=f'Global Median PRC (AUC={XG_AP_M:.3f})')


# Optional: Scatter plot the round medians to see the "center" of the fluctuations
plt.scatter(df_round_medians['Median_Recall'], df_round_medians['Median_Precision'],
            color='red', s=10, alpha=0.5, label='Individual Round Medians', zorder=5)

plt.title(f'XGBoost Stability Analysis: {rounds} Bootstrap PR Curves')
plt.xlabel('Recall')
plt.ylabel('Precision')
plt.grid(True, linestyle=':', alpha=0.6)
plt.legend(loc='upper left', bbox_to_anchor=(1, 1))
plt.tight_layout()
plt.show()

"""Experimentation"""

# --- 1. Compute Percentiles for the Ribbon ---
# all_p contains 100 rounds of precision values at each threshold
lower_p = np.percentile(all_p, 2.5, axis=0)   # 2.5th percentile (Lower Bound)
upper_p = np.percentile(all_p, 97.5, axis=0)  # 97.5th percentile (Upper Bound)

# --- 2. Update the Plotting Code ---
plt.figure(figsize=(15, 8))

# Plot the Shaded Confidence Interval (Ribbon)
# This captures where 95% of your bootstrap rounds fall
plt.fill_between(median_r, lower_p, upper_p, color='lightblue', alpha=0.3, label='95% Confidence Interval')

# Plot the Global Median (The trend line)
plt.plot(median_r, median_p, color='black', linestyle='--', lw=3, label='Median PRC')

# Plot individual rounds with very low alpha to show "noise"
for i in range(len(XG_PRC_all_rounds)):
    plt.plot(all_r[i], all_p[i], color='gray', alpha=0.05, lw=0.5)

plt.title('95% Confidence Interval for Precision-Recall Stability')
plt.xlabel('Recall')
plt.ylabel('Precision')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# 1. Define the target recall level (THIS IS BASED ON THE BOOTSTRAP RESAMPLING)
target_recall = 0.72

# 2. Interpolate precision for every single bootstrap round at exactly 0.49 recall
# We use [::-1] because np.interp requires the x-axis (recall) to be increasing
precisions_at_target = []
for i in range(len(XG_PRC_all_rounds)):
    # Sort recall and precision so recall is increasing for interpolation
    sorted_idx = np.argsort(all_r[i])
    p_interp = np.interp(target_recall, all_r[i][sorted_idx], all_p[i][sorted_idx])
    precisions_at_target.append(p_interp)

# 3. Calculate the Percentiles at that specific point
lower_val = np.percentile(precisions_at_target, 2.5)
median_val = np.percentile(precisions_at_target, 50)
upper_val = np.percentile(precisions_at_target, 97.5)

print(f"At Recall {target_recall}:")
print(f"  - Lower 95% CI Precision: {lower_val:.4f}")
print(f"  - Median Precision:       {median_val:.4f}")
print(f"  - Upper 95% CI Precision: {upper_val:.4f}")

"""### Calibration and Bootstrap Resampling"""

X_val_m_drop_2.shape

# Use Beta calibration (via 'sigmoid' or custom) to wrap your existing model
calibrated_model = CalibratedClassifierCV(model_3, method='sigmoid', cv='prefit')
calibrated_model.fit(X_val_m_drop_2, y_val_m) # Use a separate validation set

# Now use calibrated_model for your bootstrap analysis

# Plot the calibration curve for the calibrated model
# n_bins=10  to balance detail and stability
disp = CalibrationDisplay.from_estimator(
    calibrated_model,
    X_test_m_drop_2,
    y_test_m,
    n_bins=3,
    name='Sigmoid Calibration'
)

plt.title('Calibration Curve: Calibrated XGBoost')
plt.show()

fig, ax = plt.subplots(figsize=(8, 6))

# Plot the uncalibrated model (model_3)
CalibrationDisplay.from_estimator(
    model_3,
    X_test_m_drop_2,
    y_test_m,
    n_bins=10,
    name='Uncalibrated XGBoost',
    ax=ax
)

# Plot the calibrated model
CalibrationDisplay.from_estimator(
    calibrated_model,
    X_test_m_drop_2,
    y_test_m,
    n_bins=3,
    name='Calibrated (Sigmoid)',
    ax=ax
)

ax.set_title('Calibration Comparison')
plt.legend()
plt.show()

# 1. Get probabilities for the test set
raw_probs = model_3.predict_proba(X_test_m_drop_2)[:, 1]
cal_probs = calibrated_model.predict_proba(X_test_m_drop_2)[:, 1]  # USING TEST SET ?????

# 2. Calculate the Brier Score
brier_raw = brier_score_loss(y_test_m, raw_probs)
brier_cal = brier_score_loss(y_test_m, cal_probs)

print(f"Brier Score (Raw Model):        {brier_raw:.4f}")
print(f"Brier Score (Calibrated Model):   {brier_cal:.4f}")

# Calculate improvement percentage
improvement = (brier_raw - brier_cal) / brier_raw * 100
print(f"Improvement: {improvement:.2f}%")

XG_existing_thresholds = [0.0027162, 0.0028023, 0.00299955, 0.00309932, 0.00310024, 0.00329922, 0.00339915, 0.00349915,
                      0.0035, 0.0036, 0.0037, 0.0038, 0.0039, 0.0040, 0.0041, 0.0041, 0.0042, 0.0043, 0.0044, 0.0045,
                      0.0046, 0.0047, 0.0048, 0.0049, 0.0050, 0.0051, 0.0052, 0.0053, 0.0054, 0.0055, 0.0056, 0.0057,
                      0.0058, 0.0059, 0.0060, 0.0061, 0.0062, 0.0063, 0.0064, 0.0065, 0.0066, 0.0067, 0.0068, 0.0069,
                      0.0070, 0.0071, 0.0072, 0.0073, 0.0074, 0.0075, 0.0076, 0.0077, 0.0078, 0.0079, 0.0080, 0.0081,
                      0.0082, 0.0083, 0.0084, 0.0085, 0.0086, 0.0087, 0.0088, 0.0089, 0.0090, 0.0091, 0.0092, 0.0093,
                      0.0094, 0.0095, 0.0096, 0.0097, 0.0098, 0.0099, 0.010, 0.011, 0.012, 0.013, 0.014, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15,
                      0.16, 0.17, 0.18, 0.19, 0.2, 0.3, 0.35 ,0.4, 0.5, 0.6, 0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67,
                      0.68, 0.69 ,0.7, 0.8, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.961, 0.962, 0.963, 0.964, 0.965, 0.966,
                      0.967, 0.968, 0.969, 0.97, 0.98, 0.99, 0.991, 0.992, 0.993, 0.994, 0.995, 0.996, 0.997, 0.998, 0.999,
                      0.9991, 0.9992, 0.9993, 0.9994, 0.99946, 0.99947, 0.99948, 0.9995, 0.9996]

# Generating new threshold values from 0.014 to 0.1 in increments of 0.001
XG_new_thresholds = np.arange(0.014, 0.101, 0.001).tolist()

# Combine existing and new thresholds
XG_combined_thresholds = XG_existing_thresholds + XG_new_thresholds
XG_combined_thresholds.sort()

print(len(XG_combined_thresholds))

def XG_Precision_Recall_MyFunction(df, thresholds, min_positives=1):
    results = {}

    # Ensure inputs are handled correctly
    true_labels = df['True_Label'].values
    probs = df['Predicted_Probability'].values

    for threshold in thresholds:
        # Vectorized prediction based on threshold
        preds = (probs >= threshold).astype(int)

        TP = ((true_labels == 1) & (preds == 1)).sum()
        FP = ((true_labels == 0) & (preds == 1)).sum()
        FN = ((true_labels == 1) & (preds == 0)).sum()

        num_predicted_positives = TP + FP

        if num_predicted_positives >= min_positives:
            precision = TP / num_predicted_positives
            recall = TP / (TP + FN) if (TP + FN) != 0 else 0
        else:
            # Handling edge cases where no positives are predicted
            precision = 1.0
            recall = 0.0

        results[threshold] = {
            'TP': TP, 'FP': FP, 'FN': FN,
            'Precision': precision, 'Recall': recall
        }

    return results

xgb3 = xgb.XGBClassifier(
    objective='binary:logistic',
    scale_pos_weight=2,
    colsample_bytree=0.6,
    learning_rate=0.5,
    max_depth=9,
    n_estimators=501,
    reg_alpha=0.3,
    reg_lambda=10,
    use_label_encoder=False,
    eval_metric='logloss',
    min_child_weight = 1,
    gamma = 0.3,
    random_state=42
)

# --- 1. Model & Pipeline Definition ---
#xgb3 = xgb.XGBClassifier(
#    objective='binary:logistic', scale_pos_weight=2.5, colsample_bytree=0.6,
 #   learning_rate=0.5, max_depth=9, n_estimators=501, reg_alpha=0.3,
 #   reg_lambda=25, subsample=0.9, gamma=0.3, eval_metric='logloss', random_state=42
#)

# --- 1. Model & Pipeline Definition ---
base_pipeline = Pipeline(steps=[
    ('categorical_encoding', CategoricalEncoder(reduce_df=False)),
    ('classifier', xgb3)
])

# --- 2. Data Preparation ---
Features_to_Use = ['WarehouseToHome', 'SatisfactionScore', 'OrderAmountHikeFromlastYear', 'OrderCount', 'CashbackAmount']
X_train_slim = X_train_m_drop_2[Features_to_Use].copy()
X_test_slim = X_test_m_drop_2[Features_to_Use].copy()
y_test_series = pd.Series(y_test_m.ravel(), index=X_test_slim.index)

# --- 3. Training ---
# Fit ONCE on the slim training set
base_pipeline.fit(X_train_slim, y_train_m)

# --- 4. Calibration ---
# Wrap the FITTED base_pipeline. Use cv='prefit'
calibrated_final = CalibratedClassifierCV(base_pipeline, cv='prefit', method='sigmoid')

# Ensure validation data uses the SAME features as training
X_val_slim = X_val_m_drop_2[Features_to_Use].copy()
calibrated_final.fit(X_val_slim, y_val_m)

# --- 5. Ready for Bootstrap ---
#  use calibrated_final.predict_proba() in your loop

# --- 3. Bootstrap Robustness Loop ---
rounds = 1000
true_labels_pred_list = []
custom_threshold = 0.50

print(f"Executing {rounds} bootstrap rounds...")

for i in range(rounds):
    X_resample = X_test_slim.sample(n=len(X_test_slim), replace=True, random_state=42 + i)
    # Use the .loc accessor on the pandas Series y_test_series
    y_resample = y_test_series.loc[X_resample.index]
    # Using the existing calibrated_final here
    y_prob = calibrated_final.predict_proba(X_resample)[:, 1]


    true_labels_pred_list.append(pd.DataFrame({
        'True_Label': y_resample.values,
        'Predicted_Probability': y_prob
    }))

unique_starts = len(set([df['Predicted_Probability'].iloc[0] for df in true_labels_pred_list]))
print(f"Verification: {unique_starts}/{rounds} unique rounds generated.")

# --- 4. Precision-Recall Calculations (Per Round) ---
def create_prc_df(result_dict):
    return pd.DataFrame({
        'Threshold': list(result_dict.keys()),
        'Precision': [v['Precision'] for v in result_dict.values()],
        'Recall': [v['Recall'] for v in result_dict.values()]
    })

# Generate PR DataFrames for all 50 rounds
XG_PRC_all_rounds = [create_prc_df(XG_Precision_Recall_MyFunction(df, XG_combined_thresholds))
                     for df in true_labels_pred_list]

def sort_and_calculate_auc(recall, precision):
    """
    Sorts recall values and calculates the Area Under the Curve (AUC)
    using the trapezoidal rule.
    """
    import numpy as np
    # Recall must be monotonic (sorted) for the auc function to work correctly
    sorted_indices = np.argsort(recall)
    sorted_recall = recall[sorted_indices]
    sorted_precision = precision[sorted_indices]

    return auc(sorted_recall, sorted_precision)

# --- 5. NEW: Calculate Median Precision/Recall PER ROUND ---
round_medians = []
for i, df in enumerate(XG_PRC_all_rounds):
    m_p = df['Precision'].median()
    m_r = df['Recall'].median()
    round_medians.append({'Round': i+1, 'Median_Precision': m_p, 'Median_Recall': m_r})

df_round_medians = pd.DataFrame(round_medians)
print("\n--- Summary of Medians per Round ---")
print(df_round_medians.head())

# --- 6. Final Stability Plotting ---
plt.figure(figsize=(15, 8))
cmap = plt.get_cmap('tab20')

all_p = np.array([df['Precision'].values for df in XG_PRC_all_rounds])
all_r = np.array([df['Recall'].values for df in XG_PRC_all_rounds])

# Plot individual curves
for i in range(rounds):
    plt.plot(all_r[i], all_p[i], color=cmap(i % 20), alpha=0.1, lw=0.5)

# Calculate Global Median (The black dashed line)
median_p = np.median(all_p, axis=0)
median_r = np.median(all_r, axis=0)
XG_AP_M = sort_and_calculate_auc(median_r, median_p)

plt.plot(median_r, median_p, color='black', linestyle='--', lw=3, label=f'Global Median PRC (AUC={XG_AP_M:.3f})')

# Optional: Scatter plot the round medians to see the "center" of the fluctuations
plt.scatter(df_round_medians['Median_Recall'], df_round_medians['Median_Precision'],
            color='red', s=10, alpha=0.5, label='Individual Round Medians', zorder=5)

plt.title(f'XGBoost Stability Analysis (2025): {rounds} Bootstrap PR Curves')
plt.xlabel('Recall')
plt.ylabel('Precision')
plt.grid(True, linestyle=':', alpha=0.6)
plt.legend(loc='upper left', bbox_to_anchor=(1, 1))
plt.tight_layout()
plt.show()

# --- 1. Compute Percentiles for the Ribbon ---
# all_p contains 100 rounds of precision values at each threshold
lower_p = np.percentile(all_p, 2.5, axis=0)   # 2.5th percentile (Lower Bound)
upper_p = np.percentile(all_p, 97.5, axis=0)  # 97.5th percentile (Upper Bound)

# --- 2. Update the Plotting Code ---
plt.figure(figsize=(15, 8))

# Plot the Shaded Confidence Interval (Ribbon)
# This captures where 95% of your bootstrap rounds fall
plt.fill_between(median_r, lower_p, upper_p, color='lightblue', alpha=0.3, label='95% Confidence Interval')

# Plot the Global Median (The trend line)
plt.plot(median_r, median_p, color='black', linestyle='--', lw=3,
         label=f'Median PRC (AUC = {XG_AP_M:.3f})')


# Plot individual rounds with very low alpha to show "noise"
for i in range(len(XG_PRC_all_rounds)):
    plt.plot(all_r[i], all_p[i], color='gray', alpha=0.05, lw=0.5)

plt.title('95% Confidence Interval for Precision-Recall Stability')
plt.xlabel('Recall')
plt.ylabel('Precision')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# 1. Define the target recall level (THIS IS BASED ON THE BOOTSTRAP RESAMPLING)
target_recall = 0.72

# 2. Interpolate precision for every single bootstrap round at exactly 0.49 recall
# We use [::-1] because np.interp requires the x-axis (recall) to be increasing
precisions_at_target = []
for i in range(len(XG_PRC_all_rounds)):
    # Sort recall and precision so recall is increasing for interpolation
    sorted_idx = np.argsort(all_r[i])
    p_interp = np.interp(target_recall, all_r[i][sorted_idx], all_p[i][sorted_idx])
    precisions_at_target.append(p_interp)

# 3. Calculate the Percentiles at that specific point
lower_val = np.percentile(precisions_at_target, 2.5)
median_val = np.percentile(precisions_at_target, 50)
upper_val = np.percentile(precisions_at_target, 97.5)

print(f"At Recall {target_recall}:")
print(f"  - Lower 95% CI Precision: {lower_val:.4f}")
print(f"  - Median Precision:       {median_val:.4f}")
print(f"  - Upper 95% CI Precision: {upper_val:.4f}")

"""In this section, the model has been calibrated, and 95% confidence intervals were generated based on the test data. We performed 1,000 bootstrap resamples with replacement; this process simulates variations in data difficulty to ensure the model's performance is robust. For example, at a recall of 0.7, the median precision is 0.43. We are 95% confident that the true precision of the model at this recall level lies between 0.35 and 0.51. This range provides valuable insight for stakeholders by quantifying the expected stability of the model in production

### SHAP
"""

import shap

shap.initjs()

# Now fit your calibrated model
calibrated_model.fit(X_train_m_drop_2, y_train_m.ravel())

# 1. Create a wrapper that handles the transformation AND the prediction
# This allows you to pass raw DataFrames to SHAP
def full_pipeline_predict(data):
    # If your data comes as a NumPy array (SHAP often does this internally),
    # convert it back to a DataFrame so the pipeline can read column names.
    if not isinstance(data, pd.DataFrame):
        data = pd.DataFrame(data, columns=X_train_m_drop_2.columns)

    # Use the estimator (the pipeline) stored inside the calibrated model
    return calibrated_model.predict_proba(data)[:, 1]

# 2. Initialize the Explainer using the raw (untransformed) background data
# We now use X_train_slim directly
explainer_cal = shap.Explainer(full_pipeline_predict, X_train_slim.iloc[:100])

# 3. Calculate SHAP values for Row 7 using raw data
# No manual encoder.transform() needed here
shap_values_cal = explainer_cal(X_test_slim.iloc[7:8])

# 4. Visualize
# The feature names in the plot will now be your original, human-readable names
shap.plots.waterfall(shap_values_cal[0])

# 1. Define the row index
row_num = 9

# 2. Get the explanation for just that one row
# Pass the raw DataFrame slice (explainer handles the internal pipeline logic)
shap_values_row_7 = explainer_cal(X_train_m_drop_2.iloc[row_num : row_num + 1])

# 3. Visualize using the waterfall plot
# Use [0] because the explanation object above is a list-like of one row
shap.plots.waterfall(shap_values_row_7[0])

"""Row number 7 or 8 have an overall probability of churn of 0.20. all of the features decrease the probability of churn."""

# 1. Generate explanations for a representative sample of 500 rows
# This provides a stable global view without the long wait time
shap_values_global = explainer_cal(X_train_m_drop_2.sample(500, random_state=42))

# 2. Plot the beeswarm
shap.plots.beeswarm(shap_values_global)



import shap
shap.initjs() # Required for interactive JavaScript plots

# Simply pass the row's explanation object
# Modern force_plot automatically detects base_values and feature data
shap.plots.force(shap_values_cal[0], link='logit')

# Extract values from your row 7 explanation object
base_val = shap_values_cal[0].base_values
shap_vals_array = shap_values_cal[0].values
feature_data = X_test_slim.iloc[7]

shap.force_plot(
    base_val,
    shap_vals_array,
    feature_data,
    link='logit'
)

## Check probability predictions through the model
pred_probs = calibrated_model.predict_proba(X_train_m_drop_2)[:,1]
pred_probs[row_num]

# 1. Generate an Explanation object instead of a raw array
# This is faster and bundles all metadata automatically
explanation = explainer_cal(X_train_m_drop_2)

# 2. Use the modern beeswarm plot
# This is the modern equivalent of the 'summary_plot'
shap.plots.beeswarm(explanation)

# Shows a simple horizontal bar chart of global importance
shap.summary_plot(shap_values, X_train_m_drop_2, plot_type="bar")

# Beeswarm is the best way to see Class 1 vs 0 direction globally
shap.summary_plot(shap_values, X_transformed, plot_type="dot", alpha=0.5)

"""**Summary Table for Your Plot**:
- Location	Dot Color	Meaning for Churn
- Right Side	Blue	Low values of this feature increase churn risk.
- Right Side	Red	High values of this feature increase churn risk.
- Left Side	Blue	Low values of this feature decrease churn risk.
- Left Side	Red	High values of this feature decrease churn risk.

- Cash back Amount: This feature has the highest SHAP absolute value, meaning it is the most impactful feature for the model's predictions.
The blue area towards the right indicates that low values for cash back amount increase the churn risk.
The red area (which would be towards the left) on the other hand decreases the churn risk, indicating high cash back amounts help retain customers.

- Satiscation score has the lowest SHAP absolute value. The red are towards right side increases the churn risk.
"""

# 1. Get the SHAP values for your positive class (Class 1)
# Using the explanation object from your existing explainer_cal
explanation = explainer_cal(X_train_m_drop_2)
shap_values_class1 = explanation.values

# 2. Create a list of two arrays to trick SHAP into "multiclass" mode
# Class 0's impact is just the inverse of Class 1
shap_values_multiclass = [shap_values_class1 * -1, shap_values_class1]

# 3. Generate the stacked bar chart
# This will show two colors (one for each class) for every feature
shap.summary_plot(
    shap_values_multiclass,
    X_train_m_drop_2,
    plot_type="bar",
    class_names=["Class 0", "Class 1"]
)

"""### Load saved model and make predictions on unseen/future data

Here, we'll use df_test as the unseen, future data
"""

import joblib
joblib.dump(calibrated_model, 'calibrated_model_1.sav')

import joblib
import xgboost as xgb
import sklearn
import numpy as np

# Verify the versions being baked into the file
print(f"XGBoost Version: {xgb.__version__}")
print(f"Scikit-Learn Version: {sklearn.__version__}")
print(f"NumPy Version: {np.__version__}")

# Check the type of the object you just saved
print(f"Model Object Type: {type(calibrated_model)}")

# Create a new steps list excluding the encoder
new_steps = [step for step in calibrated_model.estimator.steps if step[0] != 'categorical_encoding']

# Update the pipeline with the filtered steps
calibrated_model.estimator.steps = new_steps

calibrated_model

import xgboost as xgb
import sklearn
import pandas as pd
import joblib

print(f"XGBoost Version: {xgb.__version__}")
print(f"Scikit-Learn Version: {sklearn.__version__}")
print(f"Pandas Version: {pd.__version__}")

import joblib
from google.colab import files

# 1. Save the model to the Colab cloud environment
joblib.dump(calibrated_model, 'calibrated_model_1.sav')

# 2. Download the file from the cloud to your local machine
files.download('calibrated_model_1.sav')

import xgboost
print(xgboost.__version__)

## Save model object
joblib.dump(model, 'final_churn_model_f1_0_45.sav')

base_model = calibrated_model.estimator.named_steps['classifier']
base_model.save_model('model_native.json')
from google.colab import files
files.download('model_native.json')

import joblib

# 'model' is your trained object
#joblib.dump(calibrated_model, 'my_model.joblib')

#from google.colab import files
#files.download('calibrated_model.joblib')



from sklearn.metrics import f1_score

# Test different thresholds to find the new peak F1
thresholds = np.linspace(0.1, 0.9, 50)
f1_scores = [f1_score(y_test_m, cal_probs > t) for t in thresholds]
best_t = thresholds[np.argmax(f1_scores)]

print(f"Optimal threshold for Calibrated Model: {best_t:.2f}")

import joblib

# 1. Bypass the encoder using the variable name 'model'
# If you named your variable something else (like 'fianl_model'), use that name here instead!
model.set_params(estimator__categorical_encoding='passthrough')

# 2. Save it with compression to ensure compatibility with your local Mac/AWS
joblib.dump(model, 'final_model_v2.joblib', compress=3)

# 3. Download the new file
from google.colab import files
files.download('final_model_v2.joblib')

## Load model object
#model = joblib.load('final_churn_model_f1_0_45.sav')

X_test_m_drop_2.shape
y_test_m.shape

## Predict target probabilities
test_probs = calibrated_model.predict_proba(X_test_m_drop_2)[:,1]

## Predict target values on test data
test_preds = np.where(test_probs > 0.38, 1, 0) # Flexibility to tweak the probability threshold
#test_preds = model.predict(X_test)

# Use keyword arguments x and y
sns.boxplot(x=y_test_m.ravel(), y=test_probs)

## Test set metrics
roc_auc_score(y_test_m, test_preds)
recall_score(y_test_m, test_preds)
confusion_matrix(y_test_m, test_preds)
print(classification_report(y_test_m, test_preds))

"""A Precision of 0.69 means:
If the model points at 100 people and says "These people are going to churn,"
69 of them actually will churn (True Positives).
31 of them actually will stay (False Positives).
The difference between your Recall and Precision:
Recall (0.72) is about the Source:
"Out of all the people who really churned, how many did I catch?"
Result: You caught 72% of the total "churners."
Precision (0.69) is about the Quality of your List:
"Out of all the people I flagged, how many were actually churners?"

At Recall 0.72:
  - Lower 95% CI Precision: 0.3541
  - Median Precision:       0.4290
  - Upper 95% CI Precision: 0.5076
"""

## Predict target probabilities
test_probs = calibrated_model.predict_proba(X_test_m_drop_2)[:,1]

## Predict target values on test data
test_preds = np.where(test_probs > 0.13, 1, 0) # Flexibility to tweak the probability threshold
#test_preds = model.predict(X_test)

## Test set metrics
roc_auc_score(y_test_m, test_preds)
recall_score(y_test_m, test_preds)
confusion_matrix(y_test_m, test_preds)
print(classification_report(y_test_m, test_preds))

"""At Recall 0.82:
  - Lower 95% CI Precision: 0.3122
  - Median Precision:       0.4029
  - Upper 95% CI Precision: 0.4832
"""

## Adding predictions and their probabilities in the original test dataframe
test = X_test_m_drop_2.copy()
test['predictions'] = test_preds
test['pred_probabilities'] = test_probs
test['True_Label'] = y_test_m

test.sample(10)

"""#### Creating a list of customers who are the most likely to churn

Listing customers who have a churn probability higher than 38%. These are the ones who can be targeted immediately. Sensitivity over Caution: By lowering the threshold from 0.5 to 0.38, you are telling the model: "I would rather accidentally flag a loyal customer (False Positive) than miss someone who is actually leaving (False Negative)".
"""

high_churn_list = test[test.pred_probabilities > 0.38].sort_values(by = ['pred_probabilities'], ascending = False
                                                                 ).reset_index().drop(columns = ['index', 'True_Label', 'predictions'], axis = 1)

high_churn_list.shape
high_churn_list.sample(20)

# high_churn_list.to_csv('high_churn_list.csv', index = False)

"""#### Feature-based user segments from the above list
Based on business requirements, a prioritization matrix can be defined, wherein certain segments of customers are targeted first. These segments can be defined based on insights through data or the business teams' requirements.
E.g. Males who are an ActiveMember, have a CreditCard and are from Germany can be prioritized first because the business potentially sees the max. ROI from them

### Ending notes

#### Note on common issues with a model in production

- Data drift / Covariate shift

- Importance of incremental training

- Ensure parity between training and testing environments (model and library versions etc.)

- Tracking core business metrics

- Creation and monitoring of metrics of specific user segments

- Highlight impact to business folks : Through visualizations, Model can potentially reduce the Churn rate by 30-40% etc.

#### Future steps

 - The model can be expanded to predict when will a customer churn. This will further help sales/customer service teams to reduce churn rate by targeting the right customers at the right time
"""