通过K-MEDOIDS算法对时间序列进行聚类的实现

最近做数据挖掘相关的工作，题目是时间序列聚类研究，目前对于这方面的研究都还只是在起步阶段，被广泛使用的还是基于K-MEDOIDS的聚类，放弃K-MEANS的主要原因还是时间序列之间序列的计算难度，对于这方面我们也已经有了一定的进展，不过也还是有很多的问题。

把基于DTW与K-MEDOIDS的时间序列聚类的算法贴出来，希望对大家有些帮助吧。

这份代码是我在以前的代码的基础上直接改的，所以C和C++有些混用。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
using namespace std;
#define NA  60      /* 数据维数 */
#define K   6       /* 聚类数 */
#define Psize   36      /* 种群大小 */
#define T   30      /* 最大迭代数 */
#define ED  0.0000001   /* 结束条件 */
#define Min 1000000   /*最小值*/
#define MinCmp(a,b) (a<b?a:b)
#define INF 300000000
//记录每个点的坐标已经到K个中心点的距离
typedef struct {
    double p[NA];
    double distance[K];
}Point;
//记录整个种群聚类的相关信息
typedef struct {
    Point clu_cent[K];  /* 即cluster_center 簇类中心 */
    int cluster[K][Psize];  /* 簇类数组 */
    int cluster_num[K]; /* 簇类中一组数据的编号 */
    double fitness;     /* 样本适应度值,用于判断结束条件 */
    double old_fitness; /* 前一次迭代的适应度值 */
    double Je;      /* 所有样本的平方误差和 */
}Pop;

/* 声明函数 */
int Is_equal(int a[], int n, int b);
double Euclid1(double x, double y);
double dtw(int x, int y);
void input_data();
void Init_center();
void calculate_distance();
void Make_new_cluster();
void Make_new_center();
void output_info(int flag);
Point all_data[Psize];      /* 数据大小 */
Pop pop;
/************************************************
 * 从外部文件导入数据，对于没有数据文件将报错， *
 * 数据文件的格式根据 NA 决定，例如NA = 4时，测 *
 * 试数据为四维，则test.data 为：     *
 *  1   2   3   4       *
 *  1.0 1.2 1.3 1.4     *
 *      ......              *
 *      ......              *
 ***********************************************/
double Dtwdistance(Point x, Point y)
    {
    double distance[NA+1][NA+1];
    double output[NA+1][NA+1];
    int i,j;
    memset(distance,0,sizeof(distance));
    memset(output,0,sizeof(output));
    for (i=1;i<=NA;i++)
        for (j=1;j<=NA;j++)
            distance[i][j]=Euclid1(x.p[i-1],y.p[j-1]);
    output[1][1]=distance[1][1];
    for (i=0;i<=NA;i++)
        {
        output[i][0]=INF;
        output[0][i]=INF;
        }
    for (i=2;i<=NA;i++)
        output[i][1]=output[i-1][1]+distance[i][1];
    for (i=2;i<=NA;i++)
        output[1][i]=output[1][i-1]+distance[1][i];
    for(i=2;i<=NA;i++)
        for(j=2;j<=NA;j++)
            output[i][j]=Mintt(Mintt(output[i-1][j-1],output[i][j-1]),output[i-1][j])+distance[i][j];
    //DP过程，计算DTW距离
    return output[NA][NA];
    }

double dtw(int a, int b)
    {
    if (a==b)
    return 0;
    return Dtwdistance(all_data[a],all_data[b]);
    }

void input_data()
{
        int i,j,temp;
        freopen("data.txt","r",stdin);
        for(i = 0; i < Psize; i++)
        {
            cin>>temp;
            for (j=0;j<NA;j++)
            {
            cin>>all_data[i].p[j];
            }
        }
}

/***************************************************
 * 随机初始化聚类质心，当随机数有相同时跳过继续执行*
 **************************************************/
void Init_center()
{
    int i, j;
    int num = 0;
    int rand_num;
    int rand_num_tmp[K];
    /* 随机产生三个0～Psize的数 */
    while(num < K){
        rand_num = rand() % Psize;
        if(!Is_equal(rand_num_tmp, num, rand_num))
            rand_num_tmp[num++] = rand_num;
    }
    freopen("output.txt","w",stdout);
    cout<<"初始化重心为:"<<endl;
    for(i = 0; i < K; i++){
    cout<<rand_num_tmp[i]<<endl;
        for(j = 0; j < NA; j++)
        {
                cout<<all_data[rand_num_tmp[i]].p[j]<<'/t';
                pop.clu_cent[i].p[j] = all_data[rand_num_tmp[i]].p[j];
        }
        cout<<endl;
    }
}
/**********************************
 * 检查数据是否有相等,相等返回1 *
 *********************************/
int Is_equal(int a[], int n, int b)
{
    int i;
    for(i = 0; i < n; i++)
        if(a[i] == b) return 1;
    return 0;
}

