cluster.py

#!/usr/bin/env python
# -*- coding: UTF-8 -*-

import sys
import math
import logging
import numpy as np

logger = logging.getLogger("dpc_cluster")

def load_paperdata(distance_f):
	'''
	Load distance from data

	Args:
		distance_f : distance file, the format is column1-index 1, column2-index 2, column3-distance
	
	Returns:
	    distances dict, max distance, min distance, max continues id
	'''
	logger.info("PROGRESS: load data")
	distances = {}
	min_dis, max_dis = sys.float_info.max, 0.0
	max_id = 0
	with open(distance_f, 'r') as fp:
		for line in fp:
			x1, x2, d = line.strip().split(' ')
			x1, x2 = int(x1), int(x2)
			max_id = max(max_id, x1, x2)
			dis = float(d)
			min_dis, max_dis = min(min_dis, dis), max(max_dis, dis)
			distances[(x1, x2)] = float(d)
			distances[(x2, x1)] = float(d)
	for i in xrange(max_id):
		distances[(i, i)] = 0.0
	logger.info("PROGRESS: load end")
	return distances, max_dis, min_dis, max_id

def select_dc(max_id, max_dis, min_dis, distances, auto = False):
	'''
	Select the local density threshold, default is the method used in paper, auto is `autoselect_dc`

	Args:
		max_id    : max continues id
		max_dis   : max distance for all points
		min_dis   : min distance for all points
		distances : distance dict
		auto      : use auto dc select or not
	
	Returns:
	    dc that local density threshold
	'''
	logger.info("PROGRESS: select dc")
	if auto:
		return autoselect_dc(max_id, max_dis, min_dis, distances)
	percent = 2.0
	position = int(max_id * (max_id + 1) / 2 * percent / 100)
	dc = sorted(distances.values())[position * 2 + max_id]
	logger.info("PROGRESS: dc - " + str(dc))
	return dc

def autoselect_dc(max_id, max_dis, min_dis, distances):
	'''
	Auto select the local density threshold that let average neighbor is 1-2 percent of all nodes.

	Args:
		max_id    : max continues id
		max_dis   : max distance for all points
		min_dis   : min distance for all points
		distances : distance dict
	
	Returns:
	    dc that local density threshold
	'''
	dc = (max_dis + min_dis) / 2

	while True:
		nneighs = sum([1 for v in distances.values() if v < dc]) / max_id ** 2
		if nneighs >= 0.01 and nneighs <= 0.002:
			break
		# binary search
		if nneighs < 0.01:
			min_dis = dc
		else:
			max_dis = dc
		dc = (max_dis + min_dis) / 2
		if max_dis - min_dis < 0.0001:
			break
	return dc

def local_density(max_id, distances, dc, guass=True, cutoff=False):
	'''
	Compute all points' local density

	Args:
		max_id    : max continues id
		distances : distance dict
		gauss     : use guass func or not(can't use together with cutoff)
		cutoff    : use cutoff func or not(can't use together with guass)
	
	Returns:
	    local density vector that index is the point index that start from 1
	'''
	assert guass and cutoff == False and guass or cutoff == True
	logger.info("PROGRESS: compute local density")
	guass_func = lambda dij, dc : math.exp(- (dij / dc) ** 2)
	cutoff_func = lambda dij, dc: 1 if dij < dc else 0
	func = guass and guass_func or cutoff_func
	rho = [-1] + [0] * max_id
	for i in xrange(1, max_id):
		for j in xrange(i + 1, max_id + 1):
			rho[i] += func(distances[(i, j)], dc)
			rho[j] += func(distances[(i, j)], dc)
		if i % (max_id / 10) == 0:
			logger.info("PROGRESS: at index #%i" % (i))
	return np.array(rho, np.float32)

