data.head(2)

data.tail(2)

data.shape

data.info()

data.describe()

data_example = data.loc['2012-01-01 00:15:00':'2012-01-03 00:00:00'][['MT_001','MT_002']]
data_example.plot()
plt.show()

# Check number of days
print(data2011.shape[0]/96)
print(data2012.shape[0]/96)
print(data2013.shape[0]/96)
print(data2014.shape[0]/96)

# See number of clients with 0 demand per year
print(sum(data2011.mean()==0))
print(sum(data2012.mean()==0))
print(sum(data2013.mean()==0))
print(sum(data2014.mean()==0))

clients = data2011.columns
clients_no_demand = clients[data2013.mean()==0] # clients with 0 demand
data_13_14 = data2013.append(data2014) # appending 2013 and 2014
data_13_14 = data_13_14.drop(clients_no_demand, axis=1) # drop clients with 0 demand
print(data_13_14.shape)
print(sum(data_13_14.mean()==0)) # check that there are no clients with 0 demand

data.head(3)

datagrouped = data.groupby("hour")
average_curves = datagrouped.agg("mean")
average_curves.shape 

average_curves_norm[['MT_001','MT_002','MT_369','MT_370']].plot()
plt.show()

inertia_list = []
inertia_list2 = []
sil_scores_list = []
sse = []

for n_clusters in test_clusters:
    kmeans = KMeans(n_clusters=n_clusters, init='k-means++', max_iter=300, random_state=42)
    kmeans_model = kmeans.fit(X)
    
    labels = kmeans_model.labels_    
    msl = metrics.silhouette_score(X, labels, metric='euclidean')
    sil_scores_list.append(msl)

    inertia_list.append(kmeans_model.inertia_)  #assign inertia of cluster to list
    inertia_list2.append(getInertia(X, kmeans_model))
    
    print("The silhoutte score for {0} clusters is {1}".format(n_clusters, msl))
    #try manual calc of squared mean error
    y_kmeans = kmeans_model.predict(X)
    centroids = kmeans_model.cluster_centers_
    curr_sse = 0
    for i in range(len(X)):
        curr_center = centroids[y_kmeans[i]]
        curr_sse += (X[i, 0] - curr_center[0]) ** 2 + (X[i, 1] - curr_center[1]) **2
    sse.append(curr_sse)

plt.figure(figsize=(10, 10))
plt.plot(range(1, 8), inertia_list2)
plt.title('Elbow Method - 1')
plt.xlabel('Number of clusters')
plt.ylabel('Squared Errors')
plt.show()

plt.figure(figsize=(10, 10))
plt.plot(range(1, 8), sse)
plt.title('Elbow Method - 2')
plt.xlabel('Number of clusters')
plt.ylabel('Squared cluster errors')
plt.show()

u_labels = np.unique(label)
u_labels

print("This is the shape of X:", len(X))
print("This is the shape of the averages from original data, but normalised", average_curves_norm.shape)
print("This is the shape of the averages from original data, but not normalised", average_curves.shape)
print("The shape of the labels_ is:", model.labels_.shape)
print("The shape of the predict is:", label.shape)
print("The shape of the y_predict is:", y_predict.shape)

from scipy import stats

mapping = {}
for class_id in np.unique(label):
    mode, _ = stats.mode(label[label==class_id])
    mapping[mode[0]] = class_id

mapping

plt.figure(figsize=(15, 15))
plt.plot(X[y_pred==0, 2], X[y_pred==0, 3], "yo", label="Cluster 1")
plt.plot(X[y_pred==1, 2], X[y_pred==1, 3], "bs", label="Cluster 2")
plt.plot(X[y_pred==2, 2], X[y_pred==2, 3], "g^", label="Cluster 3")
plt.xlabel("", fontsize=14)
plt.ylabel("", fontsize=14)
plt.legend(loc="upper left", fontsize=12)
plt.show()

np.unique(y_predict, return_counts=True)

cent2

group1.shape
group1

dat_plot1 = group1.drop('cluster', axis = 1)
dat_plot1.T.plot(figsize=(15, 15),legend = False)
cent1.plot(linewidth=5, color='red')
plt.title('There are 113 customers in the first cluster')
plt.show()

dat_plot2 = group2.drop('cluster', axis = 1)
dat_plot2.T.plot(figsize=(15, 15),legend = False)
cent2.plot(linewidth=5, color='red')
plt.title('There are 201 customers in the second cluster')
plt.show()

dat_plot3 = group3.drop('cluster', axis = 1)
dat_plot3.T.plot(figsize=(15, 15),legend = False)
cent3.plot(linewidth=5, color='red')
plt.title('There are 35 customers in the second cluster')
plt.show()

