#-*- coding:utf-8 -*- 
#!/usr/bin/python 
 
''''' 
k Means K均值聚類 
''' 
# 測(cè)試 
# Ｋ均值聚類   import kMeans as KM KM.kMeansTest() 
# 二分Ｋ均值聚類 import kMeans as KM KM.biKMeansTest() 
# 地理位置 二分Ｋ均值聚類 import kMeans as KM KM.clusterClubs() 
from numpy import * 
 
# 導(dǎo)入數(shù)據(jù)集 
def loadDataSet(fileName):   #  
  dataMat = []        #  
  fr = open(fileName) 
  for line in fr.readlines(): # 每一行 
    curLine = line.strip().split('\t')# 按 Tab鍵 分割成 列表 
    fltLine = map(float,curLine)   # 映射成 浮點(diǎn)型 
    dataMat.append(fltLine)      # 放入數(shù)據(jù)集里 
  return dataMat 
 
# 計(jì)算歐幾里的距離 
def distEclud(vecA, vecB): 
  return sqrt(sum(power(vecA - vecB, 2))) #la.norm(vecA-vecB) 
 
# 初始構(gòu)建質(zhì)心(隨機(jī)) 數(shù)據(jù)集 質(zhì)心個(gè)數(shù) 
def randCent(dataSet, k): 
  n = shape(dataSet)[1] # 樣本特征維度 
  centroids = mat(zeros((k,n))) # 初始化 k個(gè) 質(zhì)心 
  for j in range(n):  # 每種樣本特征 
    minJ = min(dataSet[:,j]) # 每種樣本特征最小值 需要轉(zhuǎn)換成 numpy 的mat 
    rangeJ = float(max(dataSet[:,j]) - minJ)#每種樣本特征的幅值范圍 
    centroids[:,j] = mat(minJ + rangeJ * random.rand(k,1)) 
    # 在每種樣本的最大值和最小值間隨機(jī)生成K個(gè)樣本特征值 
  return centroids 
 
# 簡(jiǎn)單k均值聚類算法  
#    數(shù)據(jù)集 中心數(shù)量  距離算法      初始聚類中心算法  
def kMeans(dataSet, k, distMeas=distEclud, createCent=randCent): 
  m = shape(dataSet)[0]       # 樣本個(gè)數(shù) 
  clusterAssment = mat(zeros((m,2)))# 樣本標(biāo)記 分配結(jié)果 第一列索引 第二列誤差 
  centroids = createCent(dataSet, k)# 初始聚類中心 
  clusterChanged = True# 設(shè)置質(zhì)心是否仍然發(fā)送變化 
  while clusterChanged: 
    clusterChanged = False 
    for i in range(m): #對(duì)每個(gè)樣本 計(jì)算最近的中心 
    # 更新 樣本所屬關(guān)系 
      minDist = inf; minIndex = -1 # 距離變量 以及 最近的中心索引 
      for j in range(k): # 對(duì)每個(gè)中心 
        distJI = distMeas(centroids[j,:],dataSet[i,:])# 計(jì)算距離 
        if distJI < minDist: 
          minDist = distJI; minIndex = j# 得到最近的 中心 索引 
      if clusterAssment[i,0] != minIndex: clusterChanged = True  
      # 所屬索引發(fā)生了變化 即質(zhì)心還在變化，還可以優(yōu)化 
      clusterAssment[i,:] = minIndex,minDist**2 # 保存 所屬索引 以及距離平方 用以計(jì)算誤差平方和 SSE 
    # 更新質(zhì)心 
    print centroids # 每次迭代打印質(zhì)心 
    for cent in range(k):#  
      ptsInClust = dataSet[nonzero(clusterAssment[:,0].A==cent)[0]]# 數(shù)組過(guò)濾 得到各個(gè)中心所屬的樣本 
      centroids[cent,:] = mean(ptsInClust, axis=0) # 按列求平均 得到新的中心 
  return centroids, clusterAssment# 返回質(zhì)心 和各個(gè)樣本分配結(jié)果 
 
def kMeansTest(k=5): 
  MyDatMat = mat(loadDataSet("testSet.txt")) 
  MyCenters, ClustAssing = kMeans(MyDatMat, k) 
 
