最近在填之前的坑,尝试用C++实现一个ISODATA的聚类分析算法,目前代码已经码完了,就慢慢的把文档给补上,记录一下自己零零碎碎做的工作;给自己立的flag画一个句号吧。

Kmeans聚类算法的原理

  • 在分析一些数据前,我们需要选取通过一些方式对数据的种类进行划分,“物以类聚,人以群分”,分类是为了进一步的分析数据的性质。
  • Kmeans聚类算法的一般步骤为
    1. 随机选取k个中心点
    2. 遍历所有样本,将数据划分到最近的中心点中(所谓近,即欧式距离近)
    3. 计算每个聚类的平均值,得到新的聚类中心
    4. 重复上述2-3步骤,直至所有中心点趋于收敛(变化幅度较小)或达到一定的迭代次数
  • kmeans的原理简单,应用也很多,但是k的选取很重要,如果没有数据的先验知识,可能需要多次尝试不同的k值进行迭代计算,观察效果。

ISODATA聚类算法的原理

  • 可以看作Kmeans的一种衍生,既然我们有时候只能大致确定k的范围,而不能确定k,那么可以加入分裂、合并机制,通过更自动化的机制,灵活的确定k的值

ISODATA算法的基本步骤

  1. 参数设置,需要设置的参数包括:
    • nc : 初始聚类中心个数
    • c : 预期的聚类个数
    • tn : 每一类中允许的样本最少数目(少于该数目的聚类可能会被删除)
    • te : 类内相对标准差上限(超过该上限的聚类,可能会被分裂)
    • tc : 聚类中心点之间的最小距离
    • nt : 每次迭代中最多可以合并的次数
    • ns : 最多迭代次数
  2. 初始化
    1. 在所有样本中,随机选取nc个不重复样本作为聚类中心,之后根据距离最小法判断所有样本属于哪一种聚类
    2. 对第一次初始化后的聚类,检测是否符合tn条件,样本数目小于tn的类别被取消,重新根据上一步的距离最小法进行分配
  3. 计算分类后的聚类参数,包含以下参数:
    • 聚类中心
    • 各类中样本到聚类中心的平均距离
    • 各个样本到其所属类别中心的总体平均距离
  4. 根据当前的状态选择下一步的行为,进行分裂、合并或停止
    1. 若迭代次数达到要求,则停止
    2. 若当前聚类数量小于期望数量的一般,则进行分裂检测,判断当前是否需要分裂,若需要则执行分裂行为
    3. 若聚类数量大于期望数量的两倍,则进行合并检测,判断是否需要合并,若需要则执行合并行为
    4. 若聚类数量在期望聚类数量的1/2到2倍之间,奇数次迭代则执行分裂检测,偶数次迭代则进行合并检测

关键步骤原理

  • 分裂检测步骤
    1. 计算各个类别中,各个样本到类别中心的标准差,标准差为一矢量,求出每一类内标准差向量中的最大分量
    2. 若标准差向量中的最大分量大于te,同时满足以下条件之一,则对该类别进行分裂
      • 聚类个数小于期望个数一半
      • 聚类内部平均距离大于总体平均距离且聚类样本数量大于tn两倍
    3. 分裂操作:
      1. 在类内标准差最大的维度上分裂聚类中心点
      2. 利用距离最小法将当前聚类的样本点分配到两个聚类中
      3. 更新平均距离等聚类参数
  • 合并检测步骤
    1. 计算聚类中心两两之间的距离,距离小于tc的聚类对,按照距离越小优先级越大的原则,进行合并
    2. 合并的总次数小于设定值,且一次迭代中,一个聚类只能被合并一次
    3. 更新聚类参数

代码实现

代码结构

  • Cluster.h/cpp : 聚类结构体
  • common.h.cpp : 公共函数,包含一些矢量的操作符重载,主要为了编程时候的便捷
  • isodata.h/cpp : isodata聚类分析算法类
  • error.h : 定义了一些出错的输出字符,为了简便,会做一些出错的检查,然后报出错误,用于提示,但是本身不会做额外错误处理

主要类实现

  • Cluster
