radmeth sample consistency between design and data matrix

src/radmeth/radmeth.cpp

@@ -23,15 +23,14 @@
 #include "smithlab_os.hpp"
 #include "smithlab_utils.hpp"
 
-#include <gsl/gsl_cdf.h>  // GSL header for chisqrd distribution
-
 #include <algorithm>
 #include <chrono>
 #include <cmath>
 #include <filesystem>
 #include <fstream>
 #include <iomanip>
 #include <iostream>
+#include <limits>
 #include <sstream>
 #include <string>
 #include <thread>
@@ -78,7 +77,7 @@ struct file_progress {
   }
 };
 
-static bool
+[[nodiscard]] static bool
 consistent_sample_names(const Regression &reg, 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 &reg, const std::string &header) {
   return true;
 }
 
+[[nodiscard]] static std::vector<std::string>
+get_sample_names_from_header(const std::string &header) {
+  std::istringstream iss(header);
+  std::string token;
+  std::vector<std::string> 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<double>::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<double>::epsilon();
+  constexpr auto fpmin = std::numeric_limits<double>::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"

src/radmeth/radmeth_model.cpp

@@ -22,6 +22,7 @@
 #include <ostream>
 #include <sstream>
 #include <string>
+#include <unordered_map>
 #include <vector>
 
 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<std::string> &ordered_names) {
+  // Build lookup: sample name -> original index
+  const auto sample_to_index = [&] {
+    std::unordered_map<std::string, std::uint32_t> 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<std::string> ord_sample_names;
+  std::vector<std::vector<std::uint8_t>> ord_matrix;
+  std::vector<std::uint32_t> 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;

src/radmeth/radmeth_model.hpp

@@ -34,10 +34,10 @@ struct cumul_counts {
 struct Design {
   std::vector<std::string> factor_names;
   std::vector<std::string> sample_names;
-  std::vector<std::vector<std::uint8_t>> matrix;
-  std::vector<std::vector<std::uint8_t>> tmatrix;
-  std::vector<std::vector<std::uint8_t>> groups;
-  std::vector<std::uint32_t> group_id;
+  std::vector<std::vector<std::uint8_t>> matrix;   // samples=rows, factors=cols
+  std::vector<std::vector<std::uint8_t>> tmatrix;  // factors=rows, samples=cols
+  std::vector<std::vector<std::uint8_t>> groups;   // combs of fact levels
+  std::vector<std::uint32_t> 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<std::string> &ordered_names);
 };
 
 std::istream &
@@ -121,6 +124,11 @@ struct Regression {
   n_samples() const {
     return design.n_samples();
   }
+
+  void
+  order_samples(const std::vector<std::string> &ordered_names) {
+    design.order_samples(ordered_names);
+  }
 };
 
 #endif