Contents
最近在填之前的坑，尝试用C++实现一个ISODATA的聚类分析算法，目前代码已经码完了，就慢慢的把文档给补上，记录一下自己零零碎碎做的工作；给自己立的flag画一个句号吧。
Kmeans聚类算法的原理
- 在分析一些数据前，我们需要选取通过一些方式对数据的种类进行划分，“物以类聚，人以群分”，分类是为了进一步的分析数据的性质。
- Kmeans聚类算法的一般步骤为
- 随机选取k个中心点
- 遍历所有样本，将数据划分到最近的中心点中（所谓近，即欧式距离近）
- 计算每个聚类的平均值，得到新的聚类中心
- 重复上述2-3步骤，直至所有中心点趋于收敛（变化幅度较小）或达到一定的迭代次数
- kmeans的原理简单，应用也很多，但是k的选取很重要，如果没有数据的先验知识，可能需要多次尝试不同的k值进行迭代计算，观察效果。
ISODATA聚类算法的原理
- 可以看作Kmeans的一种衍生，既然我们有时候只能大致确定k的范围，而不能确定k，那么可以加入分裂、合并机制，通过更自动化的机制，灵活的确定k的值
ISODATA算法的基本步骤
- 参数设置，需要设置的参数包括：
- nc : 初始聚类中心个数
- c : 预期的聚类个数
- tn : 每一类中允许的样本最少数目（少于该数目的聚类可能会被删除）
- te : 类内相对标准差上限（超过该上限的聚类，可能会被分裂）
- tc : 聚类中心点之间的最小距离
- nt : 每次迭代中最多可以合并的次数
- ns : 最多迭代次数
- 初始化
- 在所有样本中，随机选取nc个不重复样本作为聚类中心，之后根据距离最小法判断所有样本属于哪一种聚类
- 对第一次初始化后的聚类，检测是否符合tn条件，样本数目小于tn的类别被取消，重新根据上一步的距离最小法进行分配
- 计算分类后的聚类参数，包含以下参数：
- 聚类中心
- 各类中样本到聚类中心的平均距离
- 各个样本到其所属类别中心的总体平均距离
- 根据当前的状态选择下一步的行为，进行分裂、合并或停止
- 若迭代次数达到要求，则停止
- 若当前聚类数量小于期望数量的一般，则进行分裂检测，判断当前是否需要分裂，若需要则执行分裂行为
- 若聚类数量大于期望数量的两倍，则进行合并检测，判断是否需要合并，若需要则执行合并行为
- 若聚类数量在期望聚类数量的1/2到2倍之间，奇数次迭代则执行分裂检测，偶数次迭代则进行合并检测
关键步骤原理
- 分裂检测步骤
- 计算各个类别中，各个样本到类别中心的标准差，标准差为一矢量，求出每一类内标准差向量中的最大分量
- 若标准差向量中的最大分量大于te，同时满足以下条件之一，则对该类别进行分裂
- 聚类个数小于期望个数一半
- 聚类内部平均距离大于总体平均距离且聚类样本数量大于tn两倍
- 分裂操作：
- 在类内标准差最大的维度上分裂聚类中心点
- 利用距离最小法将当前聚类的样本点分配到两个聚类中
- 更新平均距离等聚类参数
- 合并检测步骤
- 计算聚类中心两两之间的距离，距离小于tc的聚类对，按照距离越小优先级越大的原则，进行合并
- 合并的总次数小于设定值，且一次迭代中，一个聚类只能被合并一次
- 更新聚类参数
代码实现
代码结构
- 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;
}
}
}
代码使用
- 只需要包含isodata一个头文件即可
- 我提供了一份数据和使用样例
- 使用前需要自定义读数据的函数，在common.cpp中，自定义read_data函数即可
- 新建isodata对象，初始化参数
- 调用成员函数run
- 具体如下所示
#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;
}
- 聚类分析的结果保存在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将聚类分析的结果绘制出来，结果如下
内容参考
最后
觉得文章对你有帮助的，欢迎扫码关注我的公众号，我会时不时的分享我的学习经验、收获。