//Cluster.h
//
// Created by Jeff on 2019/1/6 0006.
//
#ifndef ISODATA_CLUSTER_H
#define ISODATA_CLUSTER_H
#include <vector>
#include <unordered_set>
using namespace std;
struct Cluster{
    double innerMeanDis; // 类内平均距离
    vector<double> sigma; // 每个聚类的标准差
    static double allMeanDis; // 总体平均距离
    vector<double> center; // 聚类中心位置的
    unordered_set<unsigned> ids; // 从属于此聚类的样本的id,位于isodata的data中的
    Cluster():
        center{}, innerMeanDis(0), sigma(vector<double>{}){}
    explicit Cluster(vector<double> &c):
        center(c), innerMeanDis(0), sigma(vector<double>(c.size(), 0)) {}
    void add_point(int p_index);
    void clear_ids();
};
#endif //ISODATA_CLUSTER_H
//Cluster.cpp
#include "Cluster.h"
#include "common.h"
double Cluster::allMeanDis = 0;
/**
 * 向聚类中添加点
 * @param p_index 点的id
 */
void Cluster::add_point(int p_index) {
    if (ids.find(static_cast<const unsigned int &>(p_index)) != ids.end())
    {
        cout << WARN_POINT_REPEAT << endl;
        return;
    }
    ids.emplace(p_index);
}
/**
 * 清除聚类中的所有点
 */
void Cluster::clear_ids() {
    ids.clear();
}
  • ISODATA
//
// Created by Jeff on 2019/1/2 0002.
//
#ifndef ISODATA_ISODATA_H
#define ISODATA_ISODATA_H
#include <vector>
#include <iostream>
#include <string>
#include <algorithm>
#include <unordered_set>
#include <time.h>
#include <random>
#include "error.h"
#include "Cluster.h"
#include <deque>
#include <functional>
#include "common.h"
using namespace std;
// 实现ISODATA聚类算法
class isodata {
private:
    typedef function<vector<vector<double>>(void)> READFUNC;
    // 直接初始化的数据
    // 为了简单,不预留更改设置的接口,只在初始化时设置
    unsigned _c; // 预期的聚类个数
    unsigned _nc; // 初始聚类中心数目
    unsigned _tn; // 每一类中允许的样本最少数目
    double _te; // 类内相对标准差上限
    double _tc; // 聚类中心点之间的最小距离
    unsigned _nt; // 每次迭代中最多可以合并的次数
    unsigned _ns; // 最多迭代次数
    // 不直接在构造函数中初始化
    unsigned row; // 数据的行数,也就是样本个数
    unsigned col; // 数据的列数,也就是特征个数
    vector<vector<double>> data; // 待分类的数据
    deque<Cluster> clusters; // 聚类
    READFUNC read_func; // 读取数据的函数,可以自定义
    double alpha; // 分裂系数
public:
    /**
     * 构造函数
     * @param _c 预期的聚类个数
     * @param _nc 初始聚类中心数目
     * @param _tn 每一类中允许的样本最少数目
     * @param _te 类内相对标准差上限
     * @param _tc 聚类中心点之间的最小距离
     * @param _nt 每次迭代中最多可以合并的次数
     * @param _ns 最多迭代次数
     * @param func 读取数据的函数,无输入,返回二维double vector
     */
    explicit isodata(unsigned int c, unsigned int _nc, unsigned int _tn,
                     double _te, double _tc, unsigned int _nt,
                     unsigned int _ns, READFUNC func) :
                     _c(c), _nc(_nc), _tn(_tn),
                     _te(_te), _tc(_tc), _nt(_nt),
                     _ns(_ns), row(0), col(0),
                     clusters(), read_func(std::move(func)), alpha(0.3) {
    }
    void run()
    {
        setData();
        init_clusters();
        for (int i = 0; i < _ns; ++i) {
            re_assign();
            check_tn();
            update_centers();
            update_meandis();
            switch_method(i);
        }
        output();
    }
private:
    void setData();
    void init_clusters();
    pair<int, double> get_nearest_cluster(int p_index, int ignore);
    pair<int, double> get_nearest_cluster(int p_index, unordered_set<unsigned>&);
    void re_assign();
    void check_tn();
    void update_centers();
    void update_center(Cluster& cluster);
    void update_sigmas();
    void update_sigma(Cluster& cluster);
    void update_meandis();
    void check_split();
    void split(const int& c_index);
    void check_merge();
    void merge(const int& id1, const int& id2);
    void switch_method(const int& index);
    void output() const;
};
#endif //ISODATA_ISODATA_H
//isodata.cpp
void isodata::setData()
{
    data = std::move(read_func());
    if (data.size() < _c || data.size() < _tn)
    {
        cout << WARN_DATA_SIZE << endl;
        return;
    }
    row  = static_cast<unsigned int>(data.size());
    col = static_cast<unsigned int>(data[0].size());
    for (vector<double> & d : data)
    {
        if (d.size() != col)
        {
            cout << WARN_DATA_SIZE << endl;
            data.clear();
            return;
        }
    }
}
/**
 * 初始化 随机选取_nc个聚类中心
 */
