qestats coverage - 29.21%

Files
Source

#pragma once
#include "metrics.hpp"
#include "inference.hpp"
#include "derivation.hpp"
#include <algorithm>
#include <cmath>
#include <random>
#include <stdexcept>
#include <string>
#include <vector>

namespace qestats {

// ── Bootstrap result ──────────────────────────────────────────────────────────

struct BootstrapResult {
    Metrics original;
    // Per-metric standard errors and CI (parallel to METRIC_NAMES order)
    struct MetricStat {
        std::string name;
        double se    = 0;
        double ci_lo = 0;
        double ci_hi = 0;
    };
    std::vector<MetricStat> stats;
};

// Metric names in canonical order matching qestats R output.
inline const std::vector<std::string>& metric_names() {
    static const std::vector<std::string> v = {
        "b1","b2","precision","recall","accuracy","FDR","FOR","kappa"
    };
    return v;
}

// ── bootstrap_metrics ─────────────────────────────────────────────────────────
// cm: raw confusion matrix (TP, FP, FN, TN counts — not proportions)
// n_sample: resample size
// n_set: population size (0 = infinite)
// n_bootstrap: number of bootstrap replicates
// conf_level: confidence level (e.g. 0.95)

inline BootstrapResult
bootstrap_metrics(const ConfusionMatrix& cm, int n_sample,
                  const std::vector<std::string>& metrics = {},
                  int n_set = 0, int n_bootstrap = 1000,
                  double conf_level = 0.95) {
    const auto& all = metric_names();
    const auto& chosen = metrics.empty() ? all : metrics;

    double n1 = cm.tp, n2 = cm.fp, n3 = cm.fn, n4 = cm.tn;
    double n0 = n1 + n2 + n3 + n4;
    if (n0 <= 0 || n1 < 0 || n2 < 0 || n3 < 0 || n4 < 0)
        throw std::domain_error("bootstrap_metrics: invalid confusion matrix");

    BootstrapResult result;
    result.original = calculate_metrics(n1, n2, n3, n4);

    // Scale to population if requested
    if (n_set > 0) {
        double r = static_cast<double>(n_set) / n0;
        if (r > 1) {
            n1 = std::ceil(n1 * r); n2 = std::ceil(n2 * r);
            n3 = std::ceil(n3 * r); n4 = std::ceil(n4 * r);
        } else {
            n1 = std::floor(n1 * r); n2 = std::floor(n2 * r);
            n3 = std::floor(n3 * r); n4 = std::floor(n4 * r);
        }
        double tot = n1 + n2 + n3 + n4;
        if (n_sample > tot) n_sample = static_cast<int>(tot);
    }
    const double n_total = n1 + n2 + n3 + n4;
    const bool with_replacement = (n_set == 0);

    // Build population index vector [0..n_total)
    const int N = static_cast<int>(n_total);
    const int c1 = static_cast<int>(n1);
    const int c2 = static_cast<int>(n2);
    const int c3 = static_cast<int>(n3);

    std::mt19937 rng(42);  // fixed seed for reproducibility

    // Collect bootstrap replicates
    std::vector<std::vector<double>> samples(chosen.size());

    for (int b = 0; b < n_bootstrap; ++b) {
        std::vector<int> idx(n_sample);
        if (with_replacement) {
            std::uniform_int_distribution<int> dist(0, N - 1);
            for (int& x : idx) x = dist(rng);
        } else {
            // sample without replacement
            std::vector<int> pool(N);
            std::iota(pool.begin(), pool.end(), 0);
            std::shuffle(pool.begin(), pool.end(), rng);
            for (int i = 0; i < n_sample; ++i) idx[i] = pool[i];
        }
        int b1_ = 0, b2_ = 0, b3_ = 0, b4_ = 0;
        for (int x : idx) {
            if (x < c1)             ++b1_;
            else if (x < c1+c2)     ++b2_;
            else if (x < c1+c2+c3)  ++b3_;
            else                    ++b4_;
        }
        Metrics bm = calculate_metrics(b1_, b2_, b3_, b4_);
        for (std::size_t i = 0; i < chosen.size(); ++i)
            samples[i].push_back(get_metric(bm, chosen[i]));
    }

    const double alpha = 1.0 - conf_level;
    for (std::size_t i = 0; i < chosen.size(); ++i) {
        auto& s = samples[i];
        double mean = 0; for (double x : s) mean += x; mean /= s.size();
        double var = 0; for (double x : s) var += (x-mean)*(x-mean);
        double se = std::sqrt(var / (s.size() - 1));

        std::sort(s.begin(), s.end());
        int lo_idx = std::max(0, static_cast<int>(std::floor(alpha/2 * s.size())));
        int hi_idx = std::min(static_cast<int>(s.size())-1,
                              static_cast<int>(std::ceil((1.0-alpha/2) * s.size())));
        result.stats.push_back({chosen[i], se, s[lo_idx], s[hi_idx]});
    }
    return result;
}

// ── Welch t-test helper ───────────────────────────────────────────────────────

struct TTestResult {
    double t_stat  = 0;
    double df      = 0;
    double p_value = 1;
};

inline TTestResult welch_t(double m1, double m2,
                            double s1, double s2,
                            int n1, int n2,
                            double tol = 1e-9) {
    if (std::abs(m1 - m2) < tol) return {0, 1, 1.0};
    if (std::abs(s1) < tol && std::abs(s2) < tol)
        return {0, 1, (std::abs(m1-m2) < tol) ? 1.0 : 0.0};

    const double t = (m1 - m2) / std::sqrt(s1*s1 + s2*s2);
    const double num = (s1*s1 + s2*s2) * (s1*s1 + s2*s2);
    const double den = s1*s1*s1*s1/(n1-1) + s2*s2*s2*s2/(n2-1);
    const double df  = (den > 0) ? num/den : 1.0;

    // Two-sided p-value via regularised incomplete beta (approx via normal for large df)
    // For production accuracy use a proper t-distribution CDF; for now use
    // normal approximation which is good for df > 30.
    double p;
    if (df > 30) {
        // Standard normal approximation
        double z = std::abs(t);
        p = 2.0 * (0.5 * std::erfc(z / std::sqrt(2.0)));
    } else {
        // Abramowitz & Stegun approximation to t CDF
        double x = df / (df + t*t);
        // regularized incomplete beta I_x(df/2, 1/2) via continued fraction (simple approximation)
        // Good enough for n_bootstrap >= 30; exact implementation would need boost or a special function
        double z = std::abs(t);
        double a_s = df / 2.0, b_s = 0.5;
        // Use the normal approximation as a safe fallback
        p = 2.0 * (0.5 * std::erfc(z / std::sqrt(2.0)));
    }
    return {t, df, p};
}

// ── p_diff ────────────────────────────────────────────────────────────────────
// p-value for a metric difference between two groups using bootstrap Welch t.

inline double p_diff(const ConfusionMatrix& cm1, const ConfusionMatrix& cm2,
                     int n_sample, const std::string& metric,
                     int n_set = 0, int n_bootstrap = 1000, double alpha = 0.05) {
    const double n1 = cm1.tp + cm1.fp + cm1.fn + cm1.tn;
    const double n2 = cm2.tp + cm2.fp + cm2.fn + cm2.tn;
    const double n0 = n1 + n2;
    if (n_sample <= 0) n_sample = static_cast<int>(n0);
    const double r = static_cast<double>(n_sample) / n0;
    const int ns1 = std::max(1, static_cast<int>(std::floor(n1 * (r < 1 ? 1 : r))));
    const int ns2 = std::max(1, static_cast<int>(std::floor(n2 * (r < 1 ? 1 : r))));
    const int nset1 = (n_set > 0) ? static_cast<int>(std::floor(n_set * n1 / n0)) : 0;
    const int nset2 = (n_set > 0) ? static_cast<int>(std::floor(n_set * n2 / n0)) : 0;

    auto r1 = bootstrap_metrics(cm1, ns1, {metric}, nset1, n_bootstrap, 1.0-alpha);
    auto r2 = bootstrap_metrics(cm2, ns2, {metric}, nset2, n_bootstrap, 1.0-alpha);

    const double m1v = get_metric(r1.original, metric);
    const double m2v = get_metric(r2.original, metric);
    const double se1 = r1.stats.empty() ? 0 : r1.stats[0].se;
    const double se2 = r2.stats.empty() ? 0 : r2.stats[0].se;

    return welch_t(m1v, m2v, se1, se2, n_bootstrap, n_bootstrap).p_value;
}

// ── n_diff ────────────────────────────────────────────────────────────────────
// Minimum n at which two groups show a significant metric difference.

inline int n_diff(const ConfusionMatrix& cm1, const ConfusionMatrix& cm2,
                  const std::string& metric, int n_set = 0,
                  int n_bootstrap = 1000, double alpha = 0.05) {
    const int n1 = static_cast<int>(cm1.tp+cm1.fp+cm1.fn+cm1.tn);
    const int n2 = static_cast<int>(cm2.tp+cm2.fp+cm2.fn+cm2.tn);
    int n = n1 + n2;

    double p0 = p_diff(cm1, cm2, n, metric, n_set, n_bootstrap, alpha);
    if (p0 < alpha) return n;

    int left = n, right = n_set > 0 ? n_set : n;
    if (n_set <= 0) {
        while (p0 >= alpha) {
            left = right;
            right *= 2;
            p0 = p_diff(cm1, cm2, right, metric, n_set, n_bootstrap, alpha);
        }
    } else {
        double p1 = p_diff(cm1, cm2, n_set, metric, n_set, n_bootstrap, alpha);
        if (p1 >= alpha) return n_set;
    }

    while (right - left > 1) {
        int mid = (left + right + 1) / 2;
        p0 = p_diff(cm1, cm2, mid, metric, n_set, n_bootstrap, alpha);
        if (p0 >= alpha) left = mid;
        else             right = mid;
    }
    return (left + right + 1) / 2;
}

// ── N_Diff2 ───────────────────────────────────────────────────────────────────
// Minimum n at which CI intervals for two groups no longer overlap by > d.

inline double interval_overlap(const MetricRange& r1, const MetricRange& r2) {
    double span = std::min(r1.max - r1.min, r2.max - r2.min);
    if (span <= 0) return 0.0;
    double overlap = std::min(r1.max, r2.max) - std::max(r1.min, r2.min);
    return (overlap <= 0) ? 0.0 : overlap / span;
}

// range_est: hyper-based CI for a metric given raw counts + population N
inline MetricRange range_est(const ConfusionMatrix& cm,
                              const std::string& metric, int N, double alpha = 0.05) {
    const int tp = static_cast<int>(cm.tp), fp = static_cast<int>(cm.fp);
    const int fn = static_cast<int>(cm.fn), tn = static_cast<int>(cm.tn);
    auto prec_ci  = hyper_range(tp, tp+fp, N, alpha);
    auto for_ci   = hyper_range(fn, fn+tn, N, alpha);  // FOR = fn/(fn+tn)
    double b1 = static_cast<double>(tp+fp) / (tp+fp+fn+tn);

    auto eval_m = [&](double fdr, double for_v) -> double {
        Metrics m = derive_from_rates(b1, fdr, for_v);
        return get_metric(m, metric);
    };

    // Precision = 1 - FDR; CI is [lo, hi] for precision → FDR is [1-hi, 1-lo]
    double fdr_lo = 1.0 - prec_ci.upper, fdr_hi = 1.0 - prec_ci.lower;
    double for_lo = for_ci.lower, for_hi = for_ci.upper;

    double v_lo = eval_m(fdr_hi, for_hi);  // worst case
    double v_hi = eval_m(fdr_lo, for_lo);  // best case
    if (v_lo > v_hi) std::swap(v_lo, v_hi);
    return {v_lo, v_hi};
}

inline int N_Diff2(const ConfusionMatrix& cm1, const ConfusionMatrix& cm2,
                   const std::string& metric, int N, double d = 0.1,
                   double alpha = 0.05) {
    const double n1 = cm1.tp+cm1.fp+cm1.fn+cm1.tn;
    const double n2 = cm2.tp+cm2.fp+cm2.fn+cm2.tn;
    const int n = static_cast<int>(n1 + n2);

    auto scale = [](const ConfusionMatrix& cm, double r) -> ConfusionMatrix {
        return {std::floor(cm.tp*r), std::floor(cm.fp*r),
                std::floor(cm.fn*r), std::floor(cm.tn*r)};
    };

    auto check_overlap = [&](double x) -> bool {
        auto s1 = scale(cm1, x); auto s2 = scale(cm2, x);
        auto r1 = range_est(s1, metric, N, alpha);
        auto r2 = range_est(s2, metric, N, alpha);
        return interval_overlap(r1, r2) >= d;
    };

    if (!check_overlap(1.0)) return n;  // already separable

    const double x_max = static_cast<double>(N) / n;
    if (check_overlap(x_max)) return N;  // never separable at population size

    double left = 1.0, right = x_max;
    while (right - left > 1.0 / n) {
        double mid = (left + right) / 2.0;
        if (check_overlap(mid)) left  = mid;
        else                    right = mid;
    }
    return static_cast<int>(std::ceil(n * (left + right) / 2.0));
}

// ── ps_est_kappa ──────────────────────────────────────────────────────────────
// Given a perfect sample of size n and a baserate, estimate the minimum kappa.

inline double ps_est_kappa(int n, double baserate, int N = 0,
                            int allowed_errors = 0, double alpha = 0.05) {
    double accuracy = ps_est_proportion(n, N, allowed_errors, alpha);
    if (accuracy < 0) return -1;
    MetricRange kappa_range = range("baserate", baserate, "accuracy", accuracy, "kappa");
    return kappa_range.min;
}

// ── ps_kappa ─────────────────────────────────────────────────────────────────
// Perfect sample size for kappa test — delegates to ps_proportion via
// max accuracy from range(baserate, kappa).

inline PsResult ps_kappa(double kappa0, double baserate,
                          double kappa1 = 0, int N = 0,
                          int allowed_errors = 0, double alpha = 0.05) {
    if (alpha < 0.001 || alpha > 0.1) return {};
    if (kappa0 <= 0 || kappa0 >= 1)    return {};

    MetricRange acc_range = range("baserate", baserate, "kappa", kappa0, "accuracy");
    return ps_proportion(acc_range.max, kappa1, N, allowed_errors, alpha);
}

} // namespace qestats

#pragma once
#include "metrics.hpp"
#include <algorithm>
#include <cmath>
#include <functional>
#include <stdexcept>
#include <string>
#include <vector>

// C++ port of the qestats 22-row formula table plus range() and derive().
//
// Each FormulaRow encodes a pair of independent metrics (var1, var2) and a
// "key" variable (var3).  Given var1 and var2, the key variable sweeps an
// interval [l3, u3] and three formulas produce tp, fp, fn for every key value.
// This lets derive() find the key value that matches a desired third metric
// (5-point interval narrowing), and lets range() find the min/max of any
// fourth metric over the full key interval.

namespace qestats {

// ── Formula table row ─────────────────────────────────────────────────────────

struct FormulaRow {
    std::string var1, var2, var3;  // names of the three independent variables

