样本熵(SampEn)的C/C++代码实现与优化
本文不介绍什么是样本熵，具体推荐看此文https://blog.csdn.net/Cratial/article/details/79742363，写的很好，里面的示例也被我拿来测试代码写的对不对。
本文所以代码可以在此处找到https://gitcode.net/PeaZomboss/miscellaneous，文件夹231424-sampen。
前言
一开始有人找我帮忙写些C++代码，实现一些算法，然后就帮忙写了，其中就有一个样本熵的。
当时粗写了一个，后来似乎性能不太够，当然实际上性能是没问题的，虽然确实不如现在优化的，但也不会很差，原因是代码被魔改了一部分，当然这个和算法的优化没关系，就不展开了。
不过既然都要优化了，那不如在算法层面也做一些改进，让速度再快一点，于是就想了想办法提升了性能。于是在此之后就整理了一下代码，并严谨测试，然后与大家分享。
开头提到的那篇文章里有用python写的代码，不过糖比较多让人看了比较乱，还有一个matlab写的。但是不管是python还是matlab都不是拿来实际应用的，一般追求性能的模块会用C/C++，当然用Java应该也不会太差。
代码实现
首先是一般情况下很容易想到的代码：
static double step(double *X, int N, int m, double r)
{
double sum = 0;
for (int i = 0; i <= N - m; i++) {
int Bi = 0;
for (int j = 0; j <= N - m; j++) {
if (i != j) {
/* 找出最大差值 */
double D = fabs(X[i] - X[j]);
for (int k = 1; k < m; k++) {
double t = fabs(X[i + k] - X[j + k]);
if (D < t)
D = t;
}
if (D <= r)
Bi++;
}
}
sum += 1.0 * Bi / (N - m); // 会有累加误差
}
return sum / (N - m + 1);
}
double SampEn(double *X, int N, int m, double r)
{
double B = step(X, N, m, r);
if (B == 0) // 尽管大部分时候不会是0
return 0;
double A = step(X, N, m + 1, r);
if (A == 0)
return 0;
return -log(A / B);
}
这里的所有步骤基本都是按照原算法描述的内容进行编写的，所以非常好理解啊。
不过这个step函数的sum += 1.0 * Bi / (N - m);在N较大的时候累加误差会增加，而实际上这个Bi变量完全可以累加到最后在进行运算，不但减少误差还能降低运算量，所以可以改成：
static double step(double *X, int N, int m, double r)
{
int Bi = 0;
for (int i = 0; i <= N - m; i++) {
for (int j = 0; j <= N - m; j++) {
if (i != j) {
double D = fabs(X[i] - X[j]);
for (int k = 1; k < m; k++) {
double t = fabs(X[i + k] - X[j + k]);
if (D < t)
D = t;
}
if (D <= r)
Bi++;
}
}
}
return 1.0 * Bi / (N - m) / (N - m + 1);
}
这样子几乎就是按照原来的算法思路进行了一些简单处理，但提升并不明显。
代码优化
仔细看这个step()函数，那复杂度可是O(n^2)啊，而且还要跑两次，那损耗肯定惊人。但是仔细看这两次计算，只有第一次的m变成了m+1，那能不能把两次运算整合到一起呢？
注意到m+1和m的区别在于循环次数少了一次，寻找最大差值的向量多了一个，其余是没有区别的，那我们完全可以在处理m的时候顺便把m+1的情况也放一起就行了。
请看：
double FastSampEn(double *X, int N, int m, double r)
{
int Ai = 0, Bi = 0;
int LoopsSub1 = N - m; // 循环次数减一，因为算法中的表述是从1到N-m+1，而我们是从0开始的
for (int i = 0; i <= LoopsSub1; i++) {
for (int j = 0; j <= LoopsSub1; j++) {
if (i != j) {
// 这里一样找m个的最大差值
double D = fabs(X[i] - X[j]);
for (int k = 1; k < m; k++) {
double t = fabs(X[i + k] - X[j + k]);
if (D < t)
D = t;
}
if (D <= r)
Bi++;
// 对于m+1的情况，当到达m维数的边界的时候显然是不行的
// 所以我们只要限制边界情况就行了
if (i != LoopsSub1 && j != LoopsSub1) {
double t = fabs(X[i + m] - X[j + m]);
if (D < t) // 判断最后一个是不是最大值
D = t;
if (D <= r)
Ai++;
}
} // i!=j
} // j
} // i
double B = 1.0 * Bi / (N - m) / (N - m + 1);
double A = 1.0 * Ai / (N - m - 1) / (N - m);
if (B == 0 || A == 0)
return 0;
return -log(A / B);
}
这样修改以后速度可以达到原先的1.5倍，算上循环内部的开销，达到这样也是挺不错的。
至此，针对样本熵代码的基本优化就到这了。当然一定还会有更好的方案，但是碍于水平有限，目前只能如此。
其实方法是有的，比如暂存计算的结果，因为有一半是可以复用的，或者别的办法，不过这样空间复杂度会比较高，而目前的实现几乎没有额外的内存开销。
不过，考虑到实际应用情况，这个m大多数时候都是取值2，所以针对m为2的情况又专门优化了一下：
double FastSampEn_m2(double *X, int N, double r)
{
int Ai = 0, Bi = 0;
int LoopsSub1 = N - 2;
for (int i = 0; i <= LoopsSub1; i++) {
for (int j = 0; j <= LoopsSub1; j++) {
if (i != j) {
double D = fabs(X[i] - X[j]);
double t = fabs(X[i + 1] - X[j + 1]);
if (D < t)
D = t;
if (D <= r)
Bi++;
if (i != LoopsSub1 && j != LoopsSub1) {
double t = fabs(X[i + 2] - X[j + 2]);
if (D < t)
D = t;
if (D <= r)
Ai++;
}
}
}
}
double B = 1.0 * Bi / (N - 2) / (N - 1);
double A = 1.0 * Ai / (N - 3) / (N - 2);
if (B == 0 || A == 0)
return 0;
return -log(A / B);
}
这样算是进一步压榨了CPU。
测试代码
接下来为了测试上述代码是否正确，就用开头提到的那篇文章的数据进行测试看看：
std::vector<double> x;
for (int i = 0; i < 17; i++) {
x.push_back(85);
x.push_back(80);
x.push_back(89);
}
std::cout << SampEn(x.data(), x.size(), 2, 3) << '\n';
std::cout << FastSampEn(x.data(), x.size(), 2, 3) << '\n';
std::cout << FastSampEn_m2(x.data(), x.size(), 3) << '\n';
原文给出的结果是0.0008507018803128114，不过由于std::cout的舍入问题，实际输出的是0.000850702，是符合结果的。
然后再测试一下性能：
for (int i = 0; i < 10000; i++) {
x.push_back(85);
x.push_back(80);
x.push_back(89);
}
double se;
clock_t t;
t = clock();
se = SampEn(x.data(), x.size(), 2, 3);
t = clock() - t;
std::cout << se << ", time = " << t << " ms\n";
t = clock();
se = FastSampEn(x.data(), x.size(), 2, 3);
t = clock() - t;
std::cout << se << ", time = " << t << " ms\n";
t = clock();
se = FastSampEn_m2(x.data(), x.size(), 3);
t = clock() - t;
std::cout << se << ", time = " << t << " ms\n";
先是塞了30000个数据，然后用clock()函数计时，依次测试每一个函数，然后看看结果是否一致，时间差距有多少。
如果有兴趣的话，可以去用之前提到那篇文章的python代码跑跑看，30000个数据花了我半个多小时才跑完，而我的这个代码即使没经过优化也不超过一分钟。而matlab的代码我没试过，电脑上没这个软件，也不是学这个的。
完整代码在开头给出的仓库里，打开test.cpp就可以看到。
运行测试
编译器和版本：g++ (x86_64-win32-seh-rev0, Built by MinGW-W64 project) 8.5.0
测试机器CPU：AMD Ryzen 5 4600H
- 编译命令
g++ -D ALL_IN_ONE test.cpp -o test_x64_O0 -O0
0.000850702
0.000850702
0.000850702
2.21505e-09, time = 15275 ms
2.21505e-09, time = 9409 ms
2.21505e-09, time = 8020 ms
- 编译命令
g++ -D ALL_IN_ONE test.cpp -o test_x64_O1 -O1
0.000850702
0.000850702
0.000850702
2.21505e-09, time = 2512 ms
2.21505e-09, time = 1604 ms
2.21505e-09, time = 1151 ms
- 编译命令
g++ -D ALL_IN_ONE test.cpp -o test_x64_O2 -O2
0.000850702
0.000850702
0.000850702
2.21505e-09, time = 2654 ms
2.21505e-09, time = 1604 ms
2.21505e-09, time = 1143 ms
- 编译命令
g++ -D ALL_IN_ONE test.cpp -o test_x64_O3 -O3
0.000850702
0.000850702
0.000850702
2.21505e-09, time = 2656 ms
2.21505e-09, time = 1603 ms
2.21505e-09, time = 1145 ms
- 编译命令
g++ -D ALL_IN_ONE test.cpp -o test_x64_Os -Os
0.000850702
0.000850702
0.000850702
2.21505e-09, time = 3058 ms
2.21505e-09, time = 1631 ms
2.21505e-09, time = 1167 ms
可以看到开了-O1及以后差距已经不明显了，说明后面的两段代码对编译器是比较友好的。
32位的就不贴了，性能肯定是不如64位的，但是依然可以看到的是即使是开了编译器优化，性能提升也是非常显著的。
不过32位有个有意思的情况，就是用-Os优化的性能居然比-O3优化的高出一截，这点是挺有意思的，因为不管开哪个优化32位默认都是用x87 FPU进行运算的。
另一个有意思的是如果编译器指令加上-msse2 -mfpmath=sse之后，不论开-Os还是-O3性能都是差不多的了，而且和性能64位也是基本上一样的了，毕竟64位默认就用的SSE2 FPU。
以上都是用g++测试的，如果用vs的话结果也是差不多的。
编译方案
对于x86来说，使用SSE2 FPU肯定是最好的，因为上述代码涉及大量double类型浮点数的运算，所以可以获得最佳性能。至于说什么CPU支持SSE2，那就这么说吧，20年前(以2023年为基准)的新CPU大部分都支持，因为SSE2是Intel在奔腾4(2000年)加入的，而AMD也在K8(2003年)加入了。因此可以在开编译器优化的基础上选择让编译器生成SSE2的代码。至于具体怎么操作，各家编译器有自己的方法，这个去搜一下就行了。
而对于x64来说，默认就会用SSE2，所以只要开编译器优化就行了。
如果用vs的话似乎x86默认就是开启SSE2的，而且也不用管优化，直接Release模式就行了。
标签：SampEn,0.000850702,09,int,double,2.21505,样本,C++,ms From： https://www.cnblogs.com/PeaZomboss/p/17074163.html