# bisecting K-means 二分K均值算法 克服局部最優(yōu)值 
def biKmeans(dataSet, k, distMeas=distEclud): 
  m = shape(dataSet)[0]       # 樣本個(gè)數(shù) 
  clusterAssment = mat(zeros((m,2)))# 樣本標(biāo)記 分配結(jié)果 第一列索引 第二列誤差 
  centroid0 = mean(dataSet, axis=0).tolist()[0]# 創(chuàng)建一個(gè)初始質(zhì)心 
  centList =[centroid0] # 一個(gè)中心的 列表 
  for j in range(m):  # 計(jì)算初始誤差 
    clusterAssment[j,1] = distMeas(mat(centroid0), dataSet[j,:])**2#每個(gè)樣本與中心的距離平方 
  while (len(centList) < k):# 中心數(shù)倆個(gè)未達(dá)到指定中心數(shù)量 繼續(xù)迭代 
    lowestSSE = inf    # 最小的 誤差平方和 SSE 
    for i in range(len(centList)):# 對(duì)于每一個(gè)中心 
      ptsInCurrCluster = dataSet[nonzero(clusterAssment[:,0].A==i)[0],:] # 處于當(dāng)前中心的樣本點(diǎn) 
      centroidMat, splitClustAss = kMeans(ptsInCurrCluster, 2, distMeas) # 對(duì)此中心內(nèi)的點(diǎn)進(jìn)行二分類 
      # 該樣本中心 二分類之后的 誤差平方和 SSE 
    sseSplit = sum(splitClustAss[:,1]) 
      # 其他未劃分?jǐn)?shù)據(jù)集的誤差平方和 SSE 
      sseNotSplit = sum(clusterAssment[nonzero(clusterAssment[:,0].A!=i)[0],1]) 
      print "sseSplit, and notSplit: ",sseSplit,sseNotSplit 
      # 劃分后的誤差和沒(méi)有進(jìn)行劃分的數(shù)據(jù)集的誤差為本次誤差 
      if (sseSplit + sseNotSplit) < lowestSSE: # 小于上次 的 誤差  
        bestCentToSplit = i # 記錄應(yīng)該被劃分的中心 的索引 
        bestNewCents = centroidMat # 最好的新劃分出來(lái)的中心 
        bestClustAss = splitClustAss.copy()# 新中心 對(duì)于的 劃分記錄 索引(0或1)以及 誤差平方 
        lowestSSE = sseSplit + sseNotSplit # 更新總的 誤差平方和 
    # 記錄中心 劃分 數(shù)據(jù) 
    bestClustAss[nonzero(bestClustAss[:,0].A == 1)[0],0] = len(centList) # 現(xiàn)有中心數(shù)量 
    bestClustAss[nonzero(bestClustAss[:,0].A == 0)[0],0] = bestCentToSplit# 最應(yīng)該被劃分的中心 
    print 'the bestCentToSplit is: ',bestCentToSplit 
    print 'the len of bestClustAss is: ', len(bestClustAss) 
    # 將最應(yīng)該被劃分的中心 替換為 劃分后的 兩個(gè) 中心(一個(gè)替換，另一個(gè) append在最后添加) 
    centList[bestCentToSplit] = bestNewCents[0,:].tolist()[0]# 替換 
    centList.append(bestNewCents[1,:].tolist()[0])      # 添加 
    # 更新 樣本標(biāo)記 分配結(jié)果 替換 被劃分中心的記錄 
    clusterAssment[nonzero(clusterAssment[:,0].A == bestCentToSplit)[0],:]= bestClustAss 
  return mat(centList), clusterAssment 
 
def biKMeansTest(k=5): 
  MyDatMat = mat(loadDataSet("testSet.txt")) 
  MyCenters, ClustAssing = biKmeans(MyDatMat, k) 
 
