C++代码实现MATLAB中的rocmetrics函数功能

// rocmetrics.cpp
// g++ -std=c++17 -O2 -o rocmetrics rocmetrics.cpp

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <iostream>
#include <limits>
#include <numeric>
#include <random>
#include <stdexcept>
#include <string>
#include <vector>

namespace rocmetrics {

// ————————————————————
// 基础数据结构
// ————————————————————
struct ROCCurve {
std::vector<double> thresholds;
std::vector<double> fpr;
std::vector<double> tpr;
double auc = 0.0;
};

struct ConfidenceInterval {
double lower = std::numeric_limits<double>::quiet_NaN();
double upper = std::numeric_limits<double>::quiet_NaN();
double level = 0.95;
};

struct ROCMetricsOptions {
bool compute_bootstrap_ci = false;
int num_bootstrap = 1000;
double ci_level = 0.95;
unsigned int seed = 0;
bool use_nearest_neighbor = false;
std::string nan_flag = "omitan";
};

// ————————————————————
// 内部工具
// ————————————————————
namespace detail {

inline void validate_labels(const std::vector<int>& labels) {
for (int y : labels) {
if (y != 0 && y != 1)
throw std::invalid_argument("labels must be 0 (negative) or 1 (positive)");
}
}

inline ROCCurve compute_roc_points_from_sorted(
const std::vector<int>& sorted_labels,
const std::vector<double>& sorted_scores)
{
const std::size_t n = sorted_labels.size();
if (n == 0) return {};

long P = 0, N = 0;
for (int y : sorted_labels) (y == 1 ? P : N)++;
if (P == 0 || N == 0)
throw std::invalid_argument("both positive and negative samples are required");

ROCCurve curve;
long tp = 0, fp = 0;

curve.thresholds.push_back(std::numeric_limits<double>::infinity());
curve.fpr.push_back(0.0);
curve.tpr.push_back(0.0);

std::size_t i = 0;
while (i < n) {
double score = sorted_scores[i];
while (i < n && sorted_scores[i] == score) {
if (sorted_labels[i] == 1) ++tp; else ++fp;
++i;
}
curve.thresholds.push_back(score);
curve.fpr.push_back(static_cast<double>(fp) / static_cast<double>(N));
curve.tpr.push_back(static_cast<double>(tp) / static_cast<double>(P));
}

if (curve.fpr.back() < 1.0 || curve.tpr.back() < 1.0) {
curve.thresholds.push_back(–std::numeric_limits<double>::infinity());
curve.fpr.push_back(1.0);
curve.tpr.push_back(1.0);
}
return curve;
}

inline double trapezoid_auc(const std::vector<double>& fpr,
const std::vector<double>& tpr) {
double area = 0.0;
for (std::size_t i = 1; i < fpr.size(); ++i) {
double dx = fpr[i] – fpr[i – 1];
area += 0.5 * dx * (tpr[i – 1] + tpr[i]);
}
return area;
}

inline double auc_mann_whitney(const std::vector<int>& labels,
const std::vector<double>& scores) {
std::size_t n = labels.size();
std::vector<std::size_t> idx(n);
std::iota(idx.begin(), idx.end(), 0);
std::sort(idx.begin(), idx.end(),
[&](std::size_t a, std::size_t b) { return scores[a] < scores[b]; });

long P = 0, N = 0;
for (int y : labels) (y == 1 ? P : N)++;
if (P == 0 || N == 0) return std::numeric_limits<double>::quiet_NaN();

double rank_sum_pos = 0.0;
std::size_t i = 0;
while (i < n) {
std::size_t j = i;
while (j < n && scores[idx[j]] == scores[idx[i]]) ++j;
double avg_rank = 0.5 * (i + 1 + j);
for (std::size_t k = i; k < j; ++k) {
if (labels[idx[k]] == 1) rank_sum_pos += avg_rank;
}
i = j;
}

double U = rank_sum_pos – 0.5 * P * (P + 1);
return U / (static_cast<double>(P) * static_cast<double>(N));
}

inline double quantile(const std::vector<double>& sorted, double q) {
if (sorted.empty()) return std::numeric_limits<double>::quiet_NaN();
if (q <= 0.0) return sorted.front();
if (q >= 1.0) return sorted.back();
double pos = q * (sorted.size() – 1);
std::size_t lo = static_cast<std::size_t>(std::floor(pos));
std::size_t hi = static_cast<std::size_t>(std::ceil(pos));
double w = pos – lo;
return sorted[lo] * (1.0 – w) + sorted[hi] * w;
}

} // namespace detail

// ————————————————————
// 二分类 ROC
// ————————————————————
inline ROCCurve compute_binary_roc(const std::vector<int>& labels,
const std::vector<double>& scores,
const ROCMetricsOptions& opt = {}) {
if (labels.size() != scores.size())
throw std::invalid_argument("labels and scores must have same length");
detail::validate_labels(labels);

std::vector<int> y;
std::vector<double> s;
y.reserve(labels.size());
s.reserve(scores.size());
for (std::size_t i = 0; i < labels.size(); ++i) {
if (opt.nan_flag == "omitan" && !std::isfinite(scores[i])) continue;
y.push_back(labels[i]);
s.push_back(scores[i]);
}
if (y.empty())
throw std::invalid_argument("no valid samples after NaN handling");

std::vector<std::size_t> idx(y.size());
std::iota(idx.begin(), idx.end(), 0);
std::sort(idx.begin(), idx.end(),
[&](std::size_t a, std::size_t b) { return s[a] > s[b]; });

std::vector<int> sorted_y(y.size());
std::vector<double> sorted_s(y.size());
for (std::size_t i = 0; i < idx.size(); ++i) {
sorted_y[i] = y[idx[i]];
sorted_s[i] = s[idx[i]];
}

ROCCurve curve = detail::compute_roc_points_from_sorted(sorted_y, sorted_s);
curve.auc = detail::trapezoid_auc(curve.fpr, curve.tpr);
return curve;
}

// ————————————————————
// Bootstrap AUC 置信区间
// ————————————————————
inline ConfidenceInterval bootstrap_auc_ci(
const std::vector<int>& labels,
const std::vector<double>& scores,
const ROCMetricsOptions& opt = {})
{
if (labels.size() != scores.size())
throw std::invalid_argument("labels and scores must have same length");
detail::validate_labels(labels);
if (labels.empty()) throw std::invalid_argument("empty input");

std::mt19937 rng(opt.seed == 0 ? std::random_device{}() : opt.seed);
std::uniform_int_distribution<std::size_t> dist(0, labels.size() – 1);

std::vector<double> boot_auc;
boot_auc.reserve(opt.num_bootstrap);

for (int b = 0; b < opt.num_bootstrap; ++b) {
std::vector<int> by;
std::vector<double> bs;
by.reserve(labels.size());
bs.reserve(scores.size());
for (std::size_t i = 0; i < labels.size(); ++i) {
std::size_t k = dist(rng);
by.push_back(labels[k]);
bs.push_back(scores[k]);
}
long P = 0, N = 0;
for (int y : by) (y == 1 ? P : N)++;
if (P == 0 || N == 0) continue;
double a = detail::auc_mann_whitney(by, bs);
if (std::isfinite(a)) boot_auc.push_back(a);
}

if (boot_auc.empty()) return {};
std::sort(boot_auc.begin(), boot_auc.end());

double alpha = 1.0 – opt.ci_level;
ConfidenceInterval ci;
ci.level = opt.ci_level;
ci.lower = detail::quantile(boot_auc, alpha / 2.0);
ci.upper = detail::quantile(boot_auc, 1.0 – alpha / 2.0);
return ci;
}

// ————————————————————
// 多分类 OVA
// ————————————————————
struct MulticlassROCMetrics {
std::vector<int> class_names;
std::vector<ROCCurve> roc_per_class;
std::vector<double> auc_per_class;
std::vector<ConfidenceInterval> ci_per_class;
double macro_auc = 0.0;
};

inline MulticlassROCMetrics compute_rocmetrics_multiclass(
const std::vector<int>& true_labels,
const std::vector<std::vector<double>>& scores,
const std::vector<int>& class_names,
const ROCMetricsOptions& opt = {})
{
const std::size_t N = true_labels.size();
const std::size_t K = class_names.size();
if (scores.size() != N)
throw std::invalid_argument("scores row count must equal label count");
for (const auto& row : scores)
if (row.size() != K)
throw std::invalid_argument("each score row must have K elements");

MulticlassROCMetrics result;
result.class_names = class_names;
result.roc_per_class.resize(K);
result.auc_per_class.resize(K);
result.ci_per_class.resize(K);

for (std::size_t k = 0; k < K; ++k) {
int pos_class = class_names[k];
std::vector<int> binary_y(N);
std::vector<double> binary_s(N);
for (std::size_t i = 0; i < N; ++i) {
binary_y[i] = (true_labels[i] == pos_class) ? 1 : 0;
binary_s[i] = scores[i][k];
}

long P = 0, Neg = 0;
for (int y : binary_y) (y == 1 ? P : Neg)++;
if (P == 0 || Neg == 0) {
result.auc_per_class[k] = std::numeric_limits<double>::quiet_NaN();
continue;
}

result.roc_per_class[k] = compute_binary_roc(binary_y, binary_s, opt);
result.auc_per_class[k] = result.roc_per_class[k].auc;

if (opt.compute_bootstrap_ci)
result.ci_per_class[k] = bootstrap_auc_ci(binary_y, binary_s, opt);
}

double sum = 0.0;
int count = 0;
for (double a : result.auc_per_class) {
if (std::isfinite(a)) { sum += a; ++count; }
}
result.macro_auc = count > 0 ? sum / count : 0.0;
return result;
}

} // namespace rocmetrics

