Clustering using k-medoids

This is exactly the implementation found in NumPy and SciPy Recipes for Data Science on k-Medoids Clustering but with some indentation mistakes (probably due to the copy & paste process). To understand the indentation being used in the algorithm is might be worth to read the answers on the for...else syntax:

https://stackoverflow.com/questions/9979970/why-does-python-use-else-after-for-and-while-loops#9980752

A working script with the correct indentation and some improvement is found at https://github.com/letiantian/kmedoids/blob/master/kmedoids.py :

import numpy as np
import random

def kMedoids(D, k, tmax=100):
   # determine dimensions of distance matrix D
   m, n = D.shape

   if k > n:
       raise Exception('too many medoids')

   # find a set of valid initial cluster medoid indices since we
   # can't seed different clusters with two points at the same location
   valid_medoid_inds = set(range(n))
   invalid_medoid_inds = set([])
   rs,cs = np.where(D==0)
   # the rows, cols must be shuffled because we will keep the first duplicate below
   index_shuf = list(range(len(rs)))
   np.random.shuffle(index_shuf)
   rs = rs[index_shuf]
   cs = cs[index_shuf]
   for r,c in zip(rs,cs):
       # if there are two points with a distance of 0...
       # keep the first one for cluster init
       if r < c and r not in invalid_medoid_inds:
           invalid_medoid_inds.add(c)
   valid_medoid_inds = list(valid_medoid_inds - invalid_medoid_inds)

   if k > len(valid_medoid_inds):
       raise Exception('too many medoids (after removing {} duplicate points)'.format(
        len(invalid_medoid_inds)))

   # randomly initialize an array of k medoid indices
   M = np.array(valid_medoid_inds)
   np.random.shuffle(M)
   M = np.sort(M[:k])

   # create a copy of the array of medoid indices
   Mnew = np.copy(M)

   # initialize a dictionary to represent clusters
   C = {}
   for t in xrange(tmax):
       # determine clusters, i. e. arrays of data indices
       J = np.argmin(D[:,M], axis=1)
       for kappa in range(k):
           C[kappa] = np.where(J==kappa)[0]
       # update cluster medoids
       for kappa in range(k):
           J = np.mean(D[np.ix_(C[kappa],C[kappa])],axis=1)
           j = np.argmin(J)
           Mnew[kappa] = C[kappa][j]
       np.sort(Mnew)
       # check for convergence
       if np.array_equal(M, Mnew):
           break
       M = np.copy(Mnew)
   else:
       # final update of cluster memberships
       J = np.argmin(D[:,M], axis=1)
       for kappa in range(k):
           C[kappa] = np.where(J==kappa)[0]

   # return results
   return M, C

Note that the else as well as the return are at the same indentation as the for-loop.

However, I recommend to change the second inner loop

# update cluster medoids
for kappa in range(k):
       J = np.mean(D[np.ix_(C[kappa],C[kappa])],axis=1)
       j = np.argmin(J)
       Mnew[kappa] = C[kappa][j]

to

# update cluster medoids
for kappa in range(k):
       J = np.mean(D[np.ix_(C[kappa],C[kappa])],axis=1)
       # Fix for the low-idx bias by J.Nagele (2019):
       shuffled_idx = np.arange(len(J))
       np.random.shuffle(shuffled_idx) 
       j = shuffled_idx[np.argmin(J[shuffled_idx])] 
       Mnew[kappa] = C[kappa][j]

as otherwise the resulting medoids depend on the order of the input. That is because the argmin command returns always the lowest index in case multiple potential medoids have equal distances to the cluster members.