import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import r2_score, mean_absolute_error, mean_squared_error
from scipy.stats import pearsonr

df = pd.read_csv('garments_worker_productivity.csv')
print(df.shape)
df.head()

len(df[df['over_time']==0])

df.info()

df.department.value_counts()

df['department'] = df['department'].replace(['sweing'], 'sewing')
df['department'] = df['department'].replace(['finishing '], 'finishing')
df['department'] = df.department.astype('category')

df.department.value_counts()

# variable that store duplicates exactly the same as another row
dups = df.duplicated()
# count duplicates
dups.value_counts()

df.isna().sum()

finishing = df[df['department']=='finishing']
sewing = df[df['department']=='sewing']

print(finishing.shape)
sewing.shape

finishing.wip.isna().sum()

# Converting wip NaN values to zero
df.wip.fillna(0, inplace=True)

# Code preparation for plotting.
# Mapping departent to 0 and 1 as integer 
df['dept_mapped'] = df.department.map({
                                'finishing':0,
                                'sewing':1})
# convert to int dtype
df.dept_mapped  = df.dept_mapped.astype('int64')

# plt.figure(figsize=(1, 1)) 
g = sns.lmplot(x='dept_mapped', y='wip', 
               data=df, 
               x_jitter=0.15, 
               y_jitter=0.15, 
               fit_reg=False, 
               scatter_kws={'alpha':0.2})

g.set(xticks=range(2))
g.set_xticklabels(['finishing', 'sewing'])
plt.show()

# date format is MM-DD-YYYY
df['date'] = pd.to_datetime(df.date, format="%m/%d/%Y")

# add month
df['month'] = df['date'].dt.month

# Adding `work_week` columns. It's easier to analyze production data per work week.
# create a list of our conditions
conditions = [
    
    #January
    (df['month'] == 1) & (df['quarter'] == 'Quarter1'), #ww1
    (df['month'] == 1) & (df['quarter'] == 'Quarter2'), #ww2
    (df['month'] == 1) & (df['quarter'] == 'Quarter3'), #ww3
    (df['month'] == 1) & (df['quarter'] == 'Quarter4'), #ww4
    (df['month'] == 1) & (df['quarter'] == 'Quarter5'), #ww4
    
    #February
    (df['month'] == 2) & (df['quarter'] == 'Quarter1'), #ww5
    (df['month'] == 2) & (df['quarter'] == 'Quarter2'), #ww6
    (df['month'] == 2) & (df['quarter'] == 'Quarter3'), #ww7
    (df['month'] == 2) & (df['quarter'] == 'Quarter4'), #ww8
    (df['month'] == 2) & (df['quarter'] == 'Quarter5'), #ww8
    
    #March
    (df['month'] == 3) & (df['quarter'] == 'Quarter1'), #ww9
    (df['month'] == 3) & (df['quarter'] == 'Quarter2'), #ww_10
    (df['month'] == 3) & (df['quarter'] == 'Quarter3'), #ww_11
    (df['month'] == 3) & (df['quarter'] == 'Quarter4'), #ww_12
    (df['month'] == 3) & (df['quarter'] == 'Quarter5'), #ww_12
    ]

values = ['ww1','ww2','ww3','ww4','ww4',
          'ww5','ww6','ww7','ww8','ww8',
          'ww9','ww_10','ww_11','ww_12','ww_12']

df['work_week'] = np.select(conditions, values)

quarter_order = ['Quarter1', 'Quarter2', 'Quarter3', 'Quarter4', 'Quarter5' ]
df.quarter = pd.Categorical(df.quarter, quarter_order, ordered=True)
# verify quarter order
df.quarter.unique()

df.day.unique()

day_order = ['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Saturday', 'Sunday' ]
df.day = pd.Categorical(df.day, day_order, ordered=True)
# verify day order
df.day.unique()

# how many zeros are there in each columns?
# (df == 0).sum(axis=0) # uncomment to show result