/********************************************
 * 计算Psize组数据到K个质心的欧几里德距离*
 *******************************************/
void calculate_distance()
{
    int i, j;
    for(i = 0; i < Psize; i++)
        for(j = 0; j < K; j++){
            all_data[i].distance[j] = dtw(i, j);
        }
}

/************************************************
 * 此函数为欧几里德距离公式函数,此处用于计算*
 * 一组数据到对应簇中心的欧几里德距离。   *
 ***********************************************/
double Euclid1(double x, double y)
    {
    return (y-x)*(y-x);
    }
/************************
 * 将数据进行簇集归类 *
 ***********************/
void Make_new_cluster()
{
    int i, j;
    double min;

    for(i = 0; i < K; i++)       /* 初始化编号 */
        pop.cluster_num[i] = 0;
    for(i = 0; i < Psize; i++){
        int index = 0;
        min = all_data[i].distance[0];
        for(j = 1; j < K; j++){      /* 筛选到簇心欧几里德最小的 */
            if(all_data[i].distance[j] < min){
                min = all_data[i].distance[j];
                index = j;
            }
        }
        /* 划分簇集 */
        pop.cluster[index][pop.cluster_num[index]++] = i;
    }
    /* 计算所有样本的平方误差和 */
    pop.Je = 0;
    for(i = 0; i < K; i++)
        for(j = 0; j < pop.cluster_num[i]; j++){
        /* 样本到簇心的值即为其欧几里德距离 */
            pop.Je +=pow(all_data[pop.cluster[i][j]].distance[i],2);
        }
    pop.old_fitness = pop.fitness;  /* 前一次迭代适应度值 */
//  printf("old_fitness = %lf/n", pop.old_fitness);
    pop.fitness = pop.Je;   /* 所有样本平方误差和即为适应度值 */
}
/*************************************************
 * 更新簇心，即求其群类的平均距离为新的簇心  *
 ************************************************/
void Make_new_center()
{
    int i, j, t,n;
    double tmp_sum;
    int index;
    double min=Min;
    for(i = 0; i < K; i++)
    {
        tmp_sum = 0;
        for(j= 0;j< pop.cluster_num[i];j++)
        for (t=0;t<pop.cluster_num[i];t++)
        {
            tmp_sum =dtw(pop.cluster[i][j],pop.cluster[i][t]);
            if(tmp_sum<min)
                min=tmp_sum;
            index=j;
        }
        for (n=0;n<NA;n++)
            pop.clu_cent[i].p[n]=all_data[pop.cluster[i][index]].p[n];

    }
}

/********************************
 *  输出信息函数      *
 * 显示格式为：       *
 * 质心K:         *
 * NA维的质心数据     *
 * 簇类K：         *
 * NA维属于簇类K的数据      *
 *  ......          *
 *  ......          *
 *******************************/
void output_info(int flag)
{
    int i, j, n;
    for(i = 0; i < K; i++){
        if(flag == 0){
            cout<<"初始化质心:"<<i+1<<'/t'<<pop.cluster_num[i]<<endl;
            for(n = 0; n < NA; n++)
                cout<<pop.clu_cent[i].p[n]<<'/t';
            cout<<endl;
        }
        else if(flag == 1){
            cout<<"最终质心:"<<i+1<<endl;
            for(n = 0; n < NA; n++)
                cout<<pop.clu_cent[i].p[n]<<'/t';
            cout<<endl;
        }
        cout<<"簇类:" <<i + 1<<endl;
        for(j = 0; j < pop.cluster_num[i]; j++){
            cout<<pop.cluster[i][j]<<'/t';
            for(n = 0; n < NA; n++)
                cout<<all_data[pop.cluster[i][j]].p[n]<<'/t';
                cout<<endl;
            }
    }
}
/********************************
 *      主函数     *
 *******************************/
int main()
{
    int i,j;
    double differ = 1;  /* 适应度差值 */
    int flag = 0;       /* 用来标示是显示初始数据还是聚类后的数据 */
    freopen("data.txt","r",stdin);
    freopen("output.txt","w",stdout);
    input_data();       /* 导入数据 */
    Init_center();      /* 初始化质心 */
    for(i = 0; (i < T) && (differ > ED); i++){
        calculate_distance();       /* 计算欧几里德距离 */
        Make_new_cluster();     /* 产生新的聚类 */

        if(flag == 0){
            output_info(flag);  /* 输出第一次聚类后的结果 */
            flag = 1;
        }
        Make_new_center();      /* 对新的聚类产生新的质心 */
        differ = pop.old_fitness - pop.fitness; /* 判断条件 */
        differ = fabs(differ);
    //  printf("differ = %lf/n", differ);
    //  printf("old_fitness = %lf/n", pop.old_fitness);
    //  printf("fitness = %lf/n", pop.fitness);
    }
    printf("+++++++++++++++++++++++++++++++++++++++++++++++++++++++/n");
    output_info(flag);  /* 聚类后显示结果 */
    return 0;
}