corClust.cpp

#include "corClust.h"
#include <algorithm>
#include <cmath>
#include <iostream>
#include <limits>
#include <numeric>
#include <queue>
#include <stack>
#include <utility>

// HierarchyTree实现
HierarchyTree::HierarchyTree(size_t id)
    : left(nullptr), right(nullptr), count(1) {
    pre_order_list.push_back(id);
}

HierarchyTree::HierarchyTree(std::shared_ptr<HierarchyTree> left,
                             std::shared_ptr<HierarchyTree> right)
    : left(left), right(right), count(0) {
    if (left) {
        count += left->count;
        auto left_pre = left->pre_order();
        pre_order_list.insert(pre_order_list.end(), left_pre.begin(),
                              left_pre.end());
    }
    if (right) {
        count += right->count;
        auto right_pre = right->pre_order();
        pre_order_list.insert(pre_order_list.end(), right_pre.begin(),
                              right_pre.end());
    }
}

std::shared_ptr<HierarchyTree> HierarchyTree::get_left() const { return left; }

std::shared_ptr<HierarchyTree> HierarchyTree::get_right() const {
    return right;
}

std::vector<size_t> HierarchyTree::pre_order() const { return pre_order_list; }

// 精确模拟SciPy的linkage函数行为
std::vector<std::vector<double>>
linkage(const std::vector<double> &condensed_distance) {
    // 计算原始维度 n
    double d = 1.0 + 8.0 * static_cast<double>(condensed_distance.size());
    size_t n = static_cast<size_t>((1.0 + std::sqrt(d)) / 2.0);

    // 检查输入有效性
    if (n < 2) {
        std::cerr << "错误：至少需要2个点进行聚类" << std::endl;
        return {};
    }

    // 复原距离矩阵（仅使用上三角部分）
    std::vector<std::vector<double>> distances(n, std::vector<double>(n, 0.0));
    size_t k = 0;
    for (size_t i = 0; i < n; ++i) {
        distances[i][i] = 0.0;
        for (size_t j = i + 1; j < n; ++j) {
            distances[i][j] = distances[j][i] = condensed_distance[k++];
        }
    }

    // 初始化活动集群
    std::vector<size_t> active_clusters(n);
    std::iota(active_clusters.begin(), active_clusters.end(), 0);

    // 结果矩阵Z
    std::vector<std::vector<double>> Z;
    Z.reserve(n - 1);

    // 存储每个节点的大小
    std::vector<size_t> cluster_sizes(n, 1);

    // 节点ID映射
    std::vector<double> node_ids(n);
    std::iota(node_ids.begin(), node_ids.end(), 0);

    // 处理合并直到只剩一个集群
    for (size_t iter = 0; iter < n - 1; ++iter) {
        // 查找最小距离
        double min_dist = std::numeric_limits<double>::max();
        size_t min_i = 0, min_j = 0;

        for (size_t i = 0; i < active_clusters.size(); ++i) {
            size_t ci = active_clusters[i];
            for (size_t j = i + 1; j < active_clusters.size(); ++j) {
                size_t cj = active_clusters[j];
                if (distances[ci][cj] < min_dist) {
                    min_dist = distances[ci][cj];
                    min_i = i;
                    min_j = j;
                }
            }
        }

        // 获取要合并的集群
        size_t cluster1 = active_clusters[min_i];
        size_t cluster2 = active_clusters[min_j];

        // 确保顺序（SciPy约定：总是较小的ID在前）
        if (node_ids[cluster1] > node_ids[cluster2]) {
            std::swap(cluster1, cluster2);
        }

        // 创建新的集群ID
        size_t new_idx = n + iter;

        // 添加Z矩阵行
        Z.push_back({node_ids[cluster1], node_ids[cluster2], min_dist,
                     static_cast<double>(cluster_sizes[cluster1] +
                                         cluster_sizes[cluster2])});

        // 计算新的距离（平均链接）
        std::vector<double> new_distances(new_idx, 0.0);

        for (size_t i : active_clusters) {
            if (i != cluster1 && i != cluster2) {
                double w1 = static_cast<double>(cluster_sizes[cluster1]);
                double w2 = static_cast<double>(cluster_sizes[cluster2]);
                new_distances[i] = (w1 * distances[cluster1][i] +
                                    w2 * distances[cluster2][i]) /
                                   (w1 + w2);
            }
        }

        // 更新距离矩阵
        distances.push_back(new_distances);
        for (size_t i = 0; i < new_idx; ++i) {
            distances[i].push_back(new_distances[i]);
        }

        // 更新集群大小
        cluster_sizes.push_back(cluster_sizes[cluster1] +
                                cluster_sizes[cluster2]);

        // 更新节点ID映射
        node_ids.push_back(static_cast<double>(new_idx));

        // 从活动集群中删除合并的集群
        active_clusters.erase(active_clusters.begin() + std::max(min_i, min_j));
        active_clusters.erase(active_clusters.begin() + std::min(min_i, min_j));

        // 添加新集群
        active_clusters.push_back(new_idx);
    }

    return Z;
}