print('Proportion of zeros: ')
(df == 0).sum(axis=0)  / len(df) * 100

# we will not include date in this visualization
df.drop(columns=['date']).hist(figsize=(10, 7))
plt.tight_layout()
plt.show()

# select numeric columns only
df_numeric = df.select_dtypes(include=np.number).columns.tolist()
for i in df_numeric:
    print(i + ' ==> ' + str(round(df[i].skew(),2)))

df.select_dtypes(include=np.number).kurt()

team = df.groupby('team').mean() # groupby then get the mean
team['diff_in_productivity'] = team.actual_productivity - team.targeted_productivity

# limit to 2 decimal places
team = team.round(2)

team[['targeted_productivity', 'actual_productivity', 'diff_in_productivity']]

team.actual_productivity.sort_values(ascending=False).head(3)

#
x = team.index
y = team.actual_productivity

plt.figure(figsize=(8,5))
ax1 = plt.subplot()
sns.barplot(x=x, y=y, palette = 'deep')
plt.title('Productivity by Team')

for i in ax1.containers:
    ax1.bar_label(i,)

plt.tight_layout()    
plt.show()

team_dept = df.groupby(['team','department']).mean() 
team_dept['diff_productivity'] = team_dept.actual_productivity - team_dept.targeted_productivity
team_dept = team_dept.round(2)

# Show first 3 teams
team_dept[['targeted_productivity', 'actual_productivity', 'diff_productivity']].head(6)

# code preparation for visualization
team_act_prod = []
for i in team_dept.actual_productivity:
    team_act_prod.append(i)

sewing = team_act_prod[1::2]
finishing= team_act_prod[0::2]

# team_dept is a 2d list
# example code for accessing 2d indeces
# team 2 finishing dept would be:
# team_dept.actual_productivity[2][0]

#
fig, ax = plt.subplots(figsize=(12,6))

x = np.arange(1,13)
width = 0.35  # the width of the bars

rects1 = ax.bar(x - width/2, finishing, width, label= 'finishing')
rects2 = ax.bar(x + width/2, sewing, width, label='sewing')
ax.set_xticks(x)
ax.set_xlabel('Teams')
ax.set_ylabel('Actual Productivity')

def autolabel(rects):
    """Attach a text label above each bar in *rects*, displaying its height."""
    for rect in rects:
        height = rect.get_height()
        ax.annotate('{}'.format(height),
                    xy=(rect.get_x() + rect.get_width() / 2, height),
                    xytext=(0, 3),  # 3 points vertical offset
                    textcoords="offset points",
                    ha='center', va='bottom')

plt.title('Actual Productivity by Team Department')
autolabel(rects1)
autolabel(rects2)
plt.tight_layout()
plt.legend(loc='upper center')
plt.show()

from scipy.stats import ttest_ind

team1 = df[df.team == 1]
team3 = df[df.team == 3]

plt.hist(team1.actual_productivity, alpha=0.5, label='team1')
plt.hist(team3.actual_productivity,alpha=0.5, label='team3')
plt.legend()
plt.show()

# Check if STD is equal
# a ratio between 0.9 and 1.1 should suffice
ratio = np.std(team1.actual_productivity) / np.std(team3.actual_productivity)
ratio

#run the t-test here:
tstat, pval =ttest_ind(team1.actual_productivity, team3.actual_productivity)
pval

print('Average Productivity per work_week')
target_actual = df.groupby('work_week').mean()
target_actual[['targeted_productivity', 'actual_productivity']]

# setting variables
# target productivity average per workweek
tp_avg = df.groupby('work_week')['targeted_productivity'].mean()
# actual productivity average per workweek
ap_avg = df.groupby('work_week')['actual_productivity'].mean()

plt.figure(figsize=(10,5))

plt.plot(tp_avg.index, tp_avg, marker='o')
plt.plot(ap_avg.index, ap_avg, marker='o')

plt.xlabel('work weeks')
plt.ylabel('productivity score')
plt.legend(['targeted productivity', 'actual productivity'])