def min_distance(max_id, max_dis, distances, rho):
	'''
	Compute all points' min distance to the higher local density point(which is the nearest neighbor)

	Args:
		max_id    : max continues id
		max_dis   : max distance for all points
		distances : distance dict
		rho       : local density vector that index is the point index that start from 1
	
	Returns:
	    min_distance vector, nearest neighbor vector
	'''
	logger.info("PROGRESS: compute min distance to nearest higher density neigh")
	sort_rho_idx = np.argsort(-rho)
	delta, nneigh = [0.0] + [float(max_dis)] * (len(rho) - 1), [0] * len(rho)
	delta[sort_rho_idx[0]] = -1.
	for i in xrange(1, max_id):
		for j in xrange(0, i):
			old_i, old_j = sort_rho_idx[i], sort_rho_idx[j]
			if distances[(old_i, old_j)] < delta[old_i]:
				delta[old_i] = distances[(old_i, old_j)]
				nneigh[old_i] = old_j
		if i % (max_id / 10) == 0:
			logger.info("PROGRESS: at index #%i" % (i))
	delta[sort_rho_idx[0]] = max(delta)
	return np.array(delta, np.float32), np.array(nneigh, np.float32)


class DensityPeakCluster(object):
	def local_density(self, load_func, distance_f, dc = None, auto_select_dc = False):
		'''
		Just compute local density

		Args:
		    load_func     : load func to load data
		    distance_f    : distance data file
		    dc            : local density threshold, call select_dc if dc is None
		    autoselect_dc : auto select dc or not

		Returns:
		    distances dict, max distance, min distance, max index, local density vector
		'''
		assert not (dc != None and auto_select_dc)
		distances, max_dis, min_dis, max_id = load_func(distance_f)
		if dc == None:
			dc = select_dc(max_id, max_dis, min_dis, distances, auto = auto_select_dc)
		rho = local_density(max_id, distances, dc)
		return distances, max_dis, min_dis, max_id, rho, dc

	def cluster(self, load_func, distance_f, density_threshold, distance_threshold, dc = None, auto_select_dc = False):
		'''
		Cluster the data

		Args:
		    load_func          : load func to load data
		    distance_f         : distance data file
		    dc                 : local density threshold, call select_dc if dc is None
		    density_threshold  : local density threshold for choosing cluster center
		    distance_threshold : min distance threshold for choosing cluster center
		    autoselect_dc      : auto select dc or not

		Returns:
		    local density vector, min_distance vector, nearest neighbor vector
		'''	
		assert not (dc != None and auto_select_dc)
		distances, max_dis, min_dis, max_id , rho , dc= self.local_density(load_func, distance_f, dc = dc, auto_select_dc = auto_select_dc)
		delta, nneigh = min_distance(max_id, max_dis, distances, rho)
		logger.info("PROGRESS: start cluster")
		cluster, ccenter = {}, {}  #cl/icl in cluster_dp.m

		for idx, (ldensity, mdistance, nneigh_item) in enumerate(zip(rho, delta, nneigh)):
			if idx == 0: continue
			if ldensity >= density_threshold and mdistance >= distance_threshold:
				ccenter[idx] = idx
				cluster[idx] = idx
			else:
				cluster[idx] = -1

		#assignation
		ordrho=np.argsort(-rho)
		for i in range(ordrho.shape[0]-1):
			if ordrho[i] == 0: continue
			if cluster[ordrho[i]] == -1:
				cluster[ordrho[i]] = cluster[nneigh[ordrho[i]]]
			if i % (max_id / 10) == 0:
				logger.info("PROGRESS: at index #%i" % (i))
		
		#halo
		halo, bord_rho = {},{}
		for i in range(1,ordrho.shape[0]):
			halo[i] = cluster[i]
		if len(ccenter) > 0:
			for idx in ccenter.keys():
				bord_rho[idx] = 0.0
			for i in range(1,rho.shape[0]-1):
				for j in range(i + 1,rho.shape[0]):
					if cluster[i] != cluster[j] and distances[i,j] <= dc:
						rho_aver = (rho[i] + rho[j]) / 2.0
						if rho_aver > bord_rho[cluster[i]]:
							bord_rho[cluster[i]] = rho_aver
						if rho_aver > bord_rho[cluster[j]]:
							bord_rho[cluster[j]] = rho_aver
			for i in range(1,rho.shape[0]):
				if rho[i] < bord_rho[cluster[i]]:
					halo[i] = 0
		for i in range(1,rho.shape[0]):
			if halo[i] == 0:
				cluster[i] =- 1

		self.cluster, self.ccenter = cluster, ccenter
		self.distances = distances
		self.max_id = max_id
		logger.info("PROGRESS: ended")
		return rho, delta, nneigh