// 从层次聚类结果构建树
std::shared_ptr<HierarchyTree>
to_tree(const std::vector<std::vector<double>> &Z, size_t n) {
    if (Z.empty()) {
        return nullptr;
    }

    // 创建叶节点
    std::vector<std::shared_ptr<HierarchyTree>> nodes;
    for (size_t i = 0; i < n; ++i) {
        nodes.push_back(std::make_shared<HierarchyTree>(i));
    }

    // 从Z矩阵中构建树
    for (const auto &cluster : Z) {
        int left_idx = static_cast<int>(cluster[0]);
        int right_idx = static_cast<int>(cluster[1]);

        // 确保索引有效
        if (left_idx < 0 || left_idx >= static_cast<int>(nodes.size()) ||
            right_idx < 0 || right_idx >= static_cast<int>(nodes.size())) {
            std::cerr << "错误：聚类索引超出范围 [" << left_idx << ", "
                      << right_idx << "]" << std::endl;
            return nullptr;
        }

        auto new_node =
            std::make_shared<HierarchyTree>(nodes[left_idx], nodes[right_idx]);
        nodes.push_back(new_node);
    }

    // 返回根节点（最后一个节点）
    return nodes.empty() ? nullptr : nodes.back();
}

// 模拟SciPy的flat_cluster函数，直接切割层次树
std::vector<std::vector<size_t>>
flat_cluster(const std::vector<std::vector<double>> &Z, size_t n,
             size_t maxClust) {
    if (Z.empty() || n == 0) {
        return {};
    }

    // 如果n <= maxClust，则每个点都是自己的集群
    if (n <= maxClust) {
        std::vector<std::vector<size_t>> result;
        for (size_t i = 0; i < n; ++i) {
            result.push_back({i});
        }
        return result;
    }

    // 集群标签，初始化为-1
    std::vector<int> labels(n + Z.size(), -1);

    // 当前集群ID
    int next_cluster = 0;

    // 检查每个链接
    for (size_t i = 0; i < Z.size(); ++i) {
        size_t node_id = n + i;

        // 如果合并的集群总大小 <= maxClust，则创建一个新集群
        if (static_cast<size_t>(Z[i][3]) <= maxClust) {
            // 分配一个新的集群ID给当前节点
            labels[node_id] = next_cluster++;
        } else {
            // 否则，传递任何现有子集群标签
            int left_id = static_cast<int>(Z[i][0]);
            int right_id = static_cast<int>(Z[i][1]);

            if (labels[left_id] >= 0) {
                labels[node_id] = labels[left_id];
            } else if (labels[right_id] >= 0) {
                labels[node_id] = labels[right_id];
            }
        }
    }

    // 如果根节点没有标签，分配一个新标签
    if (labels[n + Z.size() - 1] < 0) {
        labels[n + Z.size() - 1] = next_cluster++;
    }

    // 进行第二遍，将标签传播到叶节点
    std::function<void(int, int)> assign_labels = [&](int node, int label) {
        if (node < static_cast<int>(n)) {
            // 叶节点
            labels[node] = label;
        } else {
            // 内部节点
            int z_index = node - n;
            int left = static_cast<int>(Z[z_index][0]);
            int right = static_cast<int>(Z[z_index][1]);

            // 如果子节点没有标签，则分配当前标签
            if (labels[left] < 0) {
                assign_labels(left, label);
            }
            if (labels[right] < 0) {
                assign_labels(right, label);
            }
        }
    };

    // 从根节点开始分配标签
    assign_labels(n + Z.size() - 1, labels[n + Z.size() - 1]);

    // 收集每个集群的成员
    std::vector<std::vector<size_t>> clusters(next_cluster);
    for (size_t i = 0; i < n; ++i) {
        if (labels[i] >= 0 && labels[i] < next_cluster) {
            clusters[labels[i]].push_back(i);
        }
    }

    // 移除空集群
    clusters.erase(
        std::remove_if(clusters.begin(), clusters.end(),
                       [](const auto &cluster) { return cluster.empty(); }),
        clusters.end());

    return clusters;
}