void isodata::init_clusters() {
    // 选取随机id
    std::default_random_engine rand(static_cast<unsigned int>(time(nullptr)));
    std::uniform_int_distribution<unsigned> rnd(0, row-1);
    unordered_set<unsigned> ids;
    while (ids.size() < _nc)
    {
        auto id = rnd(rand);
        // 避免数据之间有重复的 选取出重复的中心点坐标
        bool flg = false;
        for (auto &item : ids)
        {
            if (get_distance(data[item], data[id]) == 0.0)
                flg = true;
        }
        if (!flg)
            ids.emplace(id);
    }
    // 初始化聚类
    for (auto &id : ids)
    {
        clusters.emplace_back(Cluster{this->data[id]});
    }
}
/**
 * 根据距离聚类中心的最小距离 选取当前点所属的聚类
 * @p_index 当前点在data中的序号
 * @ignore 设置忽略序号, 缺省为-1, 主要用于删除聚类时使用
 * @return 序号+距离
 */
pair<int, double> isodata::get_nearest_cluster(int p_index, int ignore = -1) {
    if (ignore != -1 && clusters.size() == 1)
        cout << WARN_CLUSTER_SIZE_SMALL << endl;
    int c_index = 0;
    double dis(get_distance(data[p_index], clusters[c_index].center));
    for (int i = 1; i < clusters.size(); ++i) {
        if (ignore != -1 && i == ignore)
            continue;
        auto d = get_distance(data[p_index], clusters[i].center);
        if (d < dis)
        {
            dis = d;
            c_index = i;
        }
    }
    return {c_index, dis};
}
/**
 * 根据距离聚类中心的最小距离 选取当前点所属的聚类
 * @param p_index
 * @param cluster_ids 一个集合,集合内的聚类中心不做考虑
 * @return 序号+距离
 */
pair<int, double> isodata::get_nearest_cluster(int p_index, unordered_set<unsigned> &cluster_ids) {
    // 选择第一个不在cluster_ids中的聚类
    int c_index = 0;
    while (cluster_ids.find(static_cast<const unsigned int &>(c_index)) != cluster_ids.end())
        ++c_index;
    //初始一个距离
    double dis(get_distance(data[p_index], clusters[c_index].center));
    for (int i = c_index+1;
    i < clusters.size() && (cluster_ids.find(static_cast<const unsigned int &>(i)) == cluster_ids.end());
    ++i)
    {
        auto d = get_distance(data[p_index], clusters[i].center);
        if (d < dis)
        {
            dis = d;
            c_index = i;
        }
    }
    return {c_index, dis};
}
/**
 * 依据距离最小原则重新分配点
 */
void isodata::re_assign() {
    for (auto &cluster : clusters) {
        cluster.clear_ids();
    }
    for (int i = 0; i < row; ++i) {
        auto&& res = get_nearest_cluster(i);
        clusters[res.first].add_point(i);
    }
}
/**
 * 检测每个聚类中的个数是否少于_tn,如果少于则取消此类别
 */
void isodata::check_tn() {
    unordered_set<unsigned> to_erase;
    for (int i = 0; i < clusters.size(); ++i) {
        if (clusters[i].ids.size() >= _tn)
            continue;
        to_erase.emplace(i);
        auto &indexs = clusters[i].ids;
        for (unsigned int index : indexs) {
            auto &&res = get_nearest_cluster(index, to_erase);
            clusters[res.first].add_point(index);
        }
        clusters[i].ids.clear();

    }
    for (auto it = clusters.begin(); it != clusters.end();)
    {
        if (it->ids.empty())
            it = clusters.erase(it);
        else
            ++it;
    }
}
/**
 * 更新各个聚类的中心坐标
 */