plt.show()

# Checking the reason behind the crash
idle = df.groupby('work_week').sum()
idle['idle_time']

def correlation_checker(col):
    plt.figure(figsize=(6, 6)) 
    corr_matrix = df.corr()
 
    corr_target = corr_matrix[[col]].drop(labels=[col])
    corr_target = corr_target.sort_values(by=[col], ascending=False)

    sns.heatmap(corr_target, annot=True, fmt='.3', cmap='RdBu_r')
    plt.show()

correlation_checker('smv')

df_pairplot = df[['smv', 'no_of_workers', 'over_time' ]]

sns.pairplot(df_pairplot)
plt.show()

# sum all wip per work_week
df_wip = df.groupby('work_week').sum()
df_wip['wip']

plt.figure(figsize=(8,5))
ax1 = plt.subplot()
sns.barplot(x= df_wip.index ,
            y= df_wip['wip'],
            palette = 'deep')

plt.title('Work in Progress')
for i in ax1.containers:
    ax1.bar_label(i,)

plt.tight_layout()    
plt.show()

# set variables
ww = df.groupby('work_week').sum() # we use groupby then sum
ww['cummulative'] = ww.over_time.cumsum()
ww[['no_of_workers', 'over_time', 'cummulative']]

#
plt.figure(figsize=(14, 6)) 

# set variables
x = ww.index
y = ww.over_time

# plot1
# sum of overtime per work_week
ax1 = plt.subplot(1,2,1)
sns.barplot(x = x, y = y, data=df)
ax1.set_xticklabels(ax1.get_xticklabels(),rotation = 80)
plt.title('Overtime per Work Week ')

for i in ax1.containers:
    ax1.bar_label(i,)


#plot2
# cumulative sum of work_week
ax2 = plt.subplot(1,2,2)
sns.lineplot(x=x, y=ww.over_time.cumsum(), marker='o', )
plt.title('Cumulative Overtime')

plt.tight_layout()
plt.show()

from scipy.stats import pearsonr
corr, p = pearsonr(df.incentive, df.actual_productivity)
corr, p

plt.figure(figsize=(10, 5)) 

ax1 = plt.subplot(1,2,1)
ax1 =  sns.regplot(x=df.actual_productivity, y=df.incentive, scatter_kws={'alpha':0.3})

# Zoom in 
ax2 = plt.subplot(1,2,2)
plt.title('Zoom in')
plt.scatter(df.actual_productivity, df.incentive, alpha=0.3)
plt.ylim(0,200)
plt.xlabel('actual_productivity')
plt.ylabel('incentive')

plt.show()

import warnings
warnings.filterwarnings('ignore')


# replace zero value in incentive w/ np.nan
df.incentive = df['incentive'].replace(0,np.nan)

# drop np.nan values
df_x = df.dropna(subset=['incentive'], axis=0)

# impute outliers with mean value
df_x['incentive'] = df_x['incentive'].where(df_x['incentive'] < 500, 
                                            df_x.incentive.mean())

corr, p = pearsonr(df_x.incentive, df_x.actual_productivity)
corr

plt.figure(figsize=(10, 5)) 

ax1 = plt.subplot(1,2,1)
ax1 =  sns.regplot(x=df_x.actual_productivity, y=df_x.incentive, scatter_kws={'alpha':0.3})

# Zoom in 
ax2 = plt.subplot(1,2,2)
plt.title('Zoom in')
plt.scatter(df_x.actual_productivity, df_x.incentive, alpha=0.3)
# plt.xlim(0,200)
plt.xlabel('actual_productivity')
plt.ylabel('incentive')

plt.show()

# check details for this extreme incentives
df[df['incentive']>2500]

ww10 = df[
          (df.work_week == 'ww_10') & 
          (df.day == 'Monday') &
          (df.department == 'finishing')
         ] 

ww10.sort_values(by='team')

len(df[df['incentive']>500])

# import sklean library
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score
import warnings

