I am currently writing a K-means clustering algorithm in Python, and I seem to have coded myself into a corner... I begin my algorithm with data sets assigned randomly to the appropriate number of K clusters, but I always end up converging to a single cluster (convergence time depends on value of K). The only case where this doesn't happen is when K=N, my number of data sets. In that case, I keep 31 clusters (as expected - you can't be closer to anyone than yourself). The code requires specific file formatting, but I can supply some examples upon request. I've put the code at the bottom of the post (it's rather bulky), but my logic is as follows:
- Randomly assign N models (of size M) to K clusters
- Average the models in each K
- Calculate similarity (in place of euclidean distance) between N models and each K cluster
- Assign N models to "closest" similarity K.
- Repeat
For benchmarking purposes, 31 data sets initialized to 4 clusters will converge to a single cluster in 3-5 iterations.
import numpy as np
from math import pow,sqrt
################### Function for scoring the differences between clusters ############################
def Score(PROFILES):
NUM=len(PROFILES)
scores=np.zeros((NUM,NUM))
for x in range(NUM):
for y in range(x+1,NUM):
scores[x,y]=Similarity(PROFILES[x],PROFILES[y])
scores[y,x]=scores[x,y]
scores[y,y]=0
return scores
######################################################################################################
################################### Function for calculating the score 'distance' between two clusters #
def Similarity(PROFILE1,PROFILE2):
profile1=np.log(PROFILE1)
profile2=np.log(PROFILE2)
sum=0
for idx in range(len(profile1)):
sum+=pow(profile1[idx]-profile2[idx],2)/(2*(0.1**2)) # Denominator is approximate error, numerator is difference in I(q)
score=sqrt(sum)
return score
######################################################################################################
########### Determine Initial Clusters ###############################################################
def FirstGrouping(PROFILES,K):
Nprofiles=len(PROFILES)
CLUST_IDX=np.random.randint(1,K+1,size=Nprofiles) #Assign each profile a cluster number
AVERAGES=AverageCluster(PROFILES,CLUST_IDX) #Calculate cluster Average Profile
return CLUST_IDX,AVERAGES #List of cluster numbers and KxN array of cluster average profiles
######################################################################################################
########## Calculate average of profiles in each cluster #############################################
def AverageCluster(PROFILES,CLUST_IDX):
CLUSTER_AVERAGE=np.zeros((np.amax(CLUST_IDX),len(PROFILES[0]))) # Initialize a KxN array to store the I(q) (for N q-values) for the K clusters
for CLUSTER in CLUST_IDX:
PROFILES_IN_CLUST=np.where(CLUST_IDX==CLUSTER)[0] # Find the profiles belonging to CLUSTER
NUMBER_IN_CLUST=len(PROFILES_IN_CLUST) # Count the total number in the cluster
for PROFILE in PROFILES_IN_CLUST:
CLUSTER_AVERAGE[CLUSTER-1]+=PROFILES[PROFILE] #Calculate the sum of the scattering profiles
CLUSTER_AVERAGE[CLUSTER-1]=CLUSTER_AVERAGE[CLUSTER-1]/NUMBER_IN_CLUST #Divide by number in cluster to get the average
return CLUSTER_AVERAGE # KxN array
######################################################################################################
############## Find closest Center (Cluster average) for each profile ################################
def FindClosest(PROFILES,CLUST_IDX,CLUSTERS): #MxN array of individual profiles, K-array of cluster identities, KxN array of cluster average profiles
NEW_CLUSTERS=np.copy(CLUST_IDX) #Keep the old cluster index, but make a copy for editing to the new one
for profile in range(len(PROFILES)):
Similarity_to_clusters=np.zeros(len(CLUSTERS)) #Initialize a length K array to store similarity values within
for cluster in range(len(CLUSTERS)):
Similarity_to_clusters[cluster]=Similarity(PROFILES[profile],CLUSTERS[cluster])
print("Profile %i, Cluster %i - %f" % (profile,cluster,Similarity_to_clusters[cluster]))
NEW_CLUSTERS[profile]=np.where(Similarity_to_clusters==np.amin(Similarity_to_clusters))[0]+1 #Assign the profile to the closest cluster in similarity
NEW_AVERAGE=AverageCluster(PROFILES,NEW_CLUSTERS)
return NEW_CLUSTERS,NEW_AVERAGE,CLUST_IDX,CLUSTERS #new cluster dictionary, new cluster average, old cluster dictionary, old cluster average
######################################################################################################
################# Function for stopping the while loop ###############################################
def StopCondition(OLD_CLUSTER,NEW_CLUSTER,ITERATION,MAX):
if ITERATION > MAX: # Stop if you have surpassed max iteration count
print("Cluster search reached maximum number of iterations(%i)" % MAX)
return 1
print(np.array_equal(NEW_CLUSTER,OLD_CLUSTER))
if np.array_equal(NEW_CLUSTER,OLD_CLUSTER): # Stop if you have no change in cluster identities
print("Cluster search has converged! (after %i iterations)" % ITERATION)
return 1
######################################################################################################
########### Function for K-Means Clustering of Profiles ##############################################
def KMeansClustering(PROFILES,K):
#Bookkeeping variables
iteration=0
max_iteration=100
#Assign first group of clusters
old_idx=None
new_idx,new_avg=FirstGrouping(PROFILES,K)
while not StopCondition(old_idx,new_idx,iteration,max_iteration):
iteration+=1
new_idx,new_avg,old_idx,old_avg=FindClosest(PROFILES,new_idx,new_avg)
print("OLD:\n"+str(old_idx))
print("NEW:\n"+str(new_idx))
return new_idx,new_avg
######################################################################################################
################### Run the Clustering ###############################################################
num=31 # number of data sets
K=4 # number of desired clusters
for x in range(num):
if x==0:
profiles=np.genfromtxt("file."+str(x),usecols=2,skip_header=3)
else:
profiles=np.vstack((profiles,np.genfromtxt("file"+str(x),usecols=2,skip_header=3)
cluster_index,cluster_average=KMeansClustering(profiles,K)
print(cluster_index)
print("Clustering complete.")
######################################################################################################