void isodata::update_centers() {
    for (auto &cluster : clusters)
    {
        update_center(cluster);
    }
}
/**
 * 更新某个聚类的中心坐标
 */
void isodata::update_center(Cluster& cluster) {
    vector<double> sum(col, 0);
    for (auto &index : cluster.ids)
    {
        sum += data[index];
    }
    cluster.center = sum/ static_cast<double>(cluster.ids.size());
}
/**
 * 更新各个聚类的标准差
 */
void isodata::update_sigmas() {
    for ( auto& cluster : clusters) {
        update_sigma(cluster);
    }
}
/**
 * 更新单个聚类的标准差
 */
void isodata::update_sigma(Cluster& cluster) {
    cluster.sigma.clear();
    auto &sigma = cluster.sigma;
    sigma.resize(col, 0);
    for (auto& id : cluster.ids) {
        sigma += mypow(cluster.center - data[id], 2);
    }
    auto &&res(vector_sqrt(sigma));
    swap(res, sigma);
}
/**
 * 更新平均距离
 */
void isodata::update_meandis() {
    Cluster::allMeanDis = 0;
    for (auto &&cluster : clusters) {
        double dis(0.0);
        for (auto &id : cluster.ids) {
            dis += get_distance(data[id], cluster.center);
        }
        Cluster::allMeanDis += dis;
        cluster.innerMeanDis = dis/cluster.ids.size();
    }
    Cluster::allMeanDis /= row;
}
/**
 * 检测是否需要分裂,如果需要则执行分裂操作
 */
void isodata::check_split() {
    update_sigmas();
    while (true)
    {
        bool flag = false;
        for (unsigned j = 0; j < clusters.size(); ++j) {
            auto& cluster = clusters[j];
            for (int i = 0; i < col; ++i) {
                if (cluster.innerMeanDis > Cluster::allMeanDis) {
                    if (cluster.sigma[i] > _te && (cluster.ids.size() > 2 * _tn + 1 || clusters.size() < _c / 2)) {
                        flag = true;
                        split(j);
                    }
                } else {
                    if (cluster.sigma[i] > _te && clusters.size() < _c / 2) {
                        flag = true;
                        split(j);
                    }
                }
            }
        }
        if (!flag)
            break;
        else
            update_meandis();
    }
}
/**
 * 分裂第c_index个聚类
 * @param c_index
 */
void isodata::split(const int &c_index) {
    //根据标准差选取分裂维度
    auto& cluster = clusters[c_index];
    auto iter = max_element(cluster.sigma.begin(), cluster.sigma.end());
    auto pos = distance(cluster.sigma.begin(), iter);
    Cluster newcluster(cluster.center);
    //分裂
    newcluster.center[pos] -= alpha*cluster.center[pos];
    cluster.center[pos] += alpha*cluster.center[pos];
    for (const auto &id : cluster.ids) {
        auto d1 = get_distance(data[id], cluster.center);
        auto d2 = get_distance(data[id], newcluster.center);
        if (d2 < d1)
            newcluster.ids.emplace(id);
    }
    unordered_set<unsigned> dids;
    set_difference(cluster.ids.begin(), cluster.ids.end(),
                   newcluster.ids.begin(), newcluster.ids.end(),
                   inserter(dids, dids.begin()));
    swap(dids, cluster.ids);
    // 更新参数并保存新聚类
    update_center(newcluster);
    update_sigma(newcluster);
    update_center(cluster);
    update_sigma(cluster);
    clusters.emplace_back(newcluster);
}
/**
 * 检测是否需要合并,如果需要则执行合并操作
 */