// corClust实现
corClust::corClust(size_t n) : n(n), N(0) {
    // 初始化变量
    c.resize(n, 0.0);
    c_r.resize(n, 0.0);
    c_rs.resize(n, 0.0);
    C.resize(n, std::vector<double>(n, 0.0));
}

void corClust::update(const std::vector<double> &x) {
    if (x.size() != n) {
        throw std::invalid_argument("输入向量长度必须等于n");
    }

    N += 1;

    // 更新特征的累加和
    for (size_t i = 0; i < n; ++i) {
        c[i] += x[i];
    }

    // 计算残差
    std::vector<double> c_rt(n);
    for (size_t i = 0; i < n; ++i) {
        c_rt[i] = x[i] - c[i] / N;
        c_r[i] += c_rt[i];
        c_rs[i] += c_rt[i] * c_rt[i];
    }

    // 更新相关矩阵
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = 0; j < n; ++j) {
            C[i][j] += c_rt[i] * c_rt[j];
        }
    }
}

std::vector<std::vector<double>> corClust::corrDist() {
    std::vector<double> c_rs_sqrt(n);
    for (size_t i = 0; i < n; ++i) {
        c_rs_sqrt[i] = std::sqrt(c_rs[i]);
    }

    std::vector<std::vector<double>> D(n, std::vector<double>(n, 0.0));

    for (size_t i = 0; i < n; ++i) {
        for (size_t j = 0; j < n; ++j) {
            double denom = c_rs_sqrt[i] * c_rs_sqrt[j];
            // 防止除以零
            if (std::abs(denom) < 1e-100) {
                denom = 1e-100;
            }

            // 计算相关距离
            D[i][j] = 1.0 - C[i][j] / denom;

            // 修正可能出现的负值
            if (D[i][j] < 0) {
                D[i][j] = 0;
            }
        }
    }

    return D;
}

std::vector<std::vector<size_t>> corClust::cluster(size_t maxClust) {
    // 获取相关距离矩阵
    std::vector<std::vector<double>> D = corrDist();

    // 特殊情况处理：如果只有一个特征
    if (n <= 1) {
        if (n == 1) {
            return {{0}}; // 只有一个特征，直接返回
        }
        return {}; // 没有特征，返回空
    }

    // 提取上三角部分为压缩距离向量
    std::vector<double> condensed_D;
    condensed_D.reserve(n * (n - 1) / 2); // 预先分配内存
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = i + 1; j < n; ++j) {
            condensed_D.push_back(D[i][j]);
        }
    }

    // 执行层次聚类
    std::vector<std::vector<double>> Z = linkage(condensed_D);

    // 确保maxClust在合理范围内
    if (maxClust < 1) {
        maxClust = 1;
    }
    if (maxClust > n) {
        maxClust = n;
    }

    // 特殊情况：linkage返回空（可能因为错误）
    if (Z.empty()) {
        // 如果linkage失败，我们简单地将每个特征放在自己的聚类中
        std::vector<std::vector<size_t>> fallback_clusters;
        for (size_t i = 0; i < n; ++i) {
            fallback_clusters.push_back({i});
        }
        return fallback_clusters;
    }

    // 确定两种方法，以获得最好的一种
    auto tree_clusters = __breakClust__(to_tree(Z, n), maxClust);
    auto flat_clusters = flat_cluster(Z, n, maxClust);

    // Python版本使用的是层次树方法，我们在这里保持兼容
    return tree_clusters;
}

std::vector<std::vector<size_t>>
corClust::__breakClust__(const std::shared_ptr<HierarchyTree> &dendro,
                         size_t maxClust) {
    if (!dendro) {
        return {};
    }

    // 基本情况：找到一个最小聚类
    if (dendro->count <= maxClust) {
        return {dendro->pre_order()};
    }

    // 递归分解左右子树
    auto left_clusters = __breakClust__(dendro->get_left(), maxClust);
    auto right_clusters = __breakClust__(dendro->get_right(), maxClust);

    // 合并结果
    left_clusters.insert(left_clusters.end(), right_clusters.begin(),
                         right_clusters.end());
    return left_clusters;
}