// ————————————————————
// main
// ————————————————————
int main() {
// —- 二分类示例 —-
std::vector<int> labels = {0, 1, 0, 1, 1, 0, 1, 0};
std::vector<double> scores = {0.1, 0.8, 0.3, 0.7, 0.9, 0.2, 0.6, 0.4};

rocmetrics::ROCMetricsOptions opt;
opt.compute_bootstrap_ci = true;
opt.num_bootstrap = 2000;
opt.ci_level = 0.95;

auto curve = rocmetrics::compute_binary_roc(labels, scores, opt);

std::cout << "AUC = " << curve.auc << "\\n";
std::cout << "Threshold\\tFPR\\tTPR\\n";
for (std::size_t i = 0; i < curve.thresholds.size(); ++i) {
std::cout << curve.thresholds[i] << "\\t"
<< curve.fpr[i] << "\\t"
<< curve.tpr[i] << "\\n";
}

auto ci = rocmetrics::bootstrap_auc_ci(labels, scores, opt);
std::cout << "95% CI: [" << ci.lower << ", " << ci.upper << "]\\n";

// —- 多分类 OVA 示例 —-
std::vector<int> y_multi = {0, 1, 2, 0, 1, 2, 0, 1, 2, 0};
std::vector<std::vector<double>> s_multi = {
{0.7, 0.2, 0.1}, {0.1, 0.8, 0.1}, {0.1, 0.2, 0.7},
{0.6, 0.3, 0.1}, {0.2, 0.7, 0.1}, {0.1, 0.1, 0.8},
{0.8, 0.1, 0.1}, {0.1, 0.9, 0.0}, {0.0, 0.1, 0.9},
{0.5, 0.4, 0.1}
};
std::vector<int> classes = {0, 1, 2};

auto multi = rocmetrics::compute_rocmetrics_multiclass(
y_multi, s_multi, classes, opt);

std::cout << "Macro AUC = " << multi.macro_auc << "\\n";
for (std::size_t k = 0; k < classes.size(); ++k) {
std::cout << "Class " << classes[k]
<< " AUC = " << multi.auc_per_class[k] << "\\n";
}
return 0;
}