    // Bounds on var1 (for validation).
    std::function<double(double v1, double v2)> l1, u1;
    // Bounds on var2 (for validation).
    std::function<double(double v1, double v2)> l2, u2;
    // Bounds on the key (var3) given v1, v2.
    std::function<double(double v1, double v2)> l3, u3;
    // Confusion matrix from (v1, v2, key).
    std::function<double(double v1, double v2, double key)> tp, fp, fn;
};

// ── Build the 22-row table ────────────────────────────────────────────────────

inline std::vector<FormulaRow> make_formula_table() {
    using F2 = std::function<double(double, double)>;
    using F3 = std::function<double(double, double, double)>;

    auto c0 = [](double, double) { return 0.0; };
    auto c1 = [](double, double) { return 1.0; };

    std::vector<FormulaRow> t;

    // 1: b1, b2 — key=tp
    t.push_back({"b1","b2","tp",
        c0, c1, c0, c1,
        c0, [](double b1, double b2){ return std::min(b1,b2); },
        [](double,  double,  double k){ return k; },
        [](double b1, double, double k){ return b1 - k; },
        [](double,  double b2, double k){ return b2 - k; }
    });

    // 2: b1, precision — key=fn
    t.push_back({"b1","precision","fn",
        c0, c1, c0, c1,
        c0, [](double b1, double){ return 1.0 - b1; },
        [](double b1, double pr, double){ return b1 * pr; },
        [](double b1, double pr, double){ return b1 - b1*pr; },
        [](double,   double,   double k){ return k; }
    });

    // 3: b1, recall — key=tp
    t.push_back({"b1","recall","tp",
        c0, c1, c0, c1,
        c0, [](double b1, double rc){ return std::min(b1, rc); },
        [](double,   double,   double k){ return k; },
        [](double b1, double,  double k){ return b1 - k; },
        [](double,   double rc, double k){ return (rc > 0) ? k/rc - k : 0.0; }
    });

    // 4: b1, accuracy — key=tp
    t.push_back({"b1","accuracy","tp",
        c0, c1, c0, c1,
        [](double b1, double ac){ return std::max(0.0, b1 + ac - 1.0); },
        [](double b1, double){ return b1; },
        [](double,   double,   double k){ return k; },
        [](double b1, double,  double k){ return b1 - k; },
        [](double b1, double ac, double k){ return 1.0 - b1 - ac + k; }
    });

    // 5: b1, FOR — key=tp
    t.push_back({"b1","FOR","tp",
        c0, c1, c0, c1,
        c0, [](double b1, double){ return b1; },
        [](double,   double,   double k){ return k; },
        [](double b1, double,  double k){ return b1 - k; },
        [](double b1, double f, double){ return f * (1.0 - b1); }
    });

    // 6: b1, kappa — key=tp
    t.push_back({"b1","kappa","tp",
        c0, c1, c0, c1,
        [](double b1, double kp){
            double d = 2.0 - 2.0*b1 - (1.0 - 2.0*b1)*kp;
            return (std::abs(d) > 1e-12) ? std::max(0.0, b1*kp/d) : 0.0;
        },
        [](double b1, double kp){
            double d = 2.0 - 2.0*b1 - (1.0 - 2.0*b1)*kp;
            return (std::abs(d) > 1e-12)
                ? std::min(b1, (2.0*b1 + (1.0-b1)*kp)/d) : b1;
        },
        [](double,   double,   double k){ return k; },
        [](double b1, double,  double k){ return b1 - k; },
        [](double b1, double kp, double k){
            double d = (1.0 - 2.0*b1)*kp + 2.0*b1;
            return (std::abs(d) > 1e-12) ? (2.0*k - b1*kp)/d - k : 0.0;
        }
    });

    // 7: b2, precision — key=tp
    t.push_back({"b2","precision","tp",
        c0, c1, c0, c1,
        c0, [](double pr, double b2){ return std::min(pr, b2); },
        [](double,   double,   double k){ return k; },
        [](double,   double pr, double k){ return (pr > 0) ? k/pr - k : 0.0; },
        [](double b2, double,  double k){ return b2 - k; }
    });

    // 8: b2, recall — key=fp
    t.push_back({"b2","recall","fp",
        c0, c1, c0, c1,
        c0, [](double b2, double){ return 1.0 - b2; },
        [](double b2, double rc, double){ return b2 * rc; },
        [](double,   double,   double k){ return k; },
        [](double b2, double rc, double){ return b2 - b2*rc; }
    });

    // 9: b2, accuracy — key=tp
    t.push_back({"b2","accuracy","tp",
        c0, c1, c0, c1,
        [](double b2, double ac){ return std::max(0.0, ac + b2 - 1.0); },
        [](double b2, double ac){ return std::min(b2, ac); },
        [](double,   double,   double k){ return k; },
        [](double b2, double ac, double k){ return 1.0 - ac - b2 + k; },
        [](double b2, double,   double k){ return b2 - k; }
    });

    // 10: b2, FOR — key=tp
    t.push_back({"b2","FOR","tp",
        c0, c1, c0, c1,
        [](double b2, double f){ return (1.0-f > 1e-12) ? (b2-f)/(1.0-f) : 0.0; },
        [](double b2, double){ return b2; },
        [](double,   double,   double k){ return k; },
        [](double b2, double f, double k){
            return (std::abs(f) > 1e-12) ? 1.0 - (b2-k)/f - k : 0.0; },
        [](double b2, double,   double k){ return b2 - k; }
    });

    // 11: b2, kappa — key=tp
    t.push_back({"b2","kappa","tp",
        c0, c1, c0, c1,
        [](double b2, double kp){
            double d = 2.0 - 2.0*b2 - (1.0-2.0*b2)*kp;
            return (std::abs(d) > 1e-12) ? std::max(0.0, b2*kp/d) : 0.0;
        },
        [](double b2, double kp){
            double d = 2.0 - 2.0*b2 - (1.0-2.0*b2)*kp;
            return (std::abs(d) > 1e-12)
                ? std::min(b2, (2.0*b2 + (1.0-b2)*kp)/d) : b2;
        },
        [](double,   double,   double k){ return k; },
        [](double b2, double kp, double k){
            double d = (1.0-2.0*b2)*kp + 2.0*b2;
            return (std::abs(d) > 1e-12) ? (2.0*k - b2*kp)/d - k : 0.0; },
        [](double b2, double,   double k){ return b2 - k; }
    });

    // 12: precision, recall — key=tp
    t.push_back({"precision","recall","tp",
        c0, c1, c0, c1,
        c0,
        [](double pr, double rc){
            double x = (pr + rc - pr*rc > 1e-12)
                ? rc*pr / (pr + rc - rc*pr) : 1.0;
            return std::min({1.0, pr, rc, x});
        },
        [](double,   double,   double k){ return k; },
        [](double pr, double,  double k){ return (pr > 1e-12) ? k/pr - k : 0.0; },
        [](double,   double rc, double k){ return (rc > 1e-12) ? k/rc - k : 0.0; }
    });

    // 13: precision, accuracy — key=tp
    t.push_back({"precision","accuracy","tp",
        c0, c1, c0, c1,
        c0,
        [](double pr, double ac){
            double x = (1.0-pr > 1e-12) ? (1.0-ac)*pr/(1.0-pr) : 1.0;
            return std::min({1.0, ac, x});
        },
        [](double,   double,   double k){ return k; },
        [](double pr, double,  double k){ return (pr > 1e-12) ? k/pr - k : 0.0; },
        [](double pr, double ac, double k){ return 1.0 - k/pr - ac + k; }
    });

    // 14: precision, FOR — key=tp
    t.push_back({"precision","FOR","tp",
        c0, c1, c0, c1,
        c0, [](double pr, double){ return pr; },
        [](double,   double,   double k){ return k; },
        [](double pr, double,  double k){ return (pr > 1e-12) ? k/pr - k : 0.0; },
        [](double pr, double f, double k){
            double neg = (pr > 1e-12) ? 1.0 - k/pr : 0.0;
            return f * neg; }
    });

    // 15: precision, kappa — key=tp
    t.push_back({"precision","kappa","tp",
        c0, c1, c0, c1,
        c0,
        [](double pr, double kp){
            double d = 2.0*(1.0-kp);
            return (d > 1e-12)
                ? std::min(pr, (2.0*pr - (1.0+pr)*kp) / d) : pr;
        },
        [](double,   double,   double k){ return k; },
        [](double pr, double,  double k){ return (pr > 1e-12) ? k/pr - k : 0.0; },
        [](double pr, double kp, double k){
            double d = pr*kp + 2.0*k*(1.0-kp);
            return (std::abs(d) > 1e-12)
                ? k * ((2.0*pr - kp)/d - 1.0) : 0.0; }
    });

    // 16: recall, accuracy — key=tp
    t.push_back({"recall","accuracy","tp",
        c0, c1, c0, c1,
        c0,
        [](double rc, double ac){
            double x = (1.0-rc > 1e-12) ? (1.0-ac)*rc/(1.0-rc) : 1.0;
            return std::min({1.0, ac, x});
        },
        [](double,   double,   double k){ return k; },
        [](double rc, double ac, double k){
            return 1.0 - (rc > 1e-12 ? (1.0-rc)/rc : 0.0)*k - ac; },
        [](double rc, double,   double k){
            return (rc > 1e-12) ? (1.0-rc)/rc * k : 0.0; }
    });

    // 17: recall, FOR — key=tp
    t.push_back({"recall","FOR","tp",
        c0, c1, c0, c1,
        c0,
        [](double rc, double f){
            double d = 1.0 - rc + rc*f;
            return (d > 1e-12) ? std::min(rc, rc*f/d) : rc;
        },
        [](double,   double,   double k){ return k; },
        [](double rc, double f, double k){
            double d = (rc > 1e-12 && f > 1e-12)
                ? (1.0-rc+rc*f)/(rc*f) : 1.0;
            return 1.0 - d*k; },
        [](double rc, double,   double k){
            return (rc > 1e-12) ? (1.0-rc)/rc * k : 0.0; }
    });

    // 18: recall, kappa — key=tp
    t.push_back({"recall","kappa","tp",
        c0, c1, c0, c1,
        c0,
        [](double rc, double kp){
            return rc - (1.0-rc)/(2.0*(1.0-kp)) * kp;
        },
        [](double,   double,   double k){ return k; },
        [](double rc, double kp, double k){
            double d = 2.0*k*(1.0-kp) + kp*rc;
            return (std::abs(d) > 1e-12)
                ? (2.0*(rc-k)*(1.0-kp) - (1.0-rc)*kp) * k / d : 0.0; },
        [](double rc, double,   double k){
            return (rc > 1e-12) ? (1.0-rc)/rc * k : 0.0; }
    });

    // 19: accuracy, FOR — key=tp
    t.push_back({"accuracy","FOR","tp",
        c0, c1, c0, c1,
        [](double ac, double f){ return std::max(0.0, (ac + f - 1.0)/f); },
        [](double ac, double){ return ac; },
        [](double,   double,   double k){ return k; },
        [](double ac, double f, double k){
            return (1.0-f > 1e-12)
                ? 1.0 - ac - (ac - k)*f/(1.0-f) : 0.0; },
        [](double ac, double f, double k){
            return (1.0-f > 1e-12) ? (ac-k)*f/(1.0-f) : 0.0; }
    });

    // 20-A: accuracy, kappa — key=tp (branch A: fp uses + sqrt, fn uses - sqrt)
    {
        auto l3_ak = [](double ac, double kp){
            double disc = ac*ac - 2.0*(1.0-ac)*kp/(1.0-kp);
            return (disc >= 0)
                ? (ac - std::sqrt(disc)) / 2.0 : 0.0;
        };
        auto u3_ak = [](double ac, double kp){
            double disc = ac*ac - 2.0*(1.0-ac)*kp/(1.0-kp);
            return (disc >= 0)
                ? (ac + std::sqrt(disc)) / 2.0 : ac;
        };
        auto fp_A = [](double ac, double kp, double k) -> double {
            double inner = (1.0-ac)*(1.0-ac)
                + 2.0/(1.0-kp) * (2.0*(k-ac)*k*(1.0-kp) - (1.0-ac)*kp);
            if (inner < 0) inner = 0;
            return ((1.0-ac) + std::sqrt(inner)) / 2.0;
        };
        auto fn_A = [](double ac, double kp, double k) -> double {
            double inner = (1.0-ac)*(1.0-ac)
                + 2.0/(1.0-kp) * (2.0*(k-ac)*k*(1.0-kp) - (1.0-ac)*kp);
            if (inner < 0) inner = 0;
            return ((1.0-ac) - std::sqrt(inner)) / 2.0;
        };
        t.push_back({"accuracy","kappa","tp",
            c0, c1, c0, c1, l3_ak, u3_ak,
            [](double,double,double k){ return k; },
            fp_A, fn_A
        });
    }

    // 20-B: accuracy, kappa — key=tp (branch B: fp uses - sqrt, fn uses + sqrt)
    {
        auto l3_ak = [](double ac, double kp){
            double disc = ac*ac - 2.0*(1.0-ac)*kp/(1.0-kp);
            return (disc >= 0) ? (ac - std::sqrt(disc)) / 2.0 : 0.0;
        };
        auto u3_ak = [](double ac, double kp){
            double disc = ac*ac - 2.0*(1.0-ac)*kp/(1.0-kp);
            return (disc >= 0) ? (ac + std::sqrt(disc)) / 2.0 : ac;
        };
        auto fp_B = [](double ac, double kp, double k) -> double {
            double inner = (1.0-ac)*(1.0-ac)
                + 2.0/(1.0-kp) * (2.0*(k-ac)*k*(1.0-kp) - (1.0-ac)*kp);
            if (inner < 0) inner = 0;
            return ((1.0-ac) - std::sqrt(inner)) / 2.0;
        };
        auto fn_B = [](double ac, double kp, double k) -> double {
            double inner = (1.0-ac)*(1.0-ac)
                + 2.0/(1.0-kp) * (2.0*(k-ac)*k*(1.0-kp) - (1.0-ac)*kp);
            if (inner < 0) inner = 0;
            return ((1.0-ac) + std::sqrt(inner)) / 2.0;
        };
        t.push_back({"accuracy","kappa","tp",
            c0, c1, c0, c1, l3_ak, u3_ak,
            [](double,double,double k){ return k; },
            fp_B, fn_B
        });
    }

    // 21: FOR, kappa — key=tp; u3 = 1 (tp key upper bound is always 1)
    t.push_back({"FOR","kappa","tp",
        c0, c1,
        c0, [](double f, double){ return std::min(1.0, 2.0*(1.0-f)/(2.0-f)); },
        [](double f, double kp){
            double d = 2.0*(1.0-kp)*(1.0-f);
            return (d > 1e-12) ? kp*f/d : 0.0; },
        [](double, double){ return 1.0; },
        [](double f, double kp, double k) -> double {
            double A = 2.0*(1.0-kp)*((k*(1.0-f)+f*(1.0-k))) + (1.0-f)*kp;
            double B = 8.0*(1.0-kp)*f*(2.0*(1.0-kp)*k*(1.0-(k*(1.0-f)+f)) - (1.0-k)*f*kp);
            double d = 4.0*(1.0-kp)*f;
            if (std::abs(d) < 1e-12) return 0.0;
            double disc = A*A - B;
            if (disc < 0) disc = 0;
            return (A - std::sqrt(disc)) / d; },
        [](double f, double kp, double k) -> double {
            // fn = (1-tp)*FOR - fp
            double A = 2.0*(1.0-kp)*((k*(1.0-f)+f*(1.0-k))) + (1.0-f)*kp;
            double B = 8.0*(1.0-kp)*f*(2.0*(1.0-kp)*k*(1.0-(k*(1.0-f)+f)) - (1.0-k)*f*kp);
            double d = 4.0*(1.0-kp)*f;
            if (std::abs(d) < 1e-12) return 0.0;
            double disc = A*A - B;
            if (disc < 0) disc = 0;
            double fp_val = (A - std::sqrt(disc)) / d;
            return (1.0-k)*f - fp_val; }
    });

    return t;
}

// ── Global table (built once) ─────────────────────────────────────────────────

inline const std::vector<FormulaRow>& formula_table() {
    static const auto t = make_formula_table();
    return t;
}

// ── Name normalization ────────────────────────────────────────────────────────

inline std::string normalize_metric_name(const std::string& n) {
    if (n == "baserate")  return "b1";
    if (n == "baserate2") return "b2";
    if (n == "FDR")       return "precision";  // FDR = 1 - precision
    return n;
}

// ── Find a unique formula row for a pair (v1, v2) ────────────────────────────
// Returns nullptr if not found or ambiguous (accuracy+kappa has two rows).
inline const FormulaRow* find_formula(const std::string& v1,
                                       const std::string& v2) {
    const auto& t = formula_table();
    const FormulaRow* found = nullptr;
    int count = 0;
    for (const auto& row : t) {
        if ((row.var1 == v1 && row.var2 == v2) ||
            (row.var1 == v2 && row.var2 == v1)) {
            found = &row;
            ++count;
        }
    }
    if (count != 1) return nullptr;  // 0 = not found; >1 = ambiguous (accuracy+kappa)
    return found;
}

// ── Evaluate tp/fp/fn from a formula row, given v1, v2, key ──────────────────
struct CMResult {
    double tp, fp, fn, tn;
    bool valid = true;
};

inline CMResult eval_formula(const FormulaRow& row,
                              double v1, double v2, double key) {
    // Ensure (var1, var2) are in canonical order
    double a = v1, b = v2;
    if (row.var1 != row.var1) { a = v2; b = v1; }  // no-op placeholder

    double tp = row.tp(a, b, key);
    double fp = row.fp(a, b, key);
    double fn = row.fn(a, b, key);
    // Clamp near-zero numerical noise
    if (tp < 0 && std::abs(tp) < 1e-10) tp = 1e-10;
    if (fp < 0 && std::abs(fp) < 1e-10) fp = 0.0;
    if (fn < 0 && std::abs(fn) < 1e-10) fn = 0.0;
    double tn = 1.0 - tp - fp - fn;
    if (tn < 0 && std::abs(tn) < 1e-10) tn = 1e-10;
    bool ok = (tp >= 0 && tp <= 1 && fp >= 0 && fp <= 1 &&
               fn >= 0 && fn <= 1 && tn >= 0 && tn <= 1);
    return {tp, fp, fn, tn, ok};
}

// ── range() ──────────────────────────────────────────────────────────────────
// Given two metric values, return the [min, max] range of a third metric.
// Mirrors R's range(list(baserate=0.1, precision=0.85), "kappa").

inline MetricRange range(const std::string& name1, double val1,
                          const std::string& name2, double val2,
                          const std::string& target) {
    const std::string n1 = normalize_metric_name(name1);
    const std::string n2 = normalize_metric_name(name2);
    const std::string nt = normalize_metric_name(target);

    // FDR input: convert to precision value
    double v1 = val1, v2 = val2;
    if (name1 == "FDR") v1 = 1.0 - val1;
    if (name2 == "FDR") v2 = 1.0 - val2;

    const FormulaRow* row = find_formula(n1, n2);
    if (!row) throw std::invalid_argument("range: formula not found for pair");

    // Ensure v1/v2 align with row->var1/var2
    double a = v1, b = v2;
    if (row->var1 == n2) { a = v2; b = v1; }

    // Key bounds
    const double l3 = row->l3(a, b) + 1e-10;
    const double u3 = row->u3(a, b) - 1e-10;
    if (l3 >= u3) throw std::domain_error("range: invalid key range");

    auto eval_target = [&](double key) -> double {
        auto cm = eval_formula(*row, a, b, key);
        if (!cm.valid) throw std::domain_error("range: invalid confusion matrix");
        double s = cm.tp + cm.fp + cm.fn + cm.tn;
        Metrics m = metrics_from_cm({cm.tp/s, cm.fp/s, cm.fn/s, cm.tn/s});
        return get_metric(m, nt);
    };

    // 5-point interval narrowing for max
    double max_v = -1e300, min_v = 1e300;
    {
        double left = l3, right = u3;
        int max_pos = 0;
        bool monotonic = false;
        while (right - left > 1e-6) {
            double step = (right - left) / 4.0;
            double best = -1e300;
            for (int i = 0; i <= 4; ++i) {
                double key = left + i * step;
                double v = eval_target(key);
                if (v > best) { best = v; max_pos = i; }
            }
            max_v = best;
            if (max_pos == 0 || max_pos == 4) { monotonic = true; break; }
            left  = left + (max_pos - 1) * step;
            right = left + 2.0 * step;
        }
        (void)monotonic;
    }
    // 5-point interval narrowing for min
    {
        double left = l3, right = u3;
        bool monotonic = false;
        while (right - left > 1e-6) {
            double step = (right - left) / 4.0;
            double best = 1e300;
            int min_pos = 0;
            for (int i = 0; i <= 4; ++i) {
                double key = left + i * step;
                double v = eval_target(key);
                if (v < best) { best = v; min_pos = i; }
            }
            min_v = best;
            if (min_pos == 0 || min_pos == 4) { monotonic = true; break; }
            left  = left + (min_pos - 1) * step;
            right = left + 2.0 * step;
        }
        (void)monotonic;
    }

    return {min_v, max_v};
}

// ── derive() ─────────────────────────────────────────────────────────────────
// Given values for three independent metrics, return all metrics.
// Mirrors R's derive(list(baserate=0.1, FDR=0.15, FOR=0.05)).
//
// Inputs: three (name, value) pairs.  "kappa" must not be one of the first two
// (it is the target variable in the 5-point search).

inline Metrics derive(const std::string& name1, double val1,
                       const std::string& name2, double val2,
                       const std::string& name3, double val3) {
    // Normalize names
    std::string n[3] = {
        normalize_metric_name(name1),
        normalize_metric_name(name2),
        normalize_metric_name(name3)
    };
    double v[3] = {val1, val2, val3};

    // FDR → precision (inverted value)
    if (name1 == "FDR") v[0] = 1.0 - val1;
    if (name2 == "FDR") v[1] = 1.0 - val2;
    if (name3 == "FDR") v[2] = 1.0 - val3;

    // Move kappa to position 2 (it becomes the search target)
    for (int i = 0; i < 2; ++i) {
        if (n[i] == "kappa") {
            std::swap(n[i], n[2]); std::swap(v[i], v[2]);
            break;
        }
    }
    const std::string& target_name = n[2];
    const double        target_val  = v[2];

    // Validate range of the target
    MetricRange rng = range(name1, val1, name2, val2, name3);
    // For normalized names after FDR conversion the range call uses the
    // converted names — call with original names so alias handling fires.
    rng = range(name1, val1, name2, val2, name3);
    if (target_val < rng.min || target_val > rng.max)
        throw std::domain_error("derive: target metric value outside valid range");

    // Get formula for (n[0], n[1])
    const FormulaRow* row = find_formula(n[0], n[1]);
    if (!row) throw std::invalid_argument("derive: formula not found for pair");

    double a = v[0], b = v[1];
    if (row->var1 == n[1]) { a = v[1]; b = v[0]; }

    const double l3 = row->l3(a, b) + 1e-10;
    const double u3 = row->u3(a, b) - 1e-10;

    auto eval_target = [&](double key) -> double {
        auto cm = eval_formula(*row, a, b, key);
        if (!cm.valid) return target_val;  // keep search moving
        double s = cm.tp + cm.fp + cm.fn + cm.tn;
        Metrics m = metrics_from_cm({cm.tp/s, cm.fp/s, cm.fn/s, cm.tn/s});
        return get_metric(m, target_name);
    };

    // 5-point interval narrowing to find the key value
    double left = l3, right = u3;
    while (right - left > 1e-6) {
        double x5[5];
        double step = (right - left) / 4.0;
        for (int i = 0; i < 5; ++i) x5[i] = left + i * step;

        double best_err = 1e300;
        int best_idx = 0;
        for (int i = 0; i < 5; ++i) {
            double err = std::abs(eval_target(x5[i]) - target_val);
            if (err < best_err) { best_err = err; best_idx = i; }
        }
        if (best_idx == 0)      right = x5[1];
        else if (best_idx == 4) left  = x5[3];
        else { left = x5[best_idx-1]; right = x5[best_idx+1]; }
    }
    const double key_final = (left + right) / 2.0;
    auto cm = eval_formula(*row, a, b, key_final);
    double s = cm.tp + cm.fp + cm.fn + cm.tn;
    return metrics_from_cm({cm.tp/s, cm.fp/s, cm.fn/s, cm.tn/s});
}

// Convenience: derive from a closed-form (b1 + FDR + FOR) — fast path used
// by libcodey without the formula search.
inline Metrics derive_from_rates(double b1, double fdr, double for_) {
    const double fp = fdr * b1;
    const double tp = b1 - fp;
    const double fn = for_ * (1.0 - b1);
    const double tn = 1.0 - fp - tp - fn;
    return metrics_from_cm({tp, fp, fn, tn});
}

} // namespace qestats