client = 'MT_022'
oneClient = data_13_14[client]
X = [] # a list of arrays, each array being a normalized curve for a day
for J in range(2*365):
    X.extend([np.array(oneClient[J*96:(J+1)*96])])#/np.mean(oneClient[J*96:(J+1)*96])])
len(X)

print("The length of one customer", len(oneClient))
print("The shape of customer file before conversion", oneClient.shape)
print("The length of X after conversion", len(X))

inertia_list = [] sil_scores_list = [] cent_list = [] for n_clusters in range_n_clusters: kmeans = KMeans(n_clusters=n_clusters, init='k-means++', max_iter=300, random_state=42) kmeans_model = kmeans.fit(X) labels = kmeans_model.labels_ msl = metrics.silhouette_score(X, labels, metric='euclidean') sil_scores_list.append(msl) inertia_list.append(kmeans_model.inertia_) #assign inertia of cluster to list # inertia_list2.append(getInertia(X, kmeans_model)) print("The silhoutte score for{0} clusters is{1}".format(n_clusters, msl)) y_kmeans = kmeans_model.predict(X) centroids = kmeans_model.cluster_centers_ cent_list.append(centroids) ############### #Plot the silouette analysis graphs X = np.array(X) for n_clusters in range_n_clusters: # Create a subplot with 1 row and 2 columns fig, (ax1, ax2) = plt.subplots(1, 2) fig.set_size_inches(18, 7) # The 1st subplot is the silhouette plot # The silhouette coefficient can range from -1, 1 but in this example all # lie within [-0.1, 1] ax1.set_xlim([-0.1, 1]) # The (n_clusters+1)*10 is for inserting blank space between silhouette # plots of individual clusters, to demarcate them clearly. ax1.set_ylim([0, len(X) + (n_clusters + 1) * 10]) # Initialize the clusterer with n_clusters value and a random generator # seed of 10 for reproducibility. clusterer = KMeans(n_clusters=n_clusters, random_state=10) cluster_labels = clusterer.fit_predict(X) # The silhouette_score gives the average value for all the samples. # This gives a perspective into the density and separation of the formed # clusters silhouette_avg = silhouette_score(X, cluster_labels) print("For n_clusters =", n_clusters, "The average silhouette_score is :", silhouette_avg) # Compute the silhouette scores for each sample sample_silhouette_values = silhouette_samples(X, cluster_labels) y_lower = 10 for i in range(n_clusters): # Aggregate the silhouette scores for samples belonging to # cluster i, and sort them ith_cluster_silhouette_values = \ sample_silhouette_values[cluster_labels == i] ith_cluster_silhouette_values.sort() size_cluster_i = ith_cluster_silhouette_values.shape[0] y_upper = y_lower + size_cluster_i color = cm.nipy_spectral(float(i) / n_clusters) ax1.fill_betweenx(np.arange(y_lower, y_upper), 0, ith_cluster_silhouette_values, facecolor=color, edgecolor=color, alpha=0.7) # Label the silhouette plots with their cluster numbers at the middle ax1.text(-0.05, y_lower + 0.5 * size_cluster_i, str(i)) # Compute the new y_lower for next plot y_lower = y_upper + 10 # 10 for the 0 samples ax1.set_title("The silhouette plot for the various clusters.") ax1.set_xlabel("The silhouette coefficient values") ax1.set_ylabel("Cluster label") # The vertical line for average silhouette score of all the values ax1.axvline(x=silhouette_avg, color="red", linestyle="--") ax1.set_yticks([]) # Clear the yaxis labels / ticks ax1.set_xticks([-0.1, 0, 0.2, 0.4, 0.6, 0.8, 1]) # 2nd Plot showing the actual clusters formed colors = cm.nipy_spectral(cluster_labels.astype(float) / n_clusters) ax2.scatter(X[:, 0], X[:, 1], marker='.', s=30, lw=0, alpha=0.7, c=colors, edgecolor='k') # Labeling the clusters centers = clusterer.cluster_centers_ # Draw white circles at cluster centers ax2.scatter(centers[:, 0], centers[:, 1], marker='o', c="white", alpha=1, s=200, edgecolor='k') for i, c in enumerate(centers): ax2.scatter(c[0], c[1], marker='$%d$' % i, alpha=1, s=50, edgecolor='k') ax2.set_title("The visualization of the clustered data.") ax2.set_xlabel("Feature space for the 1st feature") ax2.set_ylabel("Feature space for the 2nd feature") plt.suptitle(("Silhouette analysis for KMeans clustering on sample data " "with n_clusters =%d" % n_clusters), fontsize=14, fontweight='bold') plt.show()

