diff --git a/src/radmeth/radmeth.cpp b/src/radmeth/radmeth.cpp index 9ace1fb4..789304c5 100644 --- a/src/radmeth/radmeth.cpp +++ b/src/radmeth/radmeth.cpp @@ -23,8 +23,6 @@ #include "smithlab_os.hpp" #include "smithlab_utils.hpp" -#include // GSL header for chisqrd distribution - #include #include #include @@ -32,6 +30,7 @@ #include #include #include +#include #include #include #include @@ -78,7 +77,7 @@ struct file_progress { } }; -static bool +[[nodiscard]] static bool consistent_sample_names(const Regression ®, const std::string &header) { std::istringstream iss(header); auto nm_itr(std::begin(reg.design.sample_names)); @@ -90,10 +89,83 @@ consistent_sample_names(const Regression ®, const std::string &header) { return true; } +[[nodiscard]] static std::vector +get_sample_names_from_header(const std::string &header) { + std::istringstream iss(header); + std::string token; + std::vector sample_names; + while (iss >> token) + sample_names.push_back(token); + return sample_names; +} + +// Series representation for the lower incomplete gamma P(a,x) +[[nodiscard]] static double +gamma_p_series(const double a, const double x) { + constexpr auto eps = std::numeric_limits::epsilon(); + constexpr auto max_iter = 100; + double sum = 1.0 / a; + double term = sum; + for (auto n = 1; n < max_iter; ++n) { + term *= x / (a + n); + sum += term; + if (term < eps * sum) + break; + } + return sum * std::exp(-x + a * std::log(x) - std::lgamma(a)); +} + +[[nodiscard]] static inline double +safe_floor(const double x, const double floor_val) { + return std::abs(x) < floor_val ? floor_val : x; +} + +// Continued fraction representation for the upper incomplete gamma Q(a,x) +[[nodiscard]] static double +gamma_q_contfrac(const double a, const double x) { + constexpr auto epsilon = std::numeric_limits::epsilon(); + constexpr auto fpmin = std::numeric_limits::min() / epsilon; + constexpr auto max_iter = 100; + double b = x + 1 - a; + double c = 1.0 / fpmin; + double d = 1.0 / b; + double h = d; + for (auto i = 1; i < max_iter; ++i) { + const double an = -i * (i - a); + b += 2; + d = safe_floor(an * d + b, fpmin); + c = safe_floor(b + an / c, fpmin); + d = 1.0 / d; + const double delta = d * c; + h *= delta; + if (std::abs(delta - 1.0) < epsilon) + break; + } + return std::exp(-x + a * std::log(x) - std::lgamma(a)) * h; +} + +// Regularized lower incomplete gamma P(a,x) +[[nodiscard]] static double +gamma_p(const double a, const double x) { + if (x < 0 || a <= 0) + return 0.0; + if (x == 0) + return 0.0; + if (x < a + 1.0) + return gamma_p_series(a, x); + return 1.0 - gamma_q_contfrac(a, x); +} + +// chi-square CDF: P(k/2, x/2) +[[nodiscard]] static double +chi_square_cdf(const double x, const double k) { + return gamma_p(k * 0.5, x * 0.5); +} + // Given the maximum likelihood estimates of the full and reduced models, the // function outputs the p-value of the log-likelihood ratio. *Note* that it is // assumed that the reduced model has one fewer factor than the reduced model. -static double +[[nodiscard]] static double llr_test(const double null_loglik, const double full_loglik) { // The log-likelihood ratio statistic. const double log_lik_stat = -2 * (null_loglik - full_loglik); @@ -103,7 +175,8 @@ llr_test(const double null_loglik, const double full_loglik) { const std::size_t degrees_of_freedom = 1; // Log-likelihood ratio statistic has a chi-sqare distribution. - const double chisq_p = gsl_cdf_chisq_P(log_lik_stat, degrees_of_freedom); + const double chisq_p = chi_square_cdf(log_lik_stat, degrees_of_freedom); + const double p_value = 1.0 - chisq_p; return p_value; @@ -195,6 +268,7 @@ radmeth(const bool show_progress, const bool more_na_info, std::string sample_names_header; std::getline(table_file, sample_names_header); + const auto sample_names = get_sample_names_from_header(sample_names_header); if (!consistent_sample_names(alt_model, sample_names_header)) { // clang-format off const auto message = @@ -350,6 +424,19 @@ that the design matrix and the proportion table are correctly formatted. progress(table_file); } +static void +ensure_sample_order(const std::string &table_filename, Regression &alt_model, + Regression &null_model) { + std::ifstream table_file(table_filename); + if (!table_file) + throw std::runtime_error("could not open file: " + table_filename); + std::string header; + std::getline(table_file, header); + const auto sample_names = get_sample_names_from_header(header); + alt_model.order_samples(sample_names); + null_model.order_samples(sample_names); +} + int main_radmeth(int argc, char *argv[]) { try { @@ -435,6 +522,8 @@ main_radmeth(int argc, char *argv[]) { if (verbose) std::cerr << "Null model:\n" << null_model.design << '\n'; + ensure_sample_order(table_filename, alt_model, null_model); + // clang-format off if (verbose) std::cerr << "Output columns:\n" diff --git a/src/radmeth/radmeth_model.cpp b/src/radmeth/radmeth_model.cpp index 5ca694a6..0df305cc 100644 --- a/src/radmeth/radmeth_model.cpp +++ b/src/radmeth/radmeth_model.cpp @@ -22,6 +22,7 @@ #include #include #include +#include #include double Regression::tolerance = 1e-4; @@ -167,6 +168,47 @@ Design::drop_factor(const std::size_t factor_idx) { return design; } +void +Design::order_samples(const std::vector &ordered_names) { + // Build lookup: sample name -> original index + const auto sample_to_index = [&] { + std::unordered_map sample_to_index; + std::uint32_t index = 0; + for (const auto &sample : sample_names) + sample_to_index.emplace(sample, index++); + return sample_to_index; + }(); + + std::vector ord_sample_names; + std::vector> ord_matrix; + std::vector ord_group_id; + // factor names should not change + // groups should not change + // tmatrix will be changed after matrix + + const auto n_names = std::size(ordered_names); + ord_sample_names.reserve(n_names); + ord_matrix.reserve(n_names); + ord_group_id.reserve(n_names); + + for (const auto &name : ordered_names) { + const auto it = sample_to_index.find(name); + if (it == std::cend(sample_to_index)) + throw std::runtime_error("Sample not found: " + name); + + const auto original_index = it->second; + ord_sample_names.push_back(sample_names[original_index]); + ord_matrix.push_back(matrix[original_index]); + ord_group_id.push_back(group_id[original_index]); + } + sample_names = std::move(ord_sample_names); + matrix = std::move(ord_matrix); + group_id = std::move(ord_group_id); + + // update tmatrix using matrix + transpose(matrix, tmatrix); +} + std::ostream & operator<<(std::ostream &out, const Design &design) { static constexpr std::uint32_t max_samples_to_report = 20; diff --git a/src/radmeth/radmeth_model.hpp b/src/radmeth/radmeth_model.hpp index 1c42cb17..fb4b5ac9 100644 --- a/src/radmeth/radmeth_model.hpp +++ b/src/radmeth/radmeth_model.hpp @@ -34,10 +34,10 @@ struct cumul_counts { struct Design { std::vector factor_names; std::vector sample_names; - std::vector> matrix; - std::vector> tmatrix; - std::vector> groups; - std::vector group_id; + std::vector> matrix; // samples=rows, factors=cols + std::vector> tmatrix; // factors=rows, samples=cols + std::vector> groups; // combs of fact levels + std::vector group_id; // assign group to sample [[nodiscard]] std::size_t n_factors() const { @@ -56,6 +56,9 @@ struct Design { [[nodiscard]] Design drop_factor(const std::size_t factor_idx); + + void + order_samples(const std::vector &ordered_names); }; std::istream & @@ -121,6 +124,11 @@ struct Regression { n_samples() const { return design.n_samples(); } + + void + order_samples(const std::vector &ordered_names) { + design.order_samples(ordered_names); + } }; #endif