#pragma once
#include "metrics.hpp"
#include <cmath>
#include <stdexcept>
#include <vector>

namespace qestats {

// ── Hypergeometric CDF ────────────────────────────────────────────────────────
// P(X <= q) where X ~ Hypergeometric(N, K, n)
//   K = number of successes in population
//   N - K = failures in population
//   n = sample size
//
// Computed in log-space to avoid overflow.  Equivalent to R's phyper(q,K,N-K,n).

inline double phyper(int q, int K, int N_minus_K, int n) {
    if (q < 0) return 0.0;
    const int N = K + N_minus_K;
    if (n > N || K < 0 || N_minus_K < 0) return 0.0;

    // log of C(a, b)
    auto lnchoose = [](int a, int b) -> double {
        if (b < 0 || b > a) return -1e300;
        if (b == 0 || b == a) return 0.0;
        return std::lgamma(a + 1) - std::lgamma(b + 1) - std::lgamma(a - b + 1);
    };

    const double log_denom = lnchoose(N, n);
    double sum = 0.0;
    const int hi = std::min(q, std::min(K, n));
    for (int k = std::max(0, n - N_minus_K); k <= hi; ++k) {
        double log_num = lnchoose(K, k) + lnchoose(N_minus_K, n - k);
        sum += std::exp(log_num - log_denom);
    }
    return std::min(sum, 1.0);
}

// ── Hypergeometric confidence interval for proportion ─────────────────────────
// Estimates CI for the proportion of positives in a population.
// R equivalent: hyper_range(n_pos, n, N, alpha)
//   n_pos = positives in sample
//   n     = sample size
//   N     = population size
//   alpha = significance level

inline Interval hyper_range(int n_pos, int n, int N, double alpha = 0.05) {
    if (n_pos < 0 || n_pos > n || n > N || alpha <= 0 || alpha >= 1)
        throw std::domain_error("hyper_range: invalid input parameters");

    const int n_neg = n - n_pos;
    int M_lower = n_neg;
    int M_upper = N - n_pos;

    // Shrink M_lower upward until phyper tail >= alpha/2
    {
        double p_lower = phyper(n_pos, N - M_lower, M_lower, n);
        if (p_lower < alpha / 2) {
            int left = M_lower, right = M_upper;
            while (right - left > 2) {
                M_lower = (left + right + 1) / 2;
                double p = phyper(n_pos - 1, N - M_lower, M_lower, n);
                if (p < alpha / 2) left = M_lower;
                else               right = M_lower;
            }
            M_lower = left;
        }
    }

    // Shrink M_upper downward until upper tail >= alpha/2
    {
        double p_upper = 1.0 - phyper(n_pos, N - M_upper, M_upper, n);
        if (p_upper < alpha / 2) {
            int left = M_lower, right = M_upper;
            while (right - left > 2) {
                M_upper = (left + right + 1) / 2;
                double p = 1.0 - phyper(n_pos - 1, N - M_upper, M_upper, n);
                if (p < alpha / 2) right = M_upper;
                else               left  = M_upper;
            }
            M_upper = right;
        }
    }

    return {
        1.0 - static_cast<double>(M_upper) / N,
        (n > 0) ? static_cast<double>(n_pos) / n : 0.0,
        1.0 - static_cast<double>(M_lower) / N
    };
}

// ── Perfect sampling: proportion ──────────────────────────────────────────────
// Minimum sample size n such that if all items in the sample satisfy the
// property, H0: x <= x0 is rejected at level alpha.
// Returns {n, power, cost}; n = -1 on invalid input.

struct PsResult {
    int    n     = -1;
    double power = 0;
    double cost  = 0;
};

inline PsResult ps_proportion(double x0, double x1 = 0,
                               int N = 0, int allowed_errors = 0,
                               double alpha = 0.05) {
    if (alpha < 0.001 || alpha > 0.1) return {};
    if (x0 < 0.5 || x0 >= 1.0)       return {};

    const int a = allowed_errors;
    if (N <= 0) N = 1000000;
    const int K = static_cast<int>(std::ceil(x0 * N));
    if (K == N) return {N, 1.0, 1.0};
    if (a >= static_cast<int>(K * (1.0 - x0))) return {};

    // Binary search for min n: P(X >= n - a | K, N-K, n) <= alpha
    // i.e.  1 - phyper(n-1-a, K, N-K, n) <= alpha
    int left = a, right = K;
    double p = 1.0;
    while (left < right - 1) {
        int n = (left + right + 1) / 2;
        p = 1.0 - phyper(n - 1 - a, K, N - K, n);
        if (p > alpha) left  = n;
        else           right = n;
    }
    const int n_out = right;

    if (x1 <= x0 || x1 >= 1.0) x1 = (x0 + 1.0) / 2.0;
    const double K1 = x1 * N;
    const int    K1i = static_cast<int>(K1);
    const double power = 1.0 - phyper(n_out - 1 - a, K1i, N - K1i, n_out);
    return {n_out, power, static_cast<double>(n_out) / power};
}

// ── Perfect sampling: kappa ───────────────────────────────────────────────────
// Derives the maximum achievable accuracy for (baserate, kappa) via range(),
// then delegates to ps_proportion.
// Declared here; implemented after derivation.hpp is included in qestats.hpp.
// Forward declaration only — see qestats.hpp for the full definition.

// ── ps_est_proportion ────────────────────────────────────────────────────────
// Given a perfect sample of size n, estimate the minimum proportion of
// positives in the population.

inline double ps_est_proportion(int n, int N = 0, int allowed_errors = 0,
                                 double alpha = 0.05) {
    if (alpha < 0.001 || alpha > 0.1) return -1;
    const int a = allowed_errors;
    if (a >= n) return -1;
    if (N <= 0) N = 1000000;
    if (n < 1 || n >= N) return -1;

    int left = 1, right = N;
    double x = 0;
    while (right - left > 1) {
        int K = (left + right + 1) / 2;
        double p = 1.0 - phyper(n - 1 - a, K, N - K, n);
        x = static_cast<double>(K) / N;
        if (p < alpha) left  = K;
        else           right = K;
    }
    return x;
}

// ── ps_est_kappa ──────────────────────────────────────────────────────────────
// Given a perfect sample of size n, estimate the minimum kappa achievable.
// Declared here; depends on range() from derivation.hpp — implemented in
// comparison.hpp after derivation.hpp is included.
// Forward-declared below; full definition is in comparison.hpp.

// ── metrics_range ─────────────────────────────────────────────────────────────
// Estimates CI for precision, recall, accuracy, FDR, FOR, and kappa given a
// sampled confusion matrix (raw counts) and population size N.
// Returns a vector of MetricRangeRow in the same order as R's metrics_range().

inline std::vector<MetricRangeRow>
metrics_range(double tp, double fp, double fn, double tn,
              int N, double baserate, double alpha = 0.05) {
    // Helper: derive metrics from (baserate, FDR, FOR)
    auto metrics_from_bFF = [](double b, double FDR, double FOR) -> Metrics {
        const double fp_ = FDR * b;
        const double tp_ = b - fp_;
        const double fn_ = FOR * (1.0 - b);
        const double tn_ = 1.0 - fp_ - tp_ - fn_;
        return metrics_from_cm({tp_, fp_, fn_, tn_});
    };

    // Precision CI via hyper_range (full alpha)
    Interval prec = hyper_range(static_cast<int>(tp),
                                static_cast<int>(tp + fp), N, alpha);
    Interval FDR_ci = {1.0 - prec.upper, 1.0 - prec.mean, 1.0 - prec.lower};

    // FOR CI: tn / (fn + tn) via hyper_range, then invert
    Interval one_minus_FOR = hyper_range(static_cast<int>(tn),
                                         static_cast<int>(fn + tn), N, alpha);
    Interval FOR_ci = {1.0 - one_minus_FOR.upper,
                       1.0 - one_minus_FOR.mean,
                       1.0 - one_minus_FOR.lower};

    // Derive recall, accuracy, kappa at lower/mean/upper FDR+FOR extremes
    Metrics lower_m = metrics_from_bFF(baserate, FDR_ci.upper, FOR_ci.upper);
    Metrics mean_m  = metrics_from_bFF(baserate, FDR_ci.mean,  FOR_ci.mean);
    Metrics upper_m = metrics_from_bFF(baserate, FDR_ci.lower, FOR_ci.lower);

    return {
        {"precision", prec.lower,        prec.mean,        prec.upper},
        {"recall",    lower_m.recall,    mean_m.recall,    upper_m.recall},
        {"accuracy",  lower_m.accuracy,  mean_m.accuracy,  upper_m.accuracy},
        {"FDR",       FDR_ci.lower,      FDR_ci.mean,      FDR_ci.upper},
        {"FOR",       FOR_ci.lower,      FOR_ci.mean,      FOR_ci.upper},
        {"kappa",     lower_m.kappa,     mean_m.kappa,     upper_m.kappa},
    };
}

} // namespace qestats