correlation_checker('no_of_workers')

# subset of our dataset
df_cluster = df[['no_of_workers', 'smv']]

# Normalise our data set
scaler = StandardScaler()
df_cluster = scaler.fit_transform(df_cluster)

# hide/show code
clusters = []
inertias = []

# Checking what is the best 'k' value
# try for 9 
for k in range(1,10):
  km = KMeans(n_clusters=k, random_state=0)
  km.fit(df_cluster)
  
  # append clusters and inertias
  inertias.append(km.inertia_)
  clusters.append(km)
  
plt.plot(range(len(inertias)), inertias, '-o')
plt.xlabel('number of clusters (k)')
plt.ylabel('inertia')
plt.show()

for i in range(1,9,1): 
    print("---------------------------------------")
    print(clusters[i])
    print("Silhouette score:",silhouette_score(df_cluster, clusters[i].predict(df_cluster)))

km = KMeans(n_clusters=2, random_state=0)

# fit and predict labels
y_km = km.fit_predict(df_cluster)
y_km

# Note:
# first column is the no_of_workers
# 2nd column is the smv
df_cluster[0:5]

no_of_workers = df_cluster[:,0]
smv = df_cluster[:,1]

# first 5 records
print(no_of_workers[0:5])
print(smv[0:5])

centers = km.cluster_centers_
centers

plt.rcParams['axes.facecolor'] = 'gray'

plt.figure(figsize=(8,5))
plt.scatter(no_of_workers, smv, c=y_km, alpha=0.5)

plt.scatter(centers[0][0], centers[0][1], marker='*', color ='purple', s=300, label ='purple centroid (0)', edgecolor = 'black')
plt.scatter(centers[1][0], centers[1][1], marker='*', color ='yellow', s=300, label ='blue centroid (1)', edgecolor = 'black')

plt.xlabel('number of workers')
plt.ylabel('smv(allocated time for a task)')
plt.legend()
plt.title('Scaled clustering')
plt.show()

# turn our predicted value first into dataframe
group_worker = pd.DataFrame(y_km)

# concatenate, assign it to a new variable
df_new = pd.concat([df, group_worker], axis=1, join='inner')

# # rename column
df_new = df_new.rename(columns={0: 'group_worker'})

plt.rcParams['axes.facecolor'] = 'gray'

sns.scatterplot(x= 'no_of_workers', y= 'smv', data=df_new, hue='group_worker', palette=['purple','yellow'])
plt.show()
plt.clf

df.no_of_style_change.value_counts()

import plotly.express as px
fig = px.pie(df, 'no_of_style_change')

fig.update_layout(title="Proportion of no_of_style_change")
fig.update_traces(textposition='inside',
                  textinfo='percent+label', showlegend=True)

fig.update_traces(pull= 0.02)
fig.show('png')

from scipy.stats import pearsonr
corr, p = pearsonr(df.no_of_style_change, df.idle_time)
corr, p

# visualize
plt.rcParams['axes.facecolor'] = 'white'
# sns.set(rc={'figure.figsize':(8,8)})

ax1 = sns.lmplot(x='no_of_style_change', y='idle_time', 
               data=df, 
               x_jitter=0.15, 
               y_jitter=0.15, 
               fit_reg=False, 
               scatter_kws={'alpha':0.5})

ax1.set(xticks=range(3))
ax1.set_xticklabels(df.no_of_style_change.unique())
ax1.set(ylim=(0,10))

# 2nd plot
ax2 = sns.lmplot(x='no_of_style_change', y='idle_time', 
               data=df, 
               x_jitter=0.15, 
               y_jitter=0.15, 
               fit_reg=False, 
               scatter_kws={'alpha':0.5})

ax2.set(xticks=range(3))
ax2.set_xticklabels(df.no_of_style_change.unique())
ax2.set(ylim=(0,.25))
plt.title('ZOOM IN')

plt.show()

change_style_idle = df[df['idle_time']>=0]