plt.figure(figsize=(5, 5)) plt.plot(range(1, 10), inertia_list) plt.title('Inertia Values') plt.xlabel('Number of clusters') plt.ylabel('Inertia') plt.show()

plt.figure(figsize=(5, 5)) plt.plot(range(1, 10), sil_scores_list,) plt.title('Silhoutte Scores') plt.xlabel('Number of clusters') plt.ylabel('Silhoutte') plt.show()

len(daysyear)

len(X)

len(y_kmeans)

date_an = pd.DataFrame(daysyear, columns=['date']) date_an.head(5)

customer1 = pd.DataFrame(y_kmeans, columns=['cluster']) customer1.head(5)

#Analyze daily / weekly cluster 1 date_analysis = pd.concat([date_an, customer1], axis=1) date1 = date_analysis.loc[date_analysis['cluster']== 0] #date3 = date_analysis.loc[date_analysis['cluster']== 2] mon1 = date1.date.str.count("mon").sum() tue1 = date1.date.str.count("tue").sum() wed1 = date1.date.str.count("wed").sum() thu1 = date1.date.str.count("thu").sum() fri1 = date1.date.str.count("fri").sum() sat1 = date1.date.str.count("sat").sum() sun1 = date1.date.str.count("sun").sum() date1_analysis ={} for var in ["mon1", "tue1", "wed1", "thu1", "fri1", "sat1", "sun1"]: date1_analysis[var] = eval(var) # Analyze daily / weekly cluster 2 date2 = date_analysis.loc[date_analysis['cluster']== 1] mon2 = date2.date.str.count("mon").sum() tue2 = date2.date.str.count("tue").sum() wed2 = date2.date.str.count("wed").sum() thu2 = date2.date.str.count("thu").sum() fri2 = date2.date.str.count("fri").sum() sat2 = date2.date.str.count("sat").sum() sun2 = date2.date.str.count("sun").sum() date2_analysis ={} for var in ["mon2", "tue2", "wed2", "thu2", "fri2", "sat2", "sun2"]: date2_analysis[var] = eval(var) #Cluster 1 monthly analysis month1 = date1.date.str.count("-01-").sum() month2 = date1.date.str.count("-02-").sum() month3 = date1.date.str.count("-03-").sum() month4 = date1.date.str.count("-04-").sum() month5 = date1.date.str.count("-05-").sum() month6 = date1.date.str.count("-06-").sum() month7 = date1.date.str.count("-07-").sum() month8 = date1.date.str.count("-08-").sum() month9 = date1.date.str.count("-09-").sum() month10 = date1.date.str.count("-10-").sum() month11 = date1.date.str.count("-11-").sum() month12 = date1.date.str.count("-12-").sum() monthly_analysis1 ={} for var in ["month1", "month2", "month3", "month4", "month5", "month6", "month7", "month8", "month9", "month10", "month11", "month12"]: monthly_analysis1[var] = eval(var) #Cluster 2 monthly analysis month1_2 = date2.date.str.count("-01-").sum() month2_2 = date2.date.str.count("-02-").sum() month3_2 = date2.date.str.count("-03-").sum() month4_2 = date2.date.str.count("-04-").sum() month5_2 = date2.date.str.count("-05-").sum() month6_2 = date2.date.str.count("-06-").sum() month7_2 = date2.date.str.count("-07-").sum() month8_2 = date2.date.str.count("-08-").sum() month9_2 = date2.date.str.count("-09-").sum() month10_2 = date2.date.str.count("-10-").sum() month11_2 = date2.date.str.count("-11-").sum() month12_2 = date2.date.str.count("-12-").sum() monthly_analysis2 ={} for var in ["month1_2", "month2_2", "month3_2", "month4_2", "month5_2", "month6_2", "month7_2", "month8_2", "month9_2", "month10_2", "month11_2", "month12_2"]: monthly_analysis2[var] = eval(var) #Extract the analysis print("These are days of the week counts from the first cluster" , date1_analysis) print("These are days of the week counts from the second cluster", date2_analysis) print("The length of cluster 1:", len(date1)) print("The length of cluster 2:", len(date2)) print("The percetage of data points that fall on the weekend and holidays for cluster 1:", ((sat1 + sun1 + 7) / len(date1) *100)) #there were 7 wkday holidays in 2013 print("The percetage of data points that fall during the week for cluster 2:", ((mon2 + tue2 + wed2 + thu2 + fri2) / len(date2) *100)) print("These are the data points for cluster 1 monthly counts", monthly_analysis1) print("These are the data points for cluster 2 monthly counts", monthly_analysis2)

monthly_analysis1

monthly_analysis2