from numpy import array, zeros, argmin, inf, ndim
import scipy.stats
import time
import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import euclidean_distances
def dtw(x, y, warp=1):
assert len(x)
assert len(y)
r, c = len(x), len(y)
D0 = zeros((r + 1, c + 1))
D0[0, 1:] = inf
D0[1:, 0] = inf
D1 = D0[1:, 1:] # view
for i in range(r):
for j in range(c):
D1[i, j] = euclidean_distances([x[i]], [y[j]])
C = D1.copy()
for i in range(r):
for j in range(c):
min_list = [D0[i, j]]
for k in range(1, warp + 1):
i_k = min(i + k, r - 1)
j_k = min(j + k, c - 1)
min_list += [D0[i_k, j], D0[i, j_k]]
D1[i, j] += min(min_list)
return D1[-1, -1] / sum(D1.shape)
# 加载数据
def loadDataSet(fileName):
data = np.loadtxt(fileName, delimiter=",", skiprows=1)
return data
def saveDataSet(fileNmae, dataSet):
np.savetxt(fileNmae, dataSet, delimiter=',')
# 为给定数据集构建一个包含K个随机质心的集合
def randCent(dataSet, k):
m, n, p = dataSet.shape
centroids = np.zeros((k, n, p))
for i in range(k):
index = int(np.random.uniform(0, m)) #
centroids[i, :, :] = dataSet[index, :, :]
return centroids
# k均值聚类
def KMeans(dataSet, k):
m = np.shape(dataSet)[0] # 行的数目
# 第一列存样本属于哪一簇
# 第二列存样本的到簇的中心点的误差
clusterAssment = np.mat(np.zeros((m, 2)))
clusterChange = True
# 第1步 初始化centroids
centroids = randCent(dataSet, k)
while clusterChange:
clusterChange = False
# 遍历所有的样本(行数)
for i in range(m):
minDist = 100000.0
minIndex = -1
# 遍历所有的质心
# 第2步 找出最近的质心
for j in range(k):
# 计算该样本到质心的欧式距离
distance = dtw(centroids[j, :, :], dataSet[i, :, :])
if distance < minDist:
minDist = distance
minIndex = j
# 第 3 步:更新每一行样本所属的簇
if clusterAssment[i, 0] != minIndex:
clusterChange = True
clusterAssment[i, :] = minIndex, minDist ** 2
# 第 4 步:更新质心
for j in range(k):
pointsInCluster = dataSet[np.nonzero(clusterAssment[:, 0].A == j)[0]] # 获取簇类所有的点
centroids[j, :, :] = np.mean(pointsInCluster, axis=0) # 对矩阵的行求均值
print("Congratulations,cluster complete!")
return centroids, clusterAssment
def showCluster(dataSet, k, centroids, clusterAssment):
m, n, p = dataSet.shape
# 绘制所有的样本
for i in range(m):
markIndex = int(clusterAssment[i, 0])
plt.plot(dataSet[i, 0, 0], dataSet[i, 0, 1])
plt.show()
if __name__ == '__main__':
start = time.clock()
dataSet = loadDataSet("D:/151617-dtm-new/mds.csv")
m, n = dataSet.shape
dataSet_1 = np.zeros((m, n / 2, 2))
for i in range(m):
for j in range(n / 2):
dataSet_1[i, j, 0] = dataSet[i, 2 * j]
dataSet_1[i, j, 1] = dataSet[i, 2 * j + 1]
m, n, p = dataSet_1.shape
k = 4
centroids, clusterAssment = KMeans(dataSet_1, k)
saveDataSet('D:/151617-dtm-new/cluster.csv', clusterAssment)
# #
showCluster(dataSet_1, k, centroids, clusterAssment)
elapsed = (time.clock() - start)
print("Time used:", elapsed)