# By looking at the below info, we see there are a lot of missing CustomerID. As we analyze customers, the rows with missing CustomerID 
# need to be dropped. There are also some NA in Description column but we ignore them as we will drop the column anyways. We also drop the duplicates.
print(df.isna().sum())
df.dropna(subset=['CustomerID'], inplace=True)
df.drop_duplicates(inplace=True)
df.reset_index(drop=True,inplace=True)

# Outliers and abnormal data check.
print('\n*** Describe Quantity & UnitPrice ***')
print(df[['Quantity','UnitPrice']].describe())
# The above description tells us that there is at least one negative number in Quantity column which shouldn't be. Let's check further.
print('\n\t\t\t\t\t\t*** Smallest 5 rows in Quantity ***')
print(df[df['Quantity']<=0].sort_values(by='Quantity').head(5))
print('\n\t\t\t\t\t\t*** Biggest 5 rows in Quantity ***')
print(df[df['Quantity']>0].sort_values(by='Quantity',ascending=False).head(5))
# From the printed info, we can tell that for rows 401131 & 401132, it might be a mistake and it was adjusted in the same day. Same for row 37511 & 37516 
df.drop([401131,401132,37511,37516],inplace=True)
# Drop the rows with negative & 0 Quantity
df = df[df.Quantity>0]

# Drop the columns not used in this case
df.drop(columns=['StockCode','Description','Country'],inplace=True)
# Create a new Total column and drop the two columns used to calculate the Total
df['Total'] = df.Quantity * df.UnitPrice
df.drop(columns=['Quantity','UnitPrice'],inplace=True)
# By looking at the data types, we can tell that we might need to change the data types of 'InvoiceDate' and 'CustomerID' columns.
print('\n*** Data types check ***')
print(df.dtypes)
# Change Invoicedate column data type to datetime then get the date part only
df.InvoiceDate = pd.to_datetime(df.InvoiceDate).apply(lambda x: x.date())
# Check if all CustomerIDs are int data type. If it is tue, change the data type.
print('\n*** Int, Float type check ***')
print('Are all the CustomerIDs Int data type?  '+str(all(df.CustomerID%1==0)))
df.CustomerID = df.CustomerID.astype('int')
# Define that if customers don't have any transactions in 180 days, we classify them as churned customer and don't include them in the analysis.
LastDate = df['InvoiceDate'].max()
dfCut = df[df['InvoiceDate'].apply(lambda x: (LastDate-x).days<180)]
# Check the final dataframe.
# print(tabulate(dfCut.head(2), headers='keys', tablefmt='html', numalign="center"))

# Recency
RFM = dfCut[['InvoiceDate','CustomerID']].groupby('CustomerID',as_index=False)['InvoiceDate'].apply(lambda x: (LastDate-x.max()).days)
# Frequency 
Frequency = dfCut[['InvoiceNo','CustomerID']].groupby('CustomerID',as_index=False).size()
# Monetary
Monetary = df[['CustomerID','Total']].groupby('CustomerID',as_index=False).sum()
RFM = RFM.merge(Frequency).merge(Monetary)
RFM.rename(columns={'InvoiceDate':'Recency','size':'Frequency','Total':'Monetary'},inplace=True)
print(RFM.describe())

# Check skewness
RFMColumns = ['Recency','Frequency','Monetary']
RFMLog = RFM[RFMColumns]
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for i, ax in zip(RFMLog.columns.to_list(),axes):    
        sns.histplot(RFMLog[i], ax=ax, kde=True, stat="density", linewidth=0)  
plt.show() 

# Remove skewness
RFMLogCopy = RFMLog.copy()   # To avoid Handling SettingWithCopyWarning, use a copy of the dataframe
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for i, ax in zip(RFMLogCopy.columns.to_list(),axes):
    small_const = 1e-10
    RFMLogCopy[i] = pd.Series(np.log(RFMLogCopy[i]+small_const))
    sns.histplot(RFMLogCopy[i], ax=ax, kde=True, stat="density", linewidth=0)  
plt.show() 

# Standardize - Combining centering and scaling using a scaler
scaler = StandardScaler()
RFMSt = pd.DataFrame(scaler.fit_transform(RFMLog))
RFMSt.columns = RFMLog.columns.to_list()
print(RFMSt.describe().round(2))

# Use Elbow method,silhouette_score, calinski_harabasz_score to determine the optimal K.
Scores = ['Inertia','SHScores','CHScores']
Inertia = []
SHScores = []
CHScores = []
for k in range(2, 11):
    model = KMeans(n_clusters=k, random_state=1)
    model.fit(RFMSt)
    pre = model.predict(RFMSt)    
    Inertia.append(model.inertia_)
    SHScores.append(silhouette_score(RFMSt,pre))
    CHScores.append(calinski_harabasz_score(RFMSt,pre))
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for Score, ax in zip(Scores,axes):    
        sns.pointplot(x=list(range(2,11)), y=globals()[Score], ax=ax)
        ax.set_xlabel('K')
        ax.set_title(f'{Score}')        
plt.show() 

# Use scatterplot3D to check the clusters
model = KMeans(n_clusters=4,random_state=1231).fit(RFMSt)
RFM['Clusters'] = model.labels_
print(RFM.head())
categoryList = []
categoryDfList =[]
[categoryList.append('category'+str(i)) for i in range(4)]
for category in categoryList:
    globals()[category+'Df'] = RFM.loc[RFM['Clusters']==int(category[-1])][RFMColumns]        
    categoryDfList.append(globals()[category+'Df'])
plt.figure(figsize=(15,15))
ax = plt.subplot(121, projection='3d')
colorList = ['y','r','g','b']
for i in range(len(categoryDfList)):    
    ax.scatter(categoryDfList[i][[RFMColumns[0]]],categoryDfList[i][[RFMColumns[1]]],categoryDfList[i][[RFMColumns[2]]],c=colorList[i])
plt.gca().set(title = '3D scatterplot for K=4 ', xlabel='Recency', ylabel='Frequency',zlabel='Monetary' )
# Format number in thousands using k.
new_ticks = [f'{int(tick/1000)}k' for tick in ax.get_zticks()]
ax.set_zticks(ax.get_zticks())  # Explicitly set the ticks to avoid warning
ax.set_zticklabels(new_ticks)
# Set X,Y,Z Axes box aspect
ax.set_box_aspect([1,1,1.11])  
plt.show()

# Get the statistics for the clusters
ClustersInfo = RFM.groupby('Clusters').agg({'Recency':'mean','Frequency':'mean','Monetary':'mean','CustomerID':'sum'})
print(ClustersInfo)