#include<stdio.h>
#include <cstdio>
#include<string>
#include<math.h>
#include<stdlib.h>
#include<set>
#include<map>
#include<vector>
#include<queue>
#include<string.h>
#include<algorithm>
#include<iostream>
#include<time.h>
using namespace std;
//Author: I-Hsin
const int maxn=100000;//原始数据量的上界
vector<double>node[maxn];//存储原始数据
int K;//簇的数目
vector<double>center[maxn];//存储均值中心
vector< vector<double> >clusters[maxn];//存储k个聚类
int number;//输入数据的数目
int dimension;//输入数据的维数
double EuclideanDistance(vector<double> a,vector<double> b)//计算两个数据之间的欧几里得距离
{
double sum=0;
for(int i=0;i<a.size();i++)
{
sum+=(a[i]-b[i])*(a[i]-b[i]);
}
sum=sqrt(sum);
return sum;
}
void ChooseSeeds()//选择k个点作为初始的聚类中心，此处用k-means++算法优化，不是随机选择
{
srand((unsigned int) time(NULL));
int idx=rand()%number;
int cnt=0;
center[cnt++]=node[idx];//记录选择的均值中心
double dis[maxn];//记录每个点距离它最近的均值中心的距离
memset(dis,0x3f,sizeof(dis));
while(cnt<K)//求出每个点与距离它最近的均值中心的距离
{
double sum=0;
for(int i=0;i<number;i++)
{
for(int j=0;j<cnt;j++)
{
if(i==j) continue;
dis[i]=min(dis[i],EuclideanDistance(node[i],center[j]));
}
sum+=dis[i];
}
for(int i=0;i<number;i++)//归一化，其后可以对应到概率
{
dis[i]/=sum;
}
double cumprobs[maxn];
cumprobs[0]=dis[0];
for(int i=1;i<number;i++)//求出概率的前缀和
{
cumprobs[i]=dis[i]+cumprobs[i-1];
}
bool used[maxn];
memset(used,true,sizeof(used));
used[idx]=false;
next:
double r= (rand()%1000)*0.001;
bool flg=true;
for(int i=0;i<number;i++)
{
if(r<cumprobs[i]&&used[i])//选择满足概率的点作为簇中心
{
idx=i;
flg=false;
used[i]=false;
break;
}
}
if(flg) goto next; //如果没有找到，重新产生随机数r继续找
center[cnt++]=node[idx];
}
}
int ClusterOfElement(vector<double>v)//根据均值中心点集判断某一点属于哪个簇
{
double dist=0x3f3f3f3f;
double tmp;
int flag=rand()%K;//属于第flag个簇
for(int i=0;i<K;i++)//找出距离该点最近的均值中心
{
tmp=EuclideanDistance(center[i],v);
if(tmp<dist)
{
dist=tmp;
flag=i;
}
}
return flag;
}
vector<double> GetMeanValue(vector<vector<double> > sample)//更新某个簇的均值
{
vector<double>tmp(dimension,0);
int clunumber=sample.size();
for (int i=0;i<clunumber;i++)//遍历某个簇中的每一个元组（向量）
{
for(int j=0;j<sample[i].size();j++)//遍历某个向量的每个分量
{
tmp[j]+=sample[i][j];//将对应维的分量相加
}
}
for(int j=0;j<dimension;j++)
tmp[j]/=clunumber;//将对应维的分量相加后取均值
return tmp;
}
double GetError()//cluster是三维的，每个cluster[i]后面跟着许多点，每一个点用vector表示
{
double eps= 0;
for (int i=0;i<K;i++)
{
for (int j=0;j<clusters[i].size();j++)
{
eps+=EuclideanDistance(clusters[i][j],center[i]);//计算某个簇中的每个元组和对应的均值中心的距离，并累加
}
}
return eps;//返回累加的误差
}
void Print()//输出聚类结果
{
for(int l=0;l<K;l++)
{
printf("第%d个簇：\n",l+1);
vector<vector<double> >tmp=clusters[l];
for(int i=0;i<tmp.size();i++)
{
printf("%d: [",i+1);
for(int j=0;j<tmp[i].size()-1;j++)
{
printf("%lf,",tmp[i][j]);
}
printf("%lf]\n",tmp[i][tmp[i].size()-1]);//用向量的形式输出每个元组
}
}
}
void KMeans()
{
ChooseSeeds();//选择初始聚类中心
int lable=0;
//将元组根据初始聚类中心分类,lable表示均值中心标号
for(int i=0;i<number;i++)
{
lable=ClusterOfElement(node[i]);
clusters[lable].push_back(node[i]);
}
double OldErr=-1;
double NewErr=GetError();
while(fabs(NewErr-OldErr) >= 1e-9) //当两次迭代误差相差不到1时迭代终止
{
for(int i =0;i<K;i++) //更新每个簇的均值中心
{
center[i] =GetMeanValue(clusters[i]);
}
OldErr=NewErr;
NewErr=GetError(); //计算新的误差
for(int i=0;i<K;i++) //清空每个簇，因为之后要重新分配
{
clusters[i].clear();
}
//根据新的均值中心获得新的簇
for(int i=0;i<number;i++)
{
lable=ClusterOfElement(node[i]);
clusters[lable].push_back(node[i]);
}
}
printf("K-Means聚类结果：\n");
Print();
}
int main()
{
printf("Author: I-Hsin\n\n");
scanf("%d %d %d",&number,&K,&dimension);//输入输入数据总数，聚类个数，数据维数
//freopen("input.txt","r",stdin);//从 input.txt 文件读入输入数据
//freopen("I2000(2).txt","r",stdin);
//freopen("Test1.txt","r",stdin);
//freopen("Test2.txt","r",stdin);
freopen("Test3.txt","r",stdin);
for(int i=0;i<number;i++)
{
for(int j=0;j<dimension;j++)
{
double x=0;
scanf("%lf",&x);
node[i].push_back(x);
}
}
KMeans();
return 0;
}
//输入样例：
//4 2 2 （键盘读入）
//1 2 （文件读入）
//2 3
//444 555
//666 777