# we use count because we need to determine how many idle_time does the production had that is >= 2
change_style_idle.groupby('no_of_style_change').count()[['idle_time']]

jan = df[df['month']==1]
feb = df[df['month']==2]
mar = df[df['month']==3]

plt.hist(jan.actual_productivity, alpha=.5, label='January')
plt.hist(feb.actual_productivity, alpha=.5, label='February')
plt.hist(mar.actual_productivity, alpha=.5, label='March')

plt.xlabel('Actual Productivity')
plt.ylabel('Count')
plt.legend()
plt.show()

feb_higher_prod = feb[feb.actual_productivity >1][['actual_productivity']]
print(len(feb_higher_prod))
feb_higher_prod.sort_values(by='actual_productivity', ascending=False).head(3)

df['month_name'] = df.month.map({
    1:'January',
    2:'February',
    3:'March'
})

from scipy.stats import f_oneway
from statsmodels.stats.multicomp import pairwise_tukeyhsd

fstat, pval = f_oneway(jan.actual_productivity,
                       feb.actual_productivity,
                       mar.actual_productivity)
pval

tukey_results =  pairwise_tukeyhsd(df.actual_productivity, df.month_name, 0.05)
print(tukey_results)

sns.boxplot(x='month_name', y='actual_productivity', data=df )
plt.show()

#define functions
def showoutliers(df, column_name = ""):
        iqr = df[column_name].quantile(.75) - df[column_name].quantile(.25)
        
        # lower whisker
        lowerbound = (df[column_name].quantile(.25)) - iqr * 1.5 
        # upper whisker
        upperbound = (df[column_name].quantile(.75)) + iqr * 1.5
        
        # dfpoints beyond lower whisker
        lowerbound_outliers = df[df[column_name] < lowerbound]
        
        # adtapoint beyond upper whisker
        higherbound_outliers = df[df[column_name] > upperbound]
        
        # outliers
        outliers = pd.concat([lowerbound_outliers,higherbound_outliers])
        return outliers
def countoutliers(df, column_name = ""):
        iqr = df[column_name].quantile(.75) - df[column_name].quantile(.25)
        lowerbound = (df[column_name].quantile(.25)) - iqr * 1.5
        upperbound = (df[column_name].quantile(.75)) + iqr * 1.5
        lowerbound_outliers = df[df[column_name] < lowerbound]
        higherbound_outliers = df[df[column_name] > upperbound]
        outliers = pd.concat([lowerbound_outliers,higherbound_outliers])
        count = len(outliers)
        return {column_name : count}
def Replace_Outliers(df_name, value, column_name = ""):
    iqr = df_name[column_name].quantile(.75) - df_name[column_name].quantile(.25)
    
    lowerbound = (df_name[column_name].quantile(.25)) - iqr * 1.5
    upperbound = (df_name[column_name].quantile(.75)) + iqr * 1.5
            
    df_name[column_name] = np.where(df_name[column_name] > upperbound, value, df_name[column_name])
    df_name[column_name] = np.where(df_name[column_name] < lowerbound, value, df_name[column_name])

# selected features only
df_n = df[[ 'targeted_productivity', 'smv', 'wip', 'over_time', 'incentive', 'no_of_workers', 'actual_productivity']]
# df_n

column_list = df_n.columns
column_list = np.array(column_list)
for i in df_n:
    print (countoutliers(df_n, i))

for i in column_list:
    col = i
    perc = countoutliers(df_n, i)[i] / len(df_n)
    print (col + ': ' + str('{:.2f}'.format(perc*100)) + '%') 

df_n.plot(kind='box', 
          subplots=True, 
          sharey=False, 
          figsize=(20, 7))
# increase spacing between subplots
plt.subplots_adjust(wspace=0.5) 

plt.tight_layout()
plt.show()

df_new.incentive.describe()

Replace_Outliers(df_new, 50, 'incentive')

# visualize
plt.figure(figsize=(4,3))
sns.boxplot(x='incentive', data=df_new)
plt.show()