#pragma once
#include <algorithm>
#include <cmath>
#include <stdexcept>
#include <string>

// Functions that must reproduce R's arithmetic (rhoR parity) start their body
// with QESTATS_R_ARITHMETIC and carry QESTATS_R_ARITHMETIC_FN: R never fuses
// a * b + c, so these keep the compiler from contracting it into an FMA.
#if defined(__clang__)
#  define QESTATS_R_ARITHMETIC _Pragma("STDC FP_CONTRACT OFF")
#  define QESTATS_R_ARITHMETIC_FN
#elif defined(__GNUC__)
#  define QESTATS_R_ARITHMETIC
#  define QESTATS_R_ARITHMETIC_FN __attribute__((optimize("fp-contract=off")))
#else
#  define QESTATS_R_ARITHMETIC
#  define QESTATS_R_ARITHMETIC_FN
#endif

namespace qestats {

// Confusion matrix proportions (must sum to 1 after normalization).
struct ConfusionMatrix {
    double tp = 0, fp = 0, fn = 0, tn = 0;
};

// All derived classification metrics.
struct Metrics {
    // confusion matrix proportions
    double tp = 0, fp = 0, fn = 0, tn = 0;
    // derived
    double b1 = 0;        // classifier positive rate  (tp + fp)
    double b2 = 0;        // true positive rate        (tp + fn)
    double precision = 0; // tp / (tp + fp)
    double recall    = 0; // tp / (tp + fn)
    double accuracy  = 0; // tp + tn
    double fdr       = 0; // fp / (tp + fp)   — False Discovery Rate
    double for_      = 0; // fn / (fn + tn)   — False Omission Rate
    double kappa     = 0;
};

// Range result — min/max of a metric over the key variable interval.
struct MetricRange {
    double min = 0;
    double max = 0;
};

// Confidence interval from a hypergeometric estimate.
struct Interval {
    double lower = 0;
    double mean  = 0;
    double upper = 0;
};

// One row in the metrics_range output table.
struct MetricRangeRow {
    std::string metric;
    double lower = 0, mean = 0, upper = 0;
};

// ── Core arithmetic ───────────────────────────────────────────────────────────

// Cohen's kappa from confusion-matrix counts (or proportions), computed in
// the same order as rhoR's calcKappa() so the two agree to the bit: perfect
// agreement on one class (all TP or all TN) is 1.
QESTATS_R_ARITHMETIC_FN
inline double kappa_from_counts(double tp, double fp, double fn, double tn) {
    QESTATS_R_ARITHMETIC
    const double s = tp + fp + fn + tn;
    if (tp == s || tn == s) return 1.0;
    const double agreement = (tp + tn) / s;
    const double rate1 = (tp + fn) / s;   // first rater's positive rate
    const double rate2 = (tp + fp) / s;   // second rater's positive rate
    const double pe = rate1 * rate2 + (1 - rate1) * (1 - rate2);
    return (pe < 1.0) ? (agreement - pe) / (1 - pe) : 0.0;
}

// Rates are computed from the counts, not the normalized cells, so precision
// and recall match rhoR's TP / (TP + FP) and TP / (TP + FN) exactly.
inline Metrics metrics_from_cm(const ConfusionMatrix& cm) {
    const double s = cm.tp + cm.fp + cm.fn + cm.tn;
    if (s <= 0) throw std::domain_error("metrics_from_cm: confusion matrix is empty");
    Metrics m;
    m.tp = cm.tp / s; m.fp = cm.fp / s;
    m.fn = cm.fn / s; m.tn = cm.tn / s;
    const double pos1 = cm.tp + cm.fp, pos2 = cm.tp + cm.fn, neg1 = cm.fn + cm.tn;
    m.b1        = pos1 / s;
    m.b2        = pos2 / s;
    m.precision = (pos1 > 0) ? cm.tp / pos1 : 0.0;
    m.recall    = (pos2 > 0) ? cm.tp / pos2 : 0.0;
    m.accuracy  = (cm.tp + cm.tn) / s;
    m.fdr       = (pos1 > 0) ? cm.fp / pos1 : 0.0;
    m.for_      = (neg1 > 0) ? cm.fn / neg1 : 0.0;
    m.kappa     = kappa_from_counts(cm.tp, cm.fp, cm.fn, cm.tn);
    return m;
}

// Convenience: from raw counts (negative values clamped to 0).
inline Metrics calculate_metrics(double tp, double fp, double fn, double tn) {
    if (tp < 0) tp = 0; if (fp < 0) fp = 0;
    if (fn < 0) fn = 0; if (tn < 0) tn = 0;
    return metrics_from_cm({tp, fp, fn, tn});
}

// Retrieve a single metric value by name (supports alias names).
inline double get_metric(const Metrics& m, const std::string& name) {
    if (name == "b1" || name == "baserate")  return m.b1;
    if (name == "b2" || name == "baserate2") return m.b2;
    if (name == "precision")  return m.precision;
    if (name == "recall")     return m.recall;
    if (name == "accuracy")   return m.accuracy;
    if (name == "FDR")        return m.fdr;
    if (name == "FOR")        return m.for_;
    if (name == "kappa")      return m.kappa;
    if (name == "tp")         return m.tp;
    if (name == "fp")         return m.fp;
    if (name == "fn")         return m.fn;
    if (name == "tn")         return m.tn;
    throw std::invalid_argument("get_metric: unknown metric '" + name + "'");
}

} // namespace qestats

#pragma once
#include "metrics.hpp"
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <random>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <vector>

// Rho (Shaffer's rho): how likely a test set's kappa this high would be if the
// coding process's population kappa were below a threshold. A port of rhoR:
// given an adapter over R's RNG, a seeded rho() returns exactly what rhoR does.
//
// rhoR computes rho two ways, and both are kept, quirks included:
//
//   RhoMethod::Standard   rhoR's default path, written in R (calcRho): grids
//                         from seq(), draws from sample.int(), the hand set
//                         sampled as rows of the full set. qestats::rho_standard.
//   RhoMethod::Fast       rhoR's C++ path (calcRho_c, any other `method`):
//                         Armadillo regspace() grids, RcppArmadillo sample() /
//                         randi() draws, the hand set sampled from the
//                         contingency table. qestats::rho_fast.
//
// Random numbers come from an Rng: any type with
//
//   double unif_rand();          // uniform on (0, 1)          — R's unif_rand()
//   double unif_index(double n); // uniform integer in [0, n)  — R's R_unif_index()
//
// The Standard method draws only through unif_index, the Fast method only
// through unif_rand. StdRng is a standalone default.