void isodata::check_merge() {
    // 自定义数据类型
    typedef pair<pair<unsigned, unsigned>, double> UNIT;
    //1 计算聚类中心两两之间的距离,保留小于_tc的,并从小到达排序
    vector<UNIT> uvec;
    for (unsigned i = 0; i < clusters.size(); ++i) {
        for (unsigned j = i+1; j < clusters.size(); ++j) {
            auto dis = get_distance(clusters[i].center, clusters[j].center);
            if (dis < _tc)
                uvec.emplace_back(UNIT({i,j}, dis));
        }
    }
    sort(uvec.begin(), uvec.end(), [](UNIT& left, UNIT& right){ return left.second < right.second;});
    //2 执行合并操作
    unordered_set<unsigned> clusterIds{};
    unsigned cnt(0);
    for (const auto &unit : uvec) {
        auto& cids = unit.first;
        auto& cdis = unit.second;
        if (clusterIds.find(cids.first) == clusterIds.end() &&
            clusterIds.find(cids.second) == clusterIds.end())
        {
            merge(cids.first, cids.second);
            clusterIds.emplace(cids.first);
            clusterIds.emplace(cids.second);
        }
        // 检测是否超过最大次数
        if (++cnt > _ns)
            break;
    }
    //3 清除被合并的聚类
    for (auto iter = clusters.begin(); iter != clusters.end();)
    {
        if (iter->ids.empty())
            iter = clusters.erase(iter);
        else
            iter++;
    }
}
/**
 * 合并操作 合并id1 id2的聚类
 * @param id1
 * @param id2
 */
void isodata::merge(const int &id1, const int &id2) {
    auto &c1 = clusters[id1];
    auto &c2 = clusters[id2];
    c1.center = (c1.center* c1.ids.size() + c2.center*c2.ids.size())/(c1.ids.size()+c2.ids.size());
    c2.ids.clear();
}
/**
 * 根据情况选择下一步操作
 */
void isodata::switch_method(const int& index) {
    if (index == _ns-1)
    {
        _tc = 0;
        check_merge();
    } else if (clusters.size() <= _c/2)
    {
        check_split();
    } else if (index % 2 == 0 || clusters.size() >= 2*_c)
    {
        check_merge();
    } else
    {
        check_split();
    }
}
/**
 * 输出聚类分析的结果
 */
void isodata::output() const {
    // 在命令行窗口打印
    cout << "Original Data Number : " << row << endl;
    cout << "Cluster Number : " << clusters.size() << endl;
    for (int i = 0; i < clusters.size(); ++i) {
        cout << "Number " << to_string(i+1) << " : " << clusters[i].ids.size() << endl;
        cout << "cluster center : " << clusters[i].center << endl;
    }
    // 输出到文件
    ofstream res(R"(E:\CPP\Clion\ISODATA\clusters.txt)");
    if (!res.is_open())
        cout << WARN_FILE_OPEN_FAIL << endl;
    res << clusters.size() << endl;
    for (int j = 0; j < clusters.size(); ++j) {
        res << j+1 << " " << clusters[j].ids.size() << endl;
        for (const auto &id : clusters[j].ids) {
            res << data[id] << endl;
        }
    }
}

代码使用

  1. 只需要包含isodata一个头文件即可
  2. 我提供了一份数据和使用样例
  3. 使用前需要自定义读数据的函数,在common.cpp中,自定义read_data函数即可
  4. 新建isodata对象,初始化参数
  5. 调用成员函数run
  6. 具体如下所示
#include "isodata.h"
#include "MyTime.h"
/**
 * 在运行前,需要自定义common.cpp中的read_data函数
 * @return
 */
int main() {
    CMyTimeWrapper c;
    c.tic();
    isodata isodata1(4, 90, 10, 90, 20, 5, 500, read_data);
    isodata1.run();
    c.tocMs();
    return 0;
}
  1. 聚类分析的结果保存在txt文件中,具体见isodata中的output函数,命令行输出的结果如下所示
Original Data Number : 927
Cluster Number : 4
Number 1 : 25
cluster center : 4.77319,16.0967
Number 2 : 604
cluster center : 30.4196,-2.65259
Number 3 : 277
cluster center : 33.988,24.7817
Number 4 : 21
cluster center : -18.4698,-33.7524
cost time: 239.088 ms

代码地址

聚类效果分析

  • 利用python将聚类分析的结果绘制出来,结果如下
    聚类分析效果图

内容参考

最后

觉得文章对你有帮助的,欢迎扫码关注我的公众号,我会时不时的分享我的学习经验、收获。