####位置數(shù)據(jù)聚類測(cè)試##### 
# 利用雅虎的服務(wù)器將地址轉(zhuǎn)換為 經(jīng)度和緯度 
import urllib 
import json 
def geoGrab(stAddress, city): 
  apiStem = 'http://where.yahooapis.com/geocode?' #  
  params = {} 
  params['flags'] = 'J'    # 設(shè)置返回類型為JSON字符串  
  params['appid'] = 'aaa0VN6k' # 注冊(cè) 帳號(hào)后獲得 http://developer.yahoo.com 
  params['location'] = '%s %s' % (stAddress, city) # 位置信息 
  url_params = urllib.urlencode(params)# 將字典轉(zhuǎn)換成可以通過(guò)ＵＲＬ進(jìn)行傳遞的字符串格式 
  yahooApi = apiStem + url_params   # 加入網(wǎng)絡(luò)地址  
  print yahooApi            # 打印 ＵＲＬ 
  c=urllib.urlopen(yahooApi)      # 打開(kāi) ＵＲＬ 
  return json.loads(c.read())     # 讀取返回的jason字符串  對(duì)位置進(jìn)行了編碼 得到經(jīng)度和緯度  
 
 
from time import sleep 
def massPlaceFind(fileName): 
  fw = open('places.txt', 'w') # 打開(kāi)位置信息文件 
  for line in open(fileName).readlines():# 每一行 
    line = line.strip() 
    lineArr = line.split('\t')# 得到列表 
    retDict = geoGrab(lineArr[1], lineArr[2])# 第二列為號(hào)牌 第三列為城市 進(jìn)行地址解碼 
    if retDict['ResultSet']['Error'] == 0: 
      lat = float(retDict['ResultSet']['Results'][0]['latitude']) #經(jīng)度 
      lng = float(retDict['ResultSet']['Results'][0]['longitude'])#緯度 
      print "%s\t%f\t%f" % (lineArr[0], lat, lng) 
      fw.write('%s\t%f\t%f\n' % (line, lat, lng)) #再寫(xiě)入到文件 
    else: print "error fetching" 
    sleep(1)#延遲1s 
  fw.close() 
 
# 返回地球表面兩點(diǎn)之間的距離 單位英里 輸入經(jīng)緯度(度) 球面余弦定理 
def distSLC(vecA, vecB):#Spherical Law of Cosines 
  a = sin(vecA[0,1]*pi/180) * sin(vecB[0,1]*pi/180) 
  b = cos(vecA[0,1]*pi/180) * cos(vecB[0,1]*pi/180) * \ 
           cos(pi * (vecB[0,0]-vecA[0,0]) /180) 
  return arccos(a + b)*6371.0 #pi in numpy 
 
 
# 位置聚類測(cè)試 畫(huà)圖可視化顯示 
import matplotlib 
import matplotlib.pyplot as plt 
 
def clusterClubs(numClust=5): 
  datList = [] # 樣本 
  for line in open('places.txt').readlines(): 
    lineArr = line.split('\t') 
    datList.append([float(lineArr[4]), float(lineArr[3])])# 保存經(jīng)緯度 
  datMat = mat(datList)# 數(shù)據(jù)集 numpy的mat類型 
  # 進(jìn)行二分K均值算法聚類 
  myCentroids, clustAssing = biKmeans(datMat, numClust, distMeas=distSLC) 
  fig = plt.figure()# 窗口 
  rect=[0.1,0.1,0.8,0.8] 
  scatterMarkers=['s', 'o', '^', '8', 'p', \ 
          'd', 'v', 'h', '>', '<'] 
  axprops = dict(xticks=[], yticks=[]) 
  ax0=fig.add_axes(rect, label='ax0', **axprops)#軸 
  imgP = plt.imread('Portland.png') # 標(biāo)注在實(shí)際的圖片上 
  ax0.imshow(imgP) 
  ax1=fig.add_axes(rect, label='ax1', frameon=False) 
  for i in range(numClust):#每一個(gè)中心 
    ptsInCurrCluster = datMat[nonzero(clustAssing[:,0].A==i)[0],:]# 屬于每個(gè)中心的樣本點(diǎn) 
    markerStyle = scatterMarkers[i % len(scatterMarkers)]# 點(diǎn)的類型 畫(huà)圖 
    # 散點(diǎn)圖 每個(gè)中心的樣本點(diǎn) 
    ax1.scatter(ptsInCurrCluster[:,0].flatten().A[0], ptsInCurrCluster[:,1].flatten().A[0], marker=markerStyle, s=90) 
  # 散 點(diǎn)圖 每個(gè)中心 
  ax1.scatter(myCentroids[:,0].flatten().A[0], myCentroids[:,1].flatten().A[0], marker='+', s=300) 
plt.show()# 顯示