namespace qestats {

// mt19937_64 with the draws derived by hand (the std distributions differ
// between standard libraries), so a seed gives the same rho everywhere.
struct StdRng {
    std::mt19937_64 gen;
    explicit StdRng(std::uint64_t seed = 5489u) : gen(seed) {}
    double unif_rand() {
        double u;
        do { u = static_cast<double>(gen() >> 11) * 0x1.0p-53; } while (u <= 0.0);
        return u;
    }
    // Rejection on the top bits, like R_unif_index.
    double unif_index(double n) {
        const auto m = static_cast<std::uint64_t>(n);
        if (m <= 1) return 0;
        int bits = 0;
        while (bits < 64 && (std::uint64_t(1) << bits) < m) ++bits;
        std::uint64_t v;
        do { v = gen() >> (64 - bits); } while (v >= m);
        return static_cast<double>(v);
    }
};

enum class RhoMethod { Standard, Fast };

// rhoR's rhoK() arguments after the observed kappa, base rate and test set
// length, with the same defaults.
struct RhoOptions {
    double test_set_baserate_inflation = 0;
    int    ocs_length          = 10000;   // OcSLength
    int    replicates          = 800;
    double scs_kappa_threshold = 0.9;     // ScSKappaThreshold
    double scs_kappa_min       = 0.40;    // ScSKappaMin
    double scs_precision_min   = 0.6;     // ScSPrecisionMin
    double scs_precision_max   = 1;       // ScSPrecisionMax
    RhoMethod method = RhoMethod::Standard;
};

struct KappaPrecision {
    double kappa = 0, precision = 0;
};

// Both rhoR paths recurse until a valid (kappa, precision) pair turns up and
// give up with this message; here the retries are capped instead.
inline constexpr int RHO_MAX_ATTEMPTS = 10000;
inline const char* const RHO_NO_VALID_SET =
    "Could not generate a valid set given ranges of kappa and precision";

// ── Standard method: rhoR's R implementation ─────────────────────────────────

namespace rho_standard {

// getR(): the recall that gives `kappa` at this base rate and precision.
QESTATS_R_ARITHMETIC_FN
inline double recall(double kappa, double baserate, double precision) {
    QESTATS_R_ARITHMETIC
    const double top = kappa * precision;
    return top / (2 * precision - 2 * baserate - kappa + 2 * baserate * kappa);
}

// checkBRPKcombo(): whether (base rate, precision, kappa) can form a set.
QESTATS_R_ARITHMETIC_FN
inline bool valid_combo(double baserate, double precision, double kappa) {
    QESTATS_R_ARITHMETIC
    const double right = (2 * baserate * kappa - 2 * baserate - kappa) / (kappa - 2);
    return precision > right;
}

// seq(from, to, length.out = n).
QESTATS_R_ARITHMETIC_FN
inline std::vector<double> seq(double from, double to, int n) {
    QESTATS_R_ARITHMETIC
    if (n <= 0) return {};
    if (n <= 2) {
        std::vector<double> v = {from, to};
        v.resize(n);
        return v;
    }
    if (from == to) return std::vector<double>(n, from);
    const double by = (to - from) / (n - 1);
    std::vector<double> v(n);
    v[0] = from;
    for (int i = 1; i < n - 1; ++i) v[i] = from + static_cast<double>(i) * by;
    v[n - 1] = to;
    return v;
}

// sample.int(n, k) without replacement, 0-based.
template <class Rng>
std::vector<std::size_t> sample_int(std::size_t n, std::size_t k, Rng& rng) {
    if (k > n)
        throw std::invalid_argument(
            "cannot take a sample larger than the population when 'replace = FALSE'");
    std::vector<std::size_t> out(k);
    if (k < 2) {
        for (std::size_t i = 0; i < k; ++i)
            out[i] = static_cast<std::size_t>(rng.unif_index(static_cast<double>(n)));
        return out;
    }
    std::vector<std::size_t> x(n);
    for (std::size_t i = 0; i < n; ++i) x[i] = i;
    for (std::size_t i = 0; i < k; ++i) {
        const auto j = static_cast<std::size_t>(rng.unif_index(static_cast<double>(n)));
        out[i] = x[j];
        x[j] = x[--n];
    }
    return out;
}

// genPKcombo(): one (kappa, precision) draw from the two grids. When the pair
// is invalid, precision is redrawn from the grid values above the minimum
// valid precision; when none are, kappa is redrawn too.
template <class Rng>
QESTATS_R_ARITHMETIC_FN
KappaPrecision kappa_precision(const std::vector<double>& kappas,
                               const std::vector<double>& precisions,
                               double baserate, Rng& rng) {
    QESTATS_R_ARITHMETIC
    for (int attempt = 0; attempt < RHO_MAX_ATTEMPTS; ++attempt) {
        const double kappa = kappas[sample_int(kappas.size(), 1, rng)[0]];
        double precision = precisions[sample_int(precisions.size(), 1, rng)[0]];
        if (valid_combo(baserate, precision, kappa)) return {kappa, precision};
        const double precision_min =
            (2 * baserate * kappa - 2 * baserate - kappa) / (kappa - 2);
        std::vector<double> above;
        for (double p : precisions) if (p > precision_min) above.push_back(p);
        if (above.empty()) continue;
        precision = above[sample_int(above.size(), 1, rng)[0]];
        return {kappa, precision};
    }
    throw std::runtime_error(RHO_NO_VALID_SET);
}

// generateKPs() with the flat distribution: kappa from a 100,000-point grid,
// precision from a 10,000-point grid.
template <class Rng>
std::vector<KappaPrecision>
generate_kps(int n, double baserate, double kappa_min, double kappa_max,
             double precision_min, double precision_max, Rng& rng) {
    const auto kappas = seq(kappa_min, kappa_max, 100000);
    const auto precisions = seq(precision_min, precision_max, 10000);
    std::vector<KappaPrecision> kps;
    kps.reserve(n > 0 ? n : 0);
    for (int i = 0; i < n; ++i) kps.push_back(kappa_precision(kappas, precisions, baserate, rng));
    return kps;
}

// prset() / pr.ct(): the full set's contingency table (counts). R's round()
// rounds halves to even.
QESTATS_R_ARITHMETIC_FN
inline ConfusionMatrix contingency_table(double precision, double recall,
                                         int length, double baserate) {
    QESTATS_R_ARITHMETIC
    const double gold1s = std::max(std::nearbyint(baserate * length), 1.0);
    const double gold0s = length - gold1s;
    const double tp = std::max(std::nearbyint(gold1s * recall), 1.0);
    const double fp = std::min(gold0s, std::max(std::nearbyint(tp * (1 - precision) / precision), 1.0));
    return {tp, fp, gold1s - tp, gold0s - fp};
}

// getHandSet(): kappa of a test set of `length` rows sampled from the full
// set laid out as prset() does (TP rows, FN rows, FP rows, TN rows). With an
// inflated base rate, ceiling(length * inflation) rows are first sampled from
// the first rater's positives.
template <class Rng>
QESTATS_R_ARITHMETIC_FN
double hand_kappa(const ConfusionMatrix& full, int length, double inflation, Rng& rng) {
    QESTATS_R_ARITHMETIC
    const auto tp = static_cast<std::size_t>(full.tp), fn = static_cast<std::size_t>(full.fn);
    const auto fp = static_cast<std::size_t>(full.fp), tn = static_cast<std::size_t>(full.tn);
    const std::size_t gold1s = tp + fn, rows = gold1s + fp + tn;
    double counts[4] = {0, 0, 0, 0};   // tp, fp, fn, tn
    auto add = [&](std::size_t row) {
        if (row < tp)                counts[0] += 1;
        else if (row < gold1s)       counts[2] += 1;
        else if (row < gold1s + fp)  counts[1] += 1;
        else                         counts[3] += 1;
    };
    const double positives = std::ceil(length * inflation);
    if (positives > static_cast<double>(gold1s))
        throw std::invalid_argument("Not enough positives in first rater to inflate to this level");
    if (positives > 0) {
        const auto k = static_cast<std::size_t>(positives);
        std::vector<bool> taken(rows, false);
        for (std::size_t row : sample_int(gold1s, k, rng)) { taken[row] = true; add(row); }
        std::vector<std::size_t> others;
        others.reserve(rows - k);
        for (std::size_t row = 0; row < rows; ++row) if (!taken[row]) others.push_back(row);
        for (std::size_t i : sample_int(others.size(), static_cast<std::size_t>(length) - k, rng))
            add(others[i]);
    } else {
        for (std::size_t row : sample_int(rows, static_cast<std::size_t>(length), rng)) add(row);
    }
    return kappa_from_counts(counts[0], counts[1], counts[2], counts[3]);
}

// mean(x, na.rm = TRUE): R sums in long double and adds a correction pass.
inline double r_mean(const std::vector<double>& x) {
    long double s = 0;
    std::size_t n = 0;
    for (double v : x) if (!std::isnan(v)) { s += v; ++n; }
    s /= static_cast<long double>(n);
    if (std::isfinite(static_cast<double>(s))) {
        long double t = 0;
        for (double v : x) if (!std::isnan(v)) t += (v - s);
        s += t / static_cast<long double>(n);
    }
    return static_cast<double>(s);
}

// getBootPvalue(): 1 below the distribution's mean, else the share of the
// distribution at or above `observed`.
inline double boot_pvalue(const std::vector<double>& distribution, double observed) {
    if (observed < r_mean(distribution)) return 1;
    std::size_t at_or_above = 0;
    for (double v : distribution) if (v >= observed) ++at_or_above;
    return static_cast<double>(at_or_above) / static_cast<double>(distribution.size());
}

// calcRho() with given (kappa, precision) pairs.
template <class Rng>
double rho(double observed_kappa, double baserate, int test_set_length,
           const std::vector<KappaPrecision>& kps, const RhoOptions& opts, Rng& rng) {
    const std::size_t replicates =
        std::min(kps.size(), static_cast<std::size_t>(std::max(opts.replicates, 0)));
    std::vector<double> kappas(replicates, 0.0);
    for (std::size_t i = 0; i < replicates; ++i) {
        const double r = recall(kps[i].kappa, baserate, kps[i].precision);
        const auto full = contingency_table(kps[i].precision, r, opts.ocs_length, baserate);
        kappas[i] = hand_kappa(full, test_set_length, opts.test_set_baserate_inflation, rng);
    }
    return boot_pvalue(kappas, observed_kappa);
}

} // namespace rho_standard

// ── Fast method: rhoR's C++ implementation ───────────────────────────────────
//
// Compiled as plain C++, like rhoR's src/rhor.cpp, so the compiler treats its
// arithmetic the way it treated rhoR's.

namespace rho_fast {

// recall(): the same formula as rho_standard::recall.
inline double recall(double kappa, double baserate, double precision) {
    const double top = kappa * precision;
    return top / (2 * precision - 2 * baserate - kappa + 2 * baserate * kappa);
}

// check_BRK_combo().
inline bool valid_combo(double baserate, double precision, double kappa) {
    const double right = ((2 * baserate * kappa) - (2 * baserate) - kappa) / (kappa - 2);
    return precision > right;
}

// Armadillo's regspace(start, delta, end).
inline std::vector<double> regspace(double start, double delta, double end) {
    if (((start < end) && (delta < 0)) || ((start > end) && (delta > 0)) || (delta == 0)) return {};
    const bool ascend = (start <= end);
    const double inc = (delta < 0) ? -delta : delta;
    const double m = (ascend ? (end - start) : (start - end)) / inc;
    const auto n = static_cast<std::size_t>(1 + static_cast<std::size_t>(std::floor(m)));
    // Two loops, as Armadillo writes them: a compiler that fuses
    // start + i * inc into an FMA then does it the same way here.
    std::vector<double> v(n);
    double* mem = v.data();
    if (ascend) {
        for (std::size_t i = 0; i < n; ++i) { mem[i] = double(start + double(i * inc)); }
    } else {
        for (std::size_t i = 0; i < n; ++i) { mem[i] = double(start - double(i * inc)); }
    }
    return v;
}

// One index from RcppArmadillo's sample(x, 1, FALSE): int(n * unif_rand()).
template <class Rng>
std::size_t sample_one(std::size_t n, Rng& rng) {
    const int n_orig = static_cast<int>(n);
    const int j = static_cast<int>(n_orig * rng.unif_rand());
    return static_cast<std::size_t>(j);
}

// Armadillo's randi<int>(distr_param(a, b)) on R's RNG (RcppArmadillo's
// Alt_R_RNG.h): runif(0, RAND_MAX) scaled to b - a + 1 values.
template <class Rng>
int randi(int a, int b, Rng& rng) {
    if (!(a <= b))
        throw std::logic_error("randi(): incorrect distribution parameters; a must be less than b");
    const double rand_max = static_cast<double>(RAND_MAX);
    const double draw = 0.0 + (rand_max - 0.0) * rng.unif_rand();
    if (a == 0 && b == RAND_MAX) return static_cast<int>(draw);
    const auto length = static_cast<unsigned long long>(b - a + 1);
    const double scale = double(length) / rand_max;
    return (std::min)(b, int(draw * scale) + a);
}

// find_valid_pk().
template <class Rng>
KappaPrecision kappa_precision(const std::vector<double>& kappas,
                               const std::vector<double>& precisions,
                               double baserate, Rng& rng) {
    for (int attempt = 0; attempt < RHO_MAX_ATTEMPTS; ++attempt) {
        const double kappa = kappas[sample_one(kappas.size(), rng)];
        double precision = precisions[sample_one(precisions.size(), rng)];
        if (!valid_combo(baserate, precision, kappa)) {
            const double precision_min =
                (2 * baserate * kappa - 2 * baserate - kappa) / (kappa - 2);
            std::vector<double> above;
            for (double p : precisions) if (p > precision_min) above.push_back(p);
            if (above.empty()) continue;
            precision = above[sample_one(above.size(), rng)];
        }
        return {kappa, precision};
    }
    throw std::runtime_error(RHO_NO_VALID_SET);
}

// generate_kp_list() with the flat distribution (distributionType 0); rhoR's
// bell-shaped option (distributionType 1) is not ported.
template <class Rng>
std::vector<KappaPrecision>
generate_kps(int n, double baserate, double kappa_min, double kappa_max,
             double precision_min, double precision_max, Rng& rng,
             long distribution_length = 10000) {
    const double kappa_step = (kappa_max - kappa_min) / (distribution_length - 1);
    const auto kappas = regspace(kappa_min, kappa_step, kappa_max);
    const double precision_step = (precision_max - precision_min) / (10000 - 1);
    const auto precisions = regspace(precision_min, precision_step, precision_max);
    if (kappas.empty() || precisions.empty())
        throw std::invalid_argument("generate_kps: empty kappa or precision range");
    std::vector<KappaPrecision> kps;
    kps.reserve(n > 0 ? n : 0);
    for (int i = 0; i < n; ++i) kps.push_back(kappa_precision(kappas, precisions, baserate, rng));
    return kps;
}

// contingency_table(): the full set's table (counts). C's round() rounds
// halves away from zero.
inline ConfusionMatrix contingency_table(double precision, double recall,
                                         int length, double baserate) {
    const int gold1s = static_cast<int>(std::max(std::round(baserate * length), 1.0));
    const int gold0s = length - gold1s;
    const int tp = static_cast<int>(std::max(std::round(gold1s * recall), 1.0));
    const int fp = static_cast<int>(std::min(static_cast<double>(gold0s),
                                             std::max(std::round(tp * (1 - precision) / precision), 1.0)));
    return {double(tp), double(fp), double(gold1s - tp), double(gold0s - fp)};
}

// sample_contingency_table(): draw n items without replacement from cells
// (a 1x2 row or a 2x2 table in column order: TP, FP, FN, TN), returning each
// draw's 0-based cell. Kept as rhoR has it: the draw spans sum + 1 values and
// the first cell takes the extra one, so an empty first cell can be drawn.
// The draw spans every cell given; draws past the third cell land in the
// fourth.
template <class Rng>
std::vector<int> sample_contingency_table(std::vector<int> cells, int n, Rng& rng) {
    if (cells.size() < 4) cells.resize(4, 0);
    std::vector<int> out(n > 0 ? n : 0);
    for (int w = 0; w < n; ++w) {
        int max_int = 0;
        for (int c : cells) max_int += c;
        const int s = randi(0, max_int, rng);
        const int cell = (s <= cells[0]) ? 0
                       : (s <= cells[0] + cells[1]) ? 1
                       : (s <= cells[0] + cells[1] + cells[2]) ? 2 : 3;
        out[w] = cell;
        cells[cell]--;
    }
    return out;
}

// getHand_ct(): the test set's table, drawn from the full set's table.
template <class Rng>
ConfusionMatrix hand_contingency_table(const ConfusionMatrix& full, int length,
                                       double inflation, Rng& rng) {
    const int positives = static_cast<int>(std::ceil(length * inflation));
    int tp = static_cast<int>(full.tp), fp = static_cast<int>(full.fp);
    int fn = static_cast<int>(full.fn), tn = static_cast<int>(full.tn);
    int counts[4] = {0, 0, 0, 0};   // column order: tp, fp, fn, tn
    if (positives > 0) {
        int pos_tp = 0, pos_fn = 0;
        for (int cell : sample_contingency_table({tp, fn}, positives, rng))
            (cell == 0 ? pos_tp : pos_fn)++;
        tp -= pos_tp;
        fn -= pos_fn;
        counts[0] += pos_tp;
        counts[2] += pos_fn;
    }
    for (int cell : sample_contingency_table({tp, fp, fn, tn}, length - positives, rng))
        counts[cell]++;
    return {double(counts[0]), double(counts[1]), double(counts[2]), double(counts[3])};
}

// kappa_ct(): kappa as rhoR's C++ computes it. Within an ulp of
// kappa_from_counts, but NaN (not 1) when one cell holds every item.
inline double kappa_ct(const ConfusionMatrix& ct) {
    const double a = ct.tp, c = ct.fn, b = ct.fp, d = ct.tn;
    const double size = a + b + c + d;
    const double p_zero = (a + d) / size;
    const double p_yes = ((a + b) / size) * ((a + c) / size);
    const double p_no = ((c + d) / size) * ((b + d) / size);
    const double p_e = (p_yes + p_no);
    return (p_zero - p_e) / (1 - p_e);
}

// Armadillo's mean(): two interleaved accumulators.
inline double arma_mean(const std::vector<double>& x) {
    const std::size_t n = x.size();
    if (n == 0) throw std::logic_error("mean(): object has no elements");
    double acc1 = 0, acc2 = 0;
    std::size_t j;
    for (j = 1; j < n; j += 2) { acc1 += x[j - 1]; acc2 += x[j]; }
    if ((j - 1) < n) acc1 += x[j - 1];
    const double mean = (acc1 + acc2) / double(n);
    if (std::isfinite(mean)) return mean;
    for (double v : x) if (!std::isfinite(v)) return mean;
    double r_mean = 0;
    for (std::size_t i = 0; i < n; ++i) r_mean = r_mean + (x[i] - r_mean) / double(i + 1);
    return r_mean;
}

// getBootPvalue_c().
inline double boot_pvalue(const std::vector<double>& distribution, double observed) {
    if (observed < arma_mean(distribution)) return 1.0;
    std::size_t at_or_above = 0;
    for (double v : distribution) if (v >= observed) ++at_or_above;
    return (at_or_above * 1.0) / distribution.size();
}

// calcRho_c() with given (kappa, precision) pairs.
template <class Rng>
double rho(double observed_kappa, double baserate, int test_set_length,
           const std::vector<KappaPrecision>& kps, const RhoOptions& opts, Rng& rng) {
    const std::size_t replicates =
        std::min(kps.size(), static_cast<std::size_t>(std::max(opts.replicates, 0)));
    std::vector<double> kappas(replicates);
    for (std::size_t i = 0; i < replicates; ++i) {
        const double r = recall(kps[i].kappa, baserate, kps[i].precision);
        const auto full = contingency_table(kps[i].precision, r, opts.ocs_length, baserate);
        kappas[i] = kappa_ct(hand_contingency_table(full, test_set_length,
                                                    opts.test_set_baserate_inflation, rng));
    }
    return boot_pvalue(kappas, observed_kappa);
}

} // namespace rho_fast

// ── rho ──────────────────────────────────────────────────────────────────────

// rhoK()'s argument checks, with its messages.
inline void check_rho_args(double observed_kappa, double baserate, const RhoOptions& o) {
    if ((observed_kappa > 1 || observed_kappa < 0) || (baserate > 1 || baserate < 0))
        throw std::invalid_argument("x and OcSBaserate must be between 0 and 1.");
    if (o.scs_kappa_threshold < 0 || o.scs_kappa_min < 0)
        throw std::invalid_argument("OcSBaserate, ScSKappaThreshold, and ScSKappaMin must be positive.");
    if (o.scs_kappa_threshold < o.scs_kappa_min)
        throw std::invalid_argument("ScSKappaThreshold must be greater than ScSKappaMin.");
    if ((o.scs_precision_min < 0 || o.scs_precision_min > 1) ||
        (o.scs_precision_max < 0 || o.scs_precision_max > 1))
        throw std::invalid_argument("ScSPrecisionMin and ScSPrecisionMax must be between 0 and 1.");
    if (o.scs_precision_max < o.scs_precision_min)
        throw std::invalid_argument("ScSPrecisionMax must be greater than ScSPrecisionMin.");
}

// The (kappa, precision) pairs rho() simulates: opts.replicates draws with
// kappa in [scs_kappa_min, scs_kappa_threshold] and precision in
// [scs_precision_min, scs_precision_max].
template <class Rng>
std::vector<KappaPrecision> rho_kps(double baserate, const RhoOptions& opts, Rng& rng) {
    if (opts.method == RhoMethod::Fast)
        return rho_fast::generate_kps(opts.replicates, baserate, opts.scs_kappa_min,
                                      opts.scs_kappa_threshold, opts.scs_precision_min,
                                      opts.scs_precision_max, rng);
    return rho_standard::generate_kps(opts.replicates, baserate, opts.scs_kappa_min,
                                      opts.scs_kappa_threshold, opts.scs_precision_min,
                                      opts.scs_precision_max, rng);
}

// rhoK(): rho for an observed kappa on a test set of `test_set_length` items
// from a population with this base rate. Only a class type is taken as an Rng,
// so an integer seed of any type picks the overload below.
template <class Rng, std::enable_if_t<std::is_class_v<Rng>, int> = 0>
double rho(double observed_kappa, double baserate, int test_set_length,
           const RhoOptions& opts, Rng& rng) {
    check_rho_args(observed_kappa, baserate, opts);
    const auto kps = rho_kps(baserate, opts, rng);
    if (opts.method == RhoMethod::Fast)
        return rho_fast::rho(observed_kappa, baserate, test_set_length, kps, opts, rng);
    return rho_standard::rho(observed_kappa, baserate, test_set_length, kps, opts, rng);
}

inline double rho(double observed_kappa, double baserate, int test_set_length,
                  const RhoOptions& opts = {}, std::uint64_t seed = 5489u) {
    StdRng rng(seed);
    return rho(observed_kappa, baserate, test_set_length, opts, rng);
}

} // namespace qestats

// [[Rcpp::depends(Rcpp)]]
// [[Rcpp::plugins(cpp17)]]
#include <Rcpp.h>
#include "qestats/qestats.hpp"

// The R binding mirrors qestats function signatures closely.
// Complex inputs/outputs use named R lists; JSON is NOT used here
// (unlike the WASM boundary) because Rcpp's native list/vector types
// are more idiomatic for R callers.

using namespace Rcpp;

// ── Helper: Metrics → R named list ───────────────────────────────────────────

static List metrics_to_rlist(const qestats::Metrics& m) {
    return List::create(
        Named("tp")        = m.tp,
        Named("fp")        = m.fp,
        Named("fn")        = m.fn,
        Named("tn")        = m.tn,
        Named("baserate")  = m.b1,
        Named("baserate2") = m.b2,
        Named("precision") = m.precision,
        Named("recall")    = m.recall,
        Named("accuracy")  = m.accuracy,
        Named("FDR")       = m.fdr,
        Named("FOR")       = m.for_,
        Named("kappa")     = m.kappa
    );
}

static List ps_to_rlist(const qestats::PsResult& r) {
    return List::create(Named("n") = r.n,
                        Named("power") = r.power,
                        Named("cost")  = r.cost);
}

// ── calculate_metrics ─────────────────────────────────────────────────────────
//' Calculate classification metrics from a confusion matrix
//' @param tp True positives
//' @param fp False positives
//' @param fn False negatives
//' @param tn True negatives
//' @return Named list of all classification metrics
//' @export
// [[Rcpp::export]]
List qestats_calculate_metrics(double tp, double fp, double fn, double tn) {
    return metrics_to_rlist(qestats::calculate_metrics(tp, fp, fn, tn));
}

// ── derive_from_rates ─────────────────────────────────────────────────────────
//' Derive all metrics from baserate + FDR + FOR (closed-form, fast)
//' @param b1       Classifier positive rate (baserate)
//' @param fdr      False Discovery Rate
//' @param for_rate False Omission Rate
//' @return Named list of all classification metrics
//' @export
// [[Rcpp::export]]
List qestats_derive_from_rates(double b1, double fdr, double for_rate) {
    return metrics_to_rlist(qestats::derive_from_rates(b1, fdr, for_rate));
}

// ── derive ────────────────────────────────────────────────────────────────────
//' Derive all metrics from three independent known metrics (general)
//' @param name1,val1  First metric name and value
//' @param name2,val2  Second metric name and value
//' @param name3,val3  Third metric name and value
//' @details Accepted names: "baserate","baserate2","b1","b2","precision",
//'   "recall","accuracy","FDR","FOR","kappa"
//' @return Named list of all classification metrics
//' @export
// [[Rcpp::export]]
List qestats_derive(const std::string& name1, double val1,
                    const std::string& name2, double val2,
                    const std::string& name3, double val3) {
    return metrics_to_rlist(qestats::derive(name1, val1, name2, val2, name3, val3));
}

// ── range ─────────────────────────────────────────────────────────────────────
//' Compute the [min, max] range of a third metric given two known metrics
//' @param name1,val1  First metric name and value
//' @param name2,val2  Second metric name and value
//' @param target      Name of the metric whose range to compute
//' @return Named list with \code{min} and \code{max}
//' @export
// [[Rcpp::export]]
List qestats_range(const std::string& name1, double val1,
                   const std::string& name2, double val2,
                   const std::string& target) {
    auto r = qestats::range(name1, val1, name2, val2, target);
    return List::create(Named("min") = r.min, Named("max") = r.max);
}

// ── phyper ────────────────────────────────────────────────────────────────────
//' Hypergeometric CDF P(X <= q) where X ~ Hypergeometric(N, K, n)
//' @param q         Upper bound
//' @param K         Successes in population
//' @param N_minus_K Failures in population
//' @param n         Sample size
//' @return Cumulative probability
//' @export
// [[Rcpp::export]]
double qestats_phyper(int q, int K, int N_minus_K, int n) {
    return qestats::phyper(q, K, N_minus_K, n);
}

// ── hyper_range ───────────────────────────────────────────────────────────────
//' Confidence interval for a proportion via hypergeometric distribution
//' @param n_pos Positives in sample
//' @param n     Sample size
//' @param N     Population size
//' @param alpha Significance level (default 0.05)
//' @return Named list with \code{lower}, \code{mean}, \code{upper}
//' @export
// [[Rcpp::export]]
List qestats_hyper_range(int n_pos, int n, int N, double alpha = 0.05) {
    auto iv = qestats::hyper_range(n_pos, n, N, alpha);
    return List::create(Named("lower") = iv.lower,
                        Named("mean")  = iv.mean,
                        Named("upper") = iv.upper);
}

// ── ps_proportion ─────────────────────────────────────────────────────────────
//' Perfect sample size for a proportion test
//' @param x0             Null hypothesis proportion (must be in [0.5, 1))
//' @param x1             Alternative proportion for power calculation (0 = auto)
//' @param N              Population size (0 = infinite)
//' @param allowed_errors Maximum allowable errors
//' @param alpha          Significance level (0.001–0.1)
//' @return Named list with \code{n}, \code{power}, \code{cost}
//' @export
// [[Rcpp::export]]
List qestats_ps_proportion(double x0, double x1 = 0, int N = 0,
                            int allowed_errors = 0, double alpha = 0.05) {
    return ps_to_rlist(qestats::ps_proportion(x0, x1, N, allowed_errors, alpha));
}

// ── ps_kappa ──────────────────────────────────────────────────────────────────
//' Perfect sample size for a kappa test
//' @param kappa0         Target kappa value (0, 1)
//' @param baserate       Base rate (proportion of positives)
//' @param kappa1         Alternative kappa for power calculation (0 = auto)
//' @param N              Population size (0 = infinite)
//' @param allowed_errors Maximum allowable errors
//' @param alpha          Significance level (0.001–0.1)
//' @return Named list with \code{n}, \code{power}, \code{cost}
//' @export
// [[Rcpp::export]]
List qestats_ps_kappa(double kappa0, double baserate, double kappa1 = 0,
                       int N = 0, int allowed_errors = 0, double alpha = 0.05) {
    return ps_to_rlist(qestats::ps_kappa(kappa0, baserate, kappa1, N, allowed_errors, alpha));
}

// ── ps_est_proportion ─────────────────────────────────────────────────────────
//' Estimate minimum proportion from a given perfect sample size
//' @param n             Perfect sample size
//' @param N             Population size (0 = infinite)
//' @param allowed_errors Maximum allowable errors
//' @param alpha          Significance level (0.001–0.1)
//' @return Lower bound for the proportion
//' @export
// [[Rcpp::export]]
double qestats_ps_est_proportion(int n, int N = 0,
                                  int allowed_errors = 0, double alpha = 0.05) {
    return qestats::ps_est_proportion(n, N, allowed_errors, alpha);
}

//' Estimate minimum kappa from a given perfect sample size
//' @param n             Perfect sample size
//' @param baserate      Classifier base rate
//' @param N             Population size (0 = infinite)
//' @param allowed_errors Maximum allowable errors
//' @param alpha          Significance level (0.001–0.1)
//' @return Lower bound for kappa
//' @export
// [[Rcpp::export]]
double qestats_ps_est_kappa(int n, double baserate, int N = 0,
                              int allowed_errors = 0, double alpha = 0.05) {
    return qestats::ps_est_kappa(n, baserate, N, allowed_errors, alpha);
}

// ── metrics_range ─────────────────────────────────────────────────────────────
//' Estimate confidence intervals for classification metrics
//' @param tp,fp,fn,tn Confusion matrix counts
//' @param N           Population size
//' @param baserate    Proportion of positives in population
//' @param alpha       Significance level (default 0.05)
//' @return Data frame with columns: metric, lower, mean, upper
//' @export
// [[Rcpp::export]]
DataFrame qestats_metrics_range(double tp, double fp, double fn, double tn,
                                 int N, double baserate, double alpha = 0.05) {
    auto rows = qestats::metrics_range(tp, fp, fn, tn, N, baserate, alpha);
    std::vector<std::string> metric;
    std::vector<double> lower, mean_v, upper;
    for (const auto& r : rows) {
        metric.push_back(r.metric);
        lower.push_back(r.lower);
        mean_v.push_back(r.mean);
        upper.push_back(r.upper);
    }
    return DataFrame::create(Named("metric") = metric,
                             Named("lower")  = lower,
                             Named("mean")   = mean_v,
                             Named("upper")  = upper);
}

// ── bootstrap_metrics ─────────────────────────────────────────────────────────
//' Bootstrap confidence intervals for classification metrics
//' @param cm         Named list with TP, FP, FN, TN counts
//' @param n_sample   Resample size
//' @param n_set      Population size (0 = infinite)
//' @param n_bootstrap Number of bootstrap replicates
//' @param conf_level Confidence level (default 0.95)
//' @return List with \code{original} (metrics list) and \code{stats} (data frame)
//' @export
// [[Rcpp::export]]
List qestats_bootstrap_metrics(List cm, int n_sample, int n_set = 0,
                                int n_bootstrap = 1000, double conf_level = 0.95) {
    qestats::ConfusionMatrix c{
        as<double>(cm["TP"]), as<double>(cm["FP"]),
        as<double>(cm["FN"]), as<double>(cm["TN"])
    };
    auto r = qestats::bootstrap_metrics(c, n_sample, {}, n_set, n_bootstrap, conf_level);

    std::vector<std::string> names;
    std::vector<double> se, ci_lo, ci_hi;
    for (const auto& s : r.stats) {
        names.push_back(s.name);
        se.push_back(s.se);
        ci_lo.push_back(s.ci_lo);
        ci_hi.push_back(s.ci_hi);
    }
    return List::create(
        Named("original") = metrics_to_rlist(r.original),
        Named("stats") = DataFrame::create(
            Named("metric") = names,
            Named("se")     = se,
            Named("ci_lo")  = ci_lo,
            Named("ci_hi")  = ci_hi)
    );
}

// ── p_diff ────────────────────────────────────────────────────────────────────
//' P-value for metric difference between two groups (bootstrap Welch t-test)
//' @param cm1,cm2    Named lists with TP, FP, FN, TN for each group
//' @param n_sample   Sample size (0 = sum of both groups)
//' @param metric     Metric name to compare
//' @param n_set      Population size (0 = infinite)
//' @param n_bootstrap Number of bootstrap replicates
//' @param alpha      Significance level
//' @return p-value (double)
//' @export
// [[Rcpp::export]]
double qestats_p_diff(List cm1, List cm2, int n_sample,
                       const std::string& metric, int n_set = 0,
                       int n_bootstrap = 1000, double alpha = 0.05) {
    qestats::ConfusionMatrix c1{
        as<double>(cm1["TP"]), as<double>(cm1["FP"]),
        as<double>(cm1["FN"]), as<double>(cm1["TN"])
    };
    qestats::ConfusionMatrix c2{
        as<double>(cm2["TP"]), as<double>(cm2["FP"]),
        as<double>(cm2["FN"]), as<double>(cm2["TN"])
    };
    return qestats::p_diff(c1, c2, n_sample, metric, n_set, n_bootstrap, alpha);
}

// ── N_Diff2 ───────────────────────────────────────────────────────────────────
//' Minimum n at which two groups' metric CIs no longer overlap
//' @param cm1,cm2  Named lists with TP, FP, FN, TN
//' @param metric   Metric name
//' @param N        Population size
//' @param d        Overlap threshold (default 0.1)
//' @param alpha    Significance level (default 0.05)
//' @return Required sample size (integer)
//' @export
// [[Rcpp::export]]
int qestats_N_Diff2(List cm1, List cm2, const std::string& metric,
                     int N, double d = 0.1, double alpha = 0.05) {
    qestats::ConfusionMatrix c1{
        as<double>(cm1["TP"]), as<double>(cm1["FP"]),
        as<double>(cm1["FN"]), as<double>(cm1["TN"])
    };
    qestats::ConfusionMatrix c2{
        as<double>(cm2["TP"]), as<double>(cm2["FP"]),
        as<double>(cm2["FN"]), as<double>(cm2["TN"])
    };
    return qestats::N_Diff2(c1, c2, metric, N, d, alpha);
}

// ── rho ───────────────────────────────────────────────────────────────────────
// R's RNG for qestats::rho: a seeded call returns what rhoR returns.
struct RRng {
    double unif_rand() { return ::unif_rand(); }
    double unif_index(double n) { return R_unif_index(n); }
};

//' Shaffer's rho for an observed kappa (rhoR's rhoK)
//' @param x                  Observed kappa of the test set
//' @param baserate           Base rate of the coded population (OcSBaserate)
//' @param test_set_length    Number of items in the test set
//' @param test_set_baserate_inflation Minimum share of positives in the test set
//' @param ocs_length         Size of each simulated population (OcSLength)
//' @param replicates         Number of simulated populations
//' @param scs_kappa_threshold Kappa threshold (ScSKappaThreshold)
//' @param scs_kappa_min      Minimum simulated kappa (ScSKappaMin)
//' @param scs_precision_min  Minimum simulated precision (ScSPrecisionMin)
//' @param scs_precision_max  Maximum simulated precision (ScSPrecisionMax)
//' @param method             "standard" (rhoR's default R path) or "fast"
//'   (rhoR's C++ path)
//' @return rho (double); uses R's random number generator
//' @export
// [[Rcpp::export]]
double qestats_rho(double x, double baserate, int test_set_length,
                   double test_set_baserate_inflation = 0, int ocs_length = 10000,
                   int replicates = 800, double scs_kappa_threshold = 0.9,
                   double scs_kappa_min = 0.40, double scs_precision_min = 0.6,
                   double scs_precision_max = 1, std::string method = "standard") {
    qestats::RhoOptions opts;
    opts.test_set_baserate_inflation = test_set_baserate_inflation;
    opts.ocs_length = ocs_length;
    opts.replicates = replicates;
    opts.scs_kappa_threshold = scs_kappa_threshold;
    opts.scs_kappa_min = scs_kappa_min;
    opts.scs_precision_min = scs_precision_min;
    opts.scs_precision_max = scs_precision_max;
    opts.method = (method == "standard") ? qestats::RhoMethod::Standard
                                         : qestats::RhoMethod::Fast;
    RRng rng;
    return qestats::rho(x, baserate, test_set_length, opts, rng);
}

1		#' @useDynLib qestats, .registration = TRUE
2		#' @importFrom Rcpp evalCpp
3		NULL
4
5		.onLoad <- function(libname, pkgname) {
6	!	invisible(NULL)
7		}

1		#pragma once
2		#include "metrics.hpp"
3		#include "inference.hpp"
4		#include "derivation.hpp"
5		#include <algorithm>
6		#include <cmath>
7		#include <random>
8		#include <stdexcept>
9		#include <string>
10		#include <vector>
11
12		namespace qestats {
13
14		// ── Bootstrap result ──────────────────────────────────────────────────────────
15
16		struct BootstrapResult {
17		Metrics original;
18		// Per-metric standard errors and CI (parallel to METRIC_NAMES order)
19		struct MetricStat {
20		std::string name;
21		double se = 0;
22		double ci_lo = 0;
23		double ci_hi = 0;
24		};
25		std::vector<MetricStat> stats;
26		};
27
28		// Metric names in canonical order matching qestats R output.
29	!	inline const std::vector<std::string>& metric_names() {
30		static const std::vector<std::string> v = {
31		"b1","b2","precision","recall","accuracy","FDR","FOR","kappa"
32		};
33	!	return v;
34		}
35
36		// ── bootstrap_metrics ─────────────────────────────────────────────────────────
37		// cm: raw confusion matrix (TP, FP, FN, TN counts — not proportions)
38		// n_sample: resample size
39		// n_set: population size (0 = infinite)
40		// n_bootstrap: number of bootstrap replicates
41		// conf_level: confidence level (e.g. 0.95)
42
43		inline BootstrapResult
44	!	bootstrap_metrics(const ConfusionMatrix& cm, int n_sample,
45		const std::vector<std::string>& metrics = {},
46		int n_set = 0, int n_bootstrap = 1000,
47		double conf_level = 0.95) {
48	!	const auto& all = metric_names();
49	!	const auto& chosen = metrics.empty() ? all : metrics;
50
51	!	double n1 = cm.tp, n2 = cm.fp, n3 = cm.fn, n4 = cm.tn;
52	!	double n0 = n1 + n2 + n3 + n4;
53	!	if (n0 <= 0 \|\| n1 < 0 \|\| n2 < 0 \|\| n3 < 0 \|\| n4 < 0)
54	!	throw std::domain_error("bootstrap_metrics: invalid confusion matrix");
55
56	!	BootstrapResult result;
57	!	result.original = calculate_metrics(n1, n2, n3, n4);
58
59		// Scale to population if requested
60	!	if (n_set > 0) {
61	!	double r = static_cast<double>(n_set) / n0;
62	!	if (r > 1) {
63	!	n1 = std::ceil(n1 * r); n2 = std::ceil(n2 * r);
64	!	n3 = std::ceil(n3 * r); n4 = std::ceil(n4 * r);
65		} else {
66	!	n1 = std::floor(n1 * r); n2 = std::floor(n2 * r);
67	!	n3 = std::floor(n3 * r); n4 = std::floor(n4 * r);
68		}
69	!	double tot = n1 + n2 + n3 + n4;
70	!	if (n_sample > tot) n_sample = static_cast<int>(tot);
71		}
72	!	const double n_total = n1 + n2 + n3 + n4;
73	!	const bool with_replacement = (n_set == 0);
74
75		// Build population index vector [0..n_total)
76	!	const int N = static_cast<int>(n_total);
77	!	const int c1 = static_cast<int>(n1);
78	!	const int c2 = static_cast<int>(n2);
79	!	const int c3 = static_cast<int>(n3);
80
81	!	std::mt19937 rng(42); // fixed seed for reproducibility
82
83		// Collect bootstrap replicates
84	!	std::vector<std::vector<double>> samples(chosen.size());
85
86	!	for (int b = 0; b < n_bootstrap; ++b) {
87	!	std::vector<int> idx(n_sample);
88	!	if (with_replacement) {
89	!	std::uniform_int_distribution<int> dist(0, N - 1);
90	!	for (int& x : idx) x = dist(rng);
91		} else {
92		// sample without replacement
93	!	std::vector<int> pool(N);
94	!	std::iota(pool.begin(), pool.end(), 0);
95	!	std::shuffle(pool.begin(), pool.end(), rng);
96	!	for (int i = 0; i < n_sample; ++i) idx[i] = pool[i];
97		}
98	!	int b1_ = 0, b2_ = 0, b3_ = 0, b4_ = 0;
99	!	for (int x : idx) {
100	!	if (x < c1) ++b1_;
101	!	else if (x < c1+c2) ++b2_;
102	!	else if (x < c1+c2+c3) ++b3_;
103	!	else ++b4_;
104		}
105	!	Metrics bm = calculate_metrics(b1_, b2_, b3_, b4_);
106	!	for (std::size_t i = 0; i < chosen.size(); ++i)
107	!	samples[i].push_back(get_metric(bm, chosen[i]));
108		}
109
110	!	const double alpha = 1.0 - conf_level;
111	!	for (std::size_t i = 0; i < chosen.size(); ++i) {
112	!	auto& s = samples[i];
113	!	double mean = 0; for (double x : s) mean += x; mean /= s.size();
114	!	double var = 0; for (double x : s) var += (x-mean)*(x-mean);
115	!	double se = std::sqrt(var / (s.size() - 1));
116
117	!	std::sort(s.begin(), s.end());
118	!	int lo_idx = std::max(0, static_cast<int>(std::floor(alpha/2 * s.size())));
119	!	int hi_idx = std::min(static_cast<int>(s.size())-1,
120	!	static_cast<int>(std::ceil((1.0-alpha/2) * s.size())));
121	!	result.stats.push_back({chosen[i], se, s[lo_idx], s[hi_idx]});
122		}
123	!	return result;
124		}
125
126		// ── Welch t-test helper ───────────────────────────────────────────────────────
127
128		struct TTestResult {
129		double t_stat = 0;
130		double df = 0;
131		double p_value = 1;
132		};
133
134	!	inline TTestResult welch_t(double m1, double m2,
135		double s1, double s2,
136		int n1, int n2,
137		double tol = 1e-9) {
138	!	if (std::abs(m1 - m2) < tol) return {0, 1, 1.0};
139	!	if (std::abs(s1) < tol && std::abs(s2) < tol)
140	!	return {0, 1, (std::abs(m1-m2) < tol) ? 1.0 : 0.0};
141
142	!	const double t = (m1 - m2) / std::sqrt(s1s1 + s2s2);
143	!	const double num = (s1s1 + s2s2) * (s1s1 + s2s2);
144	!	const double den = s1s1s1s1/(n1-1) + s2s2s2s2/(n2-1);
145	!	const double df = (den > 0) ? num/den : 1.0;
146
147		// Two-sided p-value via regularised incomplete beta (approx via normal for large df)
148		// For production accuracy use a proper t-distribution CDF; for now use
149		// normal approximation which is good for df > 30.
150		double p;
151	!	if (df > 30) {
152		// Standard normal approximation
153	!	double z = std::abs(t);
154	!	p = 2.0 * (0.5 * std::erfc(z / std::sqrt(2.0)));
155		} else {
156		// Abramowitz & Stegun approximation to t CDF
157	!	double x = df / (df + t*t);
158		// regularized incomplete beta I_x(df/2, 1/2) via continued fraction (simple approximation)
159		// Good enough for n_bootstrap >= 30; exact implementation would need boost or a special function
160	!	double z = std::abs(t);
161	!	double a_s = df / 2.0, b_s = 0.5;
162		// Use the normal approximation as a safe fallback
163	!	p = 2.0 * (0.5 * std::erfc(z / std::sqrt(2.0)));
164		}
165	!	return {t, df, p};
166		}
167
168		// ── p_diff ────────────────────────────────────────────────────────────────────
169		// p-value for a metric difference between two groups using bootstrap Welch t.
170
171	!	inline double p_diff(const ConfusionMatrix& cm1, const ConfusionMatrix& cm2,
172		int n_sample, const std::string& metric,
173		int n_set = 0, int n_bootstrap = 1000, double alpha = 0.05) {
174	!	const double n1 = cm1.tp + cm1.fp + cm1.fn + cm1.tn;
175	!	const double n2 = cm2.tp + cm2.fp + cm2.fn + cm2.tn;
176	!	const double n0 = n1 + n2;
177	!	if (n_sample <= 0) n_sample = static_cast<int>(n0);
178	!	const double r = static_cast<double>(n_sample) / n0;
179	!	const int ns1 = std::max(1, static_cast<int>(std::floor(n1 * (r < 1 ? 1 : r))));
180	!	const int ns2 = std::max(1, static_cast<int>(std::floor(n2 * (r < 1 ? 1 : r))));
181	!	const int nset1 = (n_set > 0) ? static_cast<int>(std::floor(n_set * n1 / n0)) : 0;
182	!	const int nset2 = (n_set > 0) ? static_cast<int>(std::floor(n_set * n2 / n0)) : 0;
183
184	!	auto r1 = bootstrap_metrics(cm1, ns1, {metric}, nset1, n_bootstrap, 1.0-alpha);
185	!	auto r2 = bootstrap_metrics(cm2, ns2, {metric}, nset2, n_bootstrap, 1.0-alpha);
186
187	!	const double m1v = get_metric(r1.original, metric);
188	!	const double m2v = get_metric(r2.original, metric);
189	!	const double se1 = r1.stats.empty() ? 0 : r1.stats[0].se;
190	!	const double se2 = r2.stats.empty() ? 0 : r2.stats[0].se;
191
192	!	return welch_t(m1v, m2v, se1, se2, n_bootstrap, n_bootstrap).p_value;
193		}
194
195		// ── n_diff ────────────────────────────────────────────────────────────────────
196		// Minimum n at which two groups show a significant metric difference.
197
198		inline int n_diff(const ConfusionMatrix& cm1, const ConfusionMatrix& cm2,
199		const std::string& metric, int n_set = 0,
200		int n_bootstrap = 1000, double alpha = 0.05) {
201		const int n1 = static_cast<int>(cm1.tp+cm1.fp+cm1.fn+cm1.tn);
202		const int n2 = static_cast<int>(cm2.tp+cm2.fp+cm2.fn+cm2.tn);
203		int n = n1 + n2;
204
205		double p0 = p_diff(cm1, cm2, n, metric, n_set, n_bootstrap, alpha);
206		if (p0 < alpha) return n;
207
208		int left = n, right = n_set > 0 ? n_set : n;
209		if (n_set <= 0) {
210		while (p0 >= alpha) {
211		left = right;
212		right *= 2;
213		p0 = p_diff(cm1, cm2, right, metric, n_set, n_bootstrap, alpha);
214		}
215		} else {
216		double p1 = p_diff(cm1, cm2, n_set, metric, n_set, n_bootstrap, alpha);
217		if (p1 >= alpha) return n_set;
218		}
219
220		while (right - left > 1) {
221		int mid = (left + right + 1) / 2;
222		p0 = p_diff(cm1, cm2, mid, metric, n_set, n_bootstrap, alpha);
223		if (p0 >= alpha) left = mid;
224		else right = mid;
225		}
226		return (left + right + 1) / 2;
227		}
228
229		// ── N_Diff2 ───────────────────────────────────────────────────────────────────
230		// Minimum n at which CI intervals for two groups no longer overlap by > d.
231
232	!	inline double interval_overlap(const MetricRange& r1, const MetricRange& r2) {
233	!	double span = std::min(r1.max - r1.min, r2.max - r2.min);
234	!	if (span <= 0) return 0.0;
235	!	double overlap = std::min(r1.max, r2.max) - std::max(r1.min, r2.min);
236	!	return (overlap <= 0) ? 0.0 : overlap / span;
237		}
238
239		// range_est: hyper-based CI for a metric given raw counts + population N
240	!	inline MetricRange range_est(const ConfusionMatrix& cm,
241		const std::string& metric, int N, double alpha = 0.05) {
242	!	const int tp = static_cast<int>(cm.tp), fp = static_cast<int>(cm.fp);
243	!	const int fn = static_cast<int>(cm.fn), tn = static_cast<int>(cm.tn);
244	!	auto prec_ci = hyper_range(tp, tp+fp, N, alpha);
245	!	auto for_ci = hyper_range(fn, fn+tn, N, alpha); // FOR = fn/(fn+tn)
246	!	double b1 = static_cast<double>(tp+fp) / (tp+fp+fn+tn);
247
248	!	auto eval_m = [&](double fdr, double for_v) -> double {
249	!	Metrics m = derive_from_rates(b1, fdr, for_v);
250	!	return get_metric(m, metric);
251		};
252
253		// Precision = 1 - FDR; CI is [lo, hi] for precision → FDR is [1-hi, 1-lo]
254	!	double fdr_lo = 1.0 - prec_ci.upper, fdr_hi = 1.0 - prec_ci.lower;
255	!	double for_lo = for_ci.lower, for_hi = for_ci.upper;
256
257	!	double v_lo = eval_m(fdr_hi, for_hi); // worst case
258	!	double v_hi = eval_m(fdr_lo, for_lo); // best case
259	!	if (v_lo > v_hi) std::swap(v_lo, v_hi);
260	!	return {v_lo, v_hi};
261		}
262
263	!	inline int N_Diff2(const ConfusionMatrix& cm1, const ConfusionMatrix& cm2,
264		const std::string& metric, int N, double d = 0.1,
265		double alpha = 0.05) {
266	!	const double n1 = cm1.tp+cm1.fp+cm1.fn+cm1.tn;
267	!	const double n2 = cm2.tp+cm2.fp+cm2.fn+cm2.tn;
268	!	const int n = static_cast<int>(n1 + n2);
269
270	!	auto scale = [](const ConfusionMatrix& cm, double r) -> ConfusionMatrix {
271	!	return {std::floor(cm.tpr), std::floor(cm.fpr),
272	!	std::floor(cm.fnr), std::floor(cm.tnr)};
273		};
274
275	!	auto check_overlap = [&](double x) -> bool {
276	!	auto s1 = scale(cm1, x); auto s2 = scale(cm2, x);
277	!	auto r1 = range_est(s1, metric, N, alpha);
278	!	auto r2 = range_est(s2, metric, N, alpha);
279	!	return interval_overlap(r1, r2) >= d;
280		};
281
282	!	if (!check_overlap(1.0)) return n; // already separable
283
284	!	const double x_max = static_cast<double>(N) / n;
285	!	if (check_overlap(x_max)) return N; // never separable at population size
286
287	!	double left = 1.0, right = x_max;
288	!	while (right - left > 1.0 / n) {
289	!	double mid = (left + right) / 2.0;
290	!	if (check_overlap(mid)) left = mid;
291	!	else right = mid;
292		}
293	!	return static_cast<int>(std::ceil(n * (left + right) / 2.0));
294		}
295
296		// ── ps_est_kappa ──────────────────────────────────────────────────────────────
297		// Given a perfect sample of size n and a baserate, estimate the minimum kappa.
298
299	!	inline double ps_est_kappa(int n, double baserate, int N = 0,
300		int allowed_errors = 0, double alpha = 0.05) {
301	!	double accuracy = ps_est_proportion(n, N, allowed_errors, alpha);
302	!	if (accuracy < 0) return -1;
303	!	MetricRange kappa_range = range("baserate", baserate, "accuracy", accuracy, "kappa");
304	!	return kappa_range.min;
305		}
306
307		// ── ps_kappa ─────────────────────────────────────────────────────────────────
308		// Perfect sample size for kappa test — delegates to ps_proportion via
309		// max accuracy from range(baserate, kappa).
310
311	!	inline PsResult ps_kappa(double kappa0, double baserate,
312		double kappa1 = 0, int N = 0,
313		int allowed_errors = 0, double alpha = 0.05) {
314	!	if (alpha < 0.001 \|\| alpha > 0.1) return {};
315	!	if (kappa0 <= 0 \|\| kappa0 >= 1) return {};
316
317	!	MetricRange acc_range = range("baserate", baserate, "kappa", kappa0, "accuracy");
318	!	return ps_proportion(acc_range.max, kappa1, N, allowed_errors, alpha);
319		}
320
321		} // namespace qestats

1		// [[Rcpp::depends(Rcpp)]]
2		// [[Rcpp::plugins(cpp17)]]
3		#include <Rcpp.h>
4		#include "qestats/qestats.hpp"
5
6		// The R binding mirrors qestats function signatures closely.
7		// Complex inputs/outputs use named R lists; JSON is NOT used here
8		// (unlike the WASM boundary) because Rcpp's native list/vector types
9		// are more idiomatic for R callers.
10
11		using namespace Rcpp;
12
13		// ── Helper: Metrics → R named list ───────────────────────────────────────────
14
15	3x	static List metrics_to_rlist(const qestats::Metrics& m) {
16		return List::create(
17	9x	Named("tp") = m.tp,
18	9x	Named("fp") = m.fp,
19	9x	Named("fn") = m.fn,
20	9x	Named("tn") = m.tn,
21	9x	Named("baserate") = m.b1,
22	9x	Named("baserate2") = m.b2,
23	9x	Named("precision") = m.precision,
24	9x	Named("recall") = m.recall,
25	9x	Named("accuracy") = m.accuracy,
26	9x	Named("FDR") = m.fdr,
27	9x	Named("FOR") = m.for_,
28	12x	Named("kappa") = m.kappa
29		);
30		}
31
32	!	static List ps_to_rlist(const qestats::PsResult& r) {
33	!	return List::create(Named("n") = r.n,
34	!	Named("power") = r.power,
35	!	Named("cost") = r.cost);
36		}
37
38		// ── calculate_metrics ─────────────────────────────────────────────────────────
39		//' Calculate classification metrics from a confusion matrix
40		//' @param tp True positives
41		//' @param fp False positives
42		//' @param fn False negatives
43		//' @param tn True negatives
44		//' @return Named list of all classification metrics
45		//' @export
46		// [[Rcpp::export]]
47	3x	List qestats_calculate_metrics(double tp, double fp, double fn, double tn) {
48	3x	return metrics_to_rlist(qestats::calculate_metrics(tp, fp, fn, tn));
49		}
50
51		// ── derive_from_rates ─────────────────────────────────────────────────────────
52		//' Derive all metrics from baserate + FDR + FOR (closed-form, fast)
53		//' @param b1 Classifier positive rate (baserate)
54		//' @param fdr False Discovery Rate
55		//' @param for_rate False Omission Rate
56		//' @return Named list of all classification metrics
57		//' @export
58		// [[Rcpp::export]]
59	!	List qestats_derive_from_rates(double b1, double fdr, double for_rate) {
60	!	return metrics_to_rlist(qestats::derive_from_rates(b1, fdr, for_rate));
61		}
62
63		// ── derive ────────────────────────────────────────────────────────────────────
64		//' Derive all metrics from three independent known metrics (general)
65		//' @param name1,val1 First metric name and value
66		//' @param name2,val2 Second metric name and value
67		//' @param name3,val3 Third metric name and value
68		//' @details Accepted names: "baserate","baserate2","b1","b2","precision",
69		//' "recall","accuracy","FDR","FOR","kappa"
70		//' @return Named list of all classification metrics
71		//' @export
72		// [[Rcpp::export]]
73	!	List qestats_derive(const std::string& name1, double val1,
74		const std::string& name2, double val2,
75		const std::string& name3, double val3) {
76	!	return metrics_to_rlist(qestats::derive(name1, val1, name2, val2, name3, val3));
77		}
78
79		// ── range ─────────────────────────────────────────────────────────────────────
80		//' Compute the [min, max] range of a third metric given two known metrics
81		//' @param name1,val1 First metric name and value
82		//' @param name2,val2 Second metric name and value
83		//' @param target Name of the metric whose range to compute
84		//' @return Named list with \code{min} and \code{max}
85		//' @export
86		// [[Rcpp::export]]
87	!	List qestats_range(const std::string& name1, double val1,
88		const std::string& name2, double val2,
89		const std::string& target) {
90	!	auto r = qestats::range(name1, val1, name2, val2, target);
91	!	return List::create(Named("min") = r.min, Named("max") = r.max);
92		}
93
94		// ── phyper ────────────────────────────────────────────────────────────────────
95		//' Hypergeometric CDF P(X <= q) where X ~ Hypergeometric(N, K, n)
96		//' @param q Upper bound
97		//' @param K Successes in population
98		//' @param N_minus_K Failures in population
99		//' @param n Sample size
100		//' @return Cumulative probability
101		//' @export
102		// [[Rcpp::export]]
103	!	double qestats_phyper(int q, int K, int N_minus_K, int n) {
104	!	return qestats::phyper(q, K, N_minus_K, n);
105		}
106
107		// ── hyper_range ───────────────────────────────────────────────────────────────
108		//' Confidence interval for a proportion via hypergeometric distribution
109		//' @param n_pos Positives in sample
110		//' @param n Sample size
111		//' @param N Population size
112		//' @param alpha Significance level (default 0.05)
113		//' @return Named list with \code{lower}, \code{mean}, \code{upper}
114		//' @export
115		// [[Rcpp::export]]
116	!	List qestats_hyper_range(int n_pos, int n, int N, double alpha = 0.05) {
117	!	auto iv = qestats::hyper_range(n_pos, n, N, alpha);
118	!	return List::create(Named("lower") = iv.lower,
119	!	Named("mean") = iv.mean,
120	!	Named("upper") = iv.upper);
121		}
122
123		// ── ps_proportion ─────────────────────────────────────────────────────────────
124		//' Perfect sample size for a proportion test
125		//' @param x0 Null hypothesis proportion (must be in [0.5, 1))
126		//' @param x1 Alternative proportion for power calculation (0 = auto)
127		//' @param N Population size (0 = infinite)
128		//' @param allowed_errors Maximum allowable errors
129		//' @param alpha Significance level (0.001–0.1)
130		//' @return Named list with \code{n}, \code{power}, \code{cost}
131		//' @export
132		// [[Rcpp::export]]
133	!	List qestats_ps_proportion(double x0, double x1 = 0, int N = 0,
134		int allowed_errors = 0, double alpha = 0.05) {
135	!	return ps_to_rlist(qestats::ps_proportion(x0, x1, N, allowed_errors, alpha));
136		}
137
138		// ── ps_kappa ──────────────────────────────────────────────────────────────────
139		//' Perfect sample size for a kappa test
140		//' @param kappa0 Target kappa value (0, 1)
141		//' @param baserate Base rate (proportion of positives)
142		//' @param kappa1 Alternative kappa for power calculation (0 = auto)
143		//' @param N Population size (0 = infinite)
144		//' @param allowed_errors Maximum allowable errors
145		//' @param alpha Significance level (0.001–0.1)
146		//' @return Named list with \code{n}, \code{power}, \code{cost}
147		//' @export
148		// [[Rcpp::export]]
149	!	List qestats_ps_kappa(double kappa0, double baserate, double kappa1 = 0,
150		int N = 0, int allowed_errors = 0, double alpha = 0.05) {
151	!	return ps_to_rlist(qestats::ps_kappa(kappa0, baserate, kappa1, N, allowed_errors, alpha));
152		}
153
154		// ── ps_est_proportion ─────────────────────────────────────────────────────────
155		//' Estimate minimum proportion from a given perfect sample size
156		//' @param n Perfect sample size
157		//' @param N Population size (0 = infinite)
158		//' @param allowed_errors Maximum allowable errors
159		//' @param alpha Significance level (0.001–0.1)
160		//' @return Lower bound for the proportion
161		//' @export
162		// [[Rcpp::export]]
163	!	double qestats_ps_est_proportion(int n, int N = 0,
164		int allowed_errors = 0, double alpha = 0.05) {
165	!	return qestats::ps_est_proportion(n, N, allowed_errors, alpha);
166		}
167
168		//' Estimate minimum kappa from a given perfect sample size
169		//' @param n Perfect sample size
170		//' @param baserate Classifier base rate
171		//' @param N Population size (0 = infinite)
172		//' @param allowed_errors Maximum allowable errors
173		//' @param alpha Significance level (0.001–0.1)
174		//' @return Lower bound for kappa
175		//' @export
176		// [[Rcpp::export]]
177	!	double qestats_ps_est_kappa(int n, double baserate, int N = 0,
178		int allowed_errors = 0, double alpha = 0.05) {
179	!	return qestats::ps_est_kappa(n, baserate, N, allowed_errors, alpha);
180		}
181
182		// ── metrics_range ─────────────────────────────────────────────────────────────
183		//' Estimate confidence intervals for classification metrics
184		//' @param tp,fp,fn,tn Confusion matrix counts
185		//' @param N Population size
186		//' @param baserate Proportion of positives in population
187		//' @param alpha Significance level (default 0.05)
188		//' @return Data frame with columns: metric, lower, mean, upper
189		//' @export
190		// [[Rcpp::export]]
191	!	DataFrame qestats_metrics_range(double tp, double fp, double fn, double tn,
192		int N, double baserate, double alpha = 0.05) {
193	!	auto rows = qestats::metrics_range(tp, fp, fn, tn, N, baserate, alpha);
194	!	std::vector<std::string> metric;
195	!	std::vector<double> lower, mean_v, upper;
196	!	for (const auto& r : rows) {
197	!	metric.push_back(r.metric);
198	!	lower.push_back(r.lower);
199	!	mean_v.push_back(r.mean);
200	!	upper.push_back(r.upper);
201		}
202	!	return DataFrame::create(Named("metric") = metric,
203	!	Named("lower") = lower,
204	!	Named("mean") = mean_v,
205	!	Named("upper") = upper);
206		}
207
208		// ── bootstrap_metrics ─────────────────────────────────────────────────────────
209		//' Bootstrap confidence intervals for classification metrics
210		//' @param cm Named list with TP, FP, FN, TN counts
211		//' @param n_sample Resample size
212		//' @param n_set Population size (0 = infinite)
213		//' @param n_bootstrap Number of bootstrap replicates
214		//' @param conf_level Confidence level (default 0.95)
215		//' @return List with \code{original} (metrics list) and \code{stats} (data frame)
216		//' @export
217		// [[Rcpp::export]]
218	!	List qestats_bootstrap_metrics(List cm, int n_sample, int n_set = 0,
219		int n_bootstrap = 1000, double conf_level = 0.95) {
220		qestats::ConfusionMatrix c{
221	!	as<double>(cm["TP"]), as<double>(cm["FP"]),
222	!	as<double>(cm["FN"]), as<double>(cm["TN"])
223		};
224	!	auto r = qestats::bootstrap_metrics(c, n_sample, {}, n_set, n_bootstrap, conf_level);
225
226	!	std::vector<std::string> names;
227	!	std::vector<double> se, ci_lo, ci_hi;
228	!	for (const auto& s : r.stats) {
229	!	names.push_back(s.name);
230	!	se.push_back(s.se);
231	!	ci_lo.push_back(s.ci_lo);
232	!	ci_hi.push_back(s.ci_hi);
233		}
234		return List::create(
235	!	Named("original") = metrics_to_rlist(r.original),
236	!	Named("stats") = DataFrame::create(
237	!	Named("metric") = names,
238	!	Named("se") = se,
239	!	Named("ci_lo") = ci_lo,
240	!	Named("ci_hi") = ci_hi)
241		);
242		}
243
244		// ── p_diff ────────────────────────────────────────────────────────────────────
245		//' P-value for metric difference between two groups (bootstrap Welch t-test)
246		//' @param cm1,cm2 Named lists with TP, FP, FN, TN for each group
247		//' @param n_sample Sample size (0 = sum of both groups)
248		//' @param metric Metric name to compare
249		//' @param n_set Population size (0 = infinite)
250		//' @param n_bootstrap Number of bootstrap replicates
251		//' @param alpha Significance level
252		//' @return p-value (double)
253		//' @export
254		// [[Rcpp::export]]
255	!	double qestats_p_diff(List cm1, List cm2, int n_sample,
256		const std::string& metric, int n_set = 0,
257		int n_bootstrap = 1000, double alpha = 0.05) {
258		qestats::ConfusionMatrix c1{
259	!	as<double>(cm1["TP"]), as<double>(cm1["FP"]),
260	!	as<double>(cm1["FN"]), as<double>(cm1["TN"])
261		};
262		qestats::ConfusionMatrix c2{
263	!	as<double>(cm2["TP"]), as<double>(cm2["FP"]),
264	!	as<double>(cm2["FN"]), as<double>(cm2["TN"])
265		};
266	!	return qestats::p_diff(c1, c2, n_sample, metric, n_set, n_bootstrap, alpha);
267		}
268
269		// ── N_Diff2 ───────────────────────────────────────────────────────────────────
270		//' Minimum n at which two groups' metric CIs no longer overlap
271		//' @param cm1,cm2 Named lists with TP, FP, FN, TN
272		//' @param metric Metric name
273		//' @param N Population size
274		//' @param d Overlap threshold (default 0.1)
275		//' @param alpha Significance level (default 0.05)
276		//' @return Required sample size (integer)
277		//' @export
278		// [[Rcpp::export]]
279	!	int qestats_N_Diff2(List cm1, List cm2, const std::string& metric,
280		int N, double d = 0.1, double alpha = 0.05) {
281		qestats::ConfusionMatrix c1{
282	!	as<double>(cm1["TP"]), as<double>(cm1["FP"]),
283	!	as<double>(cm1["FN"]), as<double>(cm1["TN"])
284		};
285		qestats::ConfusionMatrix c2{
286	!	as<double>(cm2["TP"]), as<double>(cm2["FP"]),
287	!	as<double>(cm2["FN"]), as<double>(cm2["TN"])
288		};
289	!	return qestats::N_Diff2(c1, c2, metric, N, d, alpha);
290		}
291
292		// ── rho ───────────────────────────────────────────────────────────────────────
293		// R's RNG for qestats::rho: a seeded call returns what rhoR returns.
294		struct RRng {
295	246829x	double unif_rand() { return ::unif_rand(); }
296	248477x	double unif_index(double n) { return R_unif_index(n); }
297		};
298
299		//' Shaffer's rho for an observed kappa (rhoR's rhoK)
300		//' @param x Observed kappa of the test set
301		//' @param baserate Base rate of the coded population (OcSBaserate)
302		//' @param test_set_length Number of items in the test set
303		//' @param test_set_baserate_inflation Minimum share of positives in the test set
304		//' @param ocs_length Size of each simulated population (OcSLength)
305		//' @param replicates Number of simulated populations
306		//' @param scs_kappa_threshold Kappa threshold (ScSKappaThreshold)
307		//' @param scs_kappa_min Minimum simulated kappa (ScSKappaMin)
308		//' @param scs_precision_min Minimum simulated precision (ScSPrecisionMin)
309		//' @param scs_precision_max Maximum simulated precision (ScSPrecisionMax)
310		//' @param method "standard" (rhoR's default R path) or "fast"
311		//' (rhoR's C++ path)
312		//' @return rho (double); uses R's random number generator
313		//' @export
314		// [[Rcpp::export]]
315	11x	double qestats_rho(double x, double baserate, int test_set_length,
316		double test_set_baserate_inflation = 0, int ocs_length = 10000,
317		int replicates = 800, double scs_kappa_threshold = 0.9,
318		double scs_kappa_min = 0.40, double scs_precision_min = 0.6,
319		double scs_precision_max = 1, std::string method = "standard") {
320	11x	qestats::RhoOptions opts;
321	11x	opts.test_set_baserate_inflation = test_set_baserate_inflation;
322	11x	opts.ocs_length = ocs_length;
323	11x	opts.replicates = replicates;
324	11x	opts.scs_kappa_threshold = scs_kappa_threshold;
325	11x	opts.scs_kappa_min = scs_kappa_min;
326	11x	opts.scs_precision_min = scs_precision_min;
327	11x	opts.scs_precision_max = scs_precision_max;
328	11x	opts.method = (method == "standard") ? qestats::RhoMethod::Standard
329		: qestats::RhoMethod::Fast;
330		RRng rng;
331	19x	return qestats::rho(x, baserate, test_set_length, opts, rng);
332		}