// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// DBSCAN.hpp — 密度ベースクラスタリング (DBSCAN)
//
// 移植元: lib++20 機械学習/クラスタリング/DBSCAN.h
// 改良点:
//   - calx Matrix<T> ベース
//   - ノイズ = -1, クラスタ = 0, 1, 2, ...
//   - デバッグ出力除去

#ifndef CALX_ML_CLUSTERING_DBSCAN_HPP
#define CALX_ML_CLUSTERING_DBSCAN_HPP

#include <math/core/matrix.hpp>
#include <cmath>
#include <vector>

namespace calx::ml {

    static constexpr int DBSCAN_NOISE = -1;

    /// DBSCAN クラスタリング
    /// @param data     (n × dim) サンプル行列
    /// @param eps      近傍半径
    /// @param minPts   コア点の最低近傍数
    /// @return         各サンプルのクラスタ番号 (ノイズ = -1)
    template<typename T>
    std::vector<int> dbscan(const Matrix<T>& data, T eps, int minPts) {
        std::size_t n = data.rows();
        std::size_t dim = data.cols();

        std::vector<int> labels(n, DBSCAN_NOISE);
        std::vector<bool> visited(n, false);
        int clusterID = 0;

        // ユークリッド距離
        auto distance = [&](std::size_t i, std::size_t j) -> T {
            T sum = T{0};
            for (std::size_t d = 0; d < dim; ++d) {
                T diff = data(i, d) - data(j, d);
                sum += diff * diff;
            }
            return std::sqrt(sum);
        };

        // 近傍探索
        auto regionQuery = [&](std::size_t p) -> std::vector<std::size_t> {
            std::vector<std::size_t> neighbors;
            for (std::size_t i = 0; i < n; ++i) {
                if (distance(p, i) <= eps) {
                    neighbors.push_back(i);
                }
            }
            return neighbors;
        };

        for (std::size_t p = 0; p < n; ++p) {
            if (visited[p]) continue;
            visited[p] = true;

            auto neighbors = regionQuery(p);

            if (static_cast<int>(neighbors.size()) < minPts) {
                labels[p] = DBSCAN_NOISE;
                continue;
            }

            // クラスタ展開
            labels[p] = clusterID;

            for (std::size_t k = 0; k < neighbors.size(); ++k) {
                std::size_t q = neighbors[k];

                if (!visited[q]) {
                    visited[q] = true;
                    auto qNeighbors = regionQuery(q);
                    if (static_cast<int>(qNeighbors.size()) >= minPts) {
                        neighbors.insert(neighbors.end(),
                                         qNeighbors.begin(), qNeighbors.end());
                    }
                }

                if (labels[q] == DBSCAN_NOISE) {
                    labels[q] = clusterID;
                }
            }

            ++clusterID;
        }

        return labels;
    }

} // namespace calx::ml

#endif // CALX_ML_CLUSTERING_DBSCAN_HPP
