From c99c3eb5ae09ed4cca91674539dd7ecd231a5fd1 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 12 Oct 2025 19:11:53 -0700 Subject: [PATCH 01/33] src/radmeth/radmeth_optimize.cpp: added an implementation for digamma functions via Chebyschev polynomial approximation for the gradient calculations. Number of coeffs to use is still to be determined. --- src/radmeth/radmeth_optimize.cpp | 84 ++++++++++++++++++++++++++++++-- 1 file changed, 81 insertions(+), 3 deletions(-) diff --git a/src/radmeth/radmeth_optimize.cpp b/src/radmeth/radmeth_optimize.cpp index 6c6e2621..d2b90bdf 100644 --- a/src/radmeth/radmeth_optimize.cpp +++ b/src/radmeth/radmeth_optimize.cpp @@ -21,10 +21,79 @@ #include #include +#include #include #include #include +/* Coefficients for the Chebyschev polynomial for the digamma function in the + range 0-1 */ +// clang-format off +static constexpr std::array psi1_cs { + -0.038057080835217922, // == -0.019028540417608961*2 + 0.491415393029387130, + -0.056815747821244730, + 0.008357821225914313, + -0.001333232857994342, + 0.000220313287069308, + -0.000037040238178456, + 0.000006283793654854, + -0.000001071263908506, + 0.000000183128394654, + -0.000000031353509361, + 0.000000005372808776, + -0.000000000921168141, + 0.000000000157981265, + -0.000000000027098646, + 0.000000000004648722, + -0.000000000000797527, + 0.000000000000136827, + -0.000000000000023475, + 0.000000000000004027, + -0.000000000000000691, + 0.000000000000000118, + -0.000000000000000020 +}; +// clang-format on + +/* Alternate set of coefficients for the Chebyschev polynomial for the digamma + function */ +// clang-format off +static constexpr std::array psi2_cs = { + -0.0204749044678185, // == -0.01023745223390925*2 + -0.0101801271534859, + 0.0000559718725387, + -0.0000012917176570, + 0.0000000572858606, + -0.0000000038213539, + 0.0000000003397434, + -0.0000000000374838, + 0.0000000000048990, + -0.0000000000007344, + 0.0000000000001233, + -0.0000000000000228, + 0.0000000000000045, + -0.0000000000000009, + 0.0000000000000002, + -0.0000000000000000 , +}; +// clang-format on + +template +[[nodiscard]] static inline double +chebyschev(const std::array &coeffs, std::uint32_t order, + const double y) { + const auto y2 = 2.0 * y; + double d = 0.0; + double dd = 0.0; + for (auto j = order; j >= 1; j--) { + const auto temp = d; + d = y2 * d - dd + coeffs[j]; + dd = temp; + } + return y * d - dd + 0.5 * coeffs[0]; +} + [[nodiscard]] static inline double logistic(const double x) { return 1.0 / (1.0 / std::exp(x) + 1.0); @@ -84,9 +153,18 @@ log_likelihood(const gsl_vector *params, Regression ®) { return ll; } -[[nodiscard]] static inline double -digamma(const double x) { - return gsl_sf_psi(x); +// digamma for x non-negative +static double +digamma(const double y) { + static constexpr auto psi_order = 7; // max=22; + static constexpr auto apsi_order = 7; // max=15; + if (y >= 2.0) { + const auto t = 8.0 / (y * y) - 1.0; + return std::log(y) - 0.5 / y + chebyschev(psi2_cs, apsi_order, t); + } + if (y < 1.0) + return -1.0 / y + chebyschev(psi1_cs, psi_order, 2.0 * y - 1.0); + return chebyschev(psi1_cs, psi_order, 2.0 * (y - 1.0) - 1.0); } static void From 72374ed15b8a85322846b15842b5839a2005d76b Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 15 Oct 2025 16:13:33 -0700 Subject: [PATCH 02/33] Makefile.am: adding additional sources --- Makefile.am | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/Makefile.am b/Makefile.am index af091f3c..e24502cd 100644 --- a/Makefile.am +++ b/Makefile.am @@ -190,7 +190,11 @@ libdnmtools_a_SOURCES += \ # optimization libdnmtools_a_SOURCES += \ src/radmeth/radmeth_optimize.hpp \ - src/radmeth/radmeth_model.hpp + src/radmeth/radmeth_optimize_nano.hpp \ + src/radmeth/radmeth_model.hpp \ + src/radmeth/radmeth_utils.hpp \ + src/radmeth/radmeth_design.hpp \ + src/radmeth/radmeth_model_nano.hpp LDADD = libdnmtools.a src/abismal/libabismal.a src/smithlab_cpp/libsmithlab_cpp.a @@ -236,9 +240,14 @@ dnmtools_SOURCES += src/amrfinder/amrtester.cpp dnmtools_SOURCES += src/radmeth/dmr.cpp dnmtools_SOURCES += src/radmeth/methdiff.cpp +dnmtools_SOURCES += src/radmeth/radmeth_utils.cpp dnmtools_SOURCES += src/radmeth/radmeth_optimize.cpp +dnmtools_SOURCES += src/radmeth/radmeth_optimize_nano.cpp +dnmtools_SOURCES += src/radmeth/radmeth_design.cpp dnmtools_SOURCES += src/radmeth/radmeth_model.cpp +dnmtools_SOURCES += src/radmeth/radmeth_model_nano.cpp dnmtools_SOURCES += src/radmeth/radmeth.cpp +dnmtools_SOURCES += src/radmeth/radmeth_nano.cpp dnmtools_SOURCES += src/radmeth/radmeth-adjust.cpp dnmtools_SOURCES += src/radmeth/radmeth-merge.cpp From 81bd7adfc0f92be14aca764769584a5c9d63508a Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 15 Oct 2025 16:14:03 -0700 Subject: [PATCH 03/33] src/dnmtools.cpp: adding the radmeth-nano command for doing DMPs on nanopore --- src/dnmtools.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/dnmtools.cpp b/src/dnmtools.cpp index b4bb81fd..4c56673c 100644 --- a/src/dnmtools.cpp +++ b/src/dnmtools.cpp @@ -111,6 +111,8 @@ main_radmeth_adjust(int argc, char *argv[]); int main_radmeth(int argc, char *argv[]); int +main_radmeth_nano(int argc, char *argv[]); +int main_radmeth_merge(int argc, char *argv[]); int main_clean_hairpins(int argc, char *argv[]); @@ -204,6 +206,7 @@ main(int argc, char *argv[]) { {{{"dmr", "identify DMRs from genomic intervals and single-CpG DM probabilities", main_dmr}, {"diff", "compute single-CpG DM probability between two methylomes", main_methdiff}, {"radmeth", "compute DM probabilities for each CpG using multiple methylomes", main_radmeth}, + {"radmeth-nano", "radmeth designed for nanopore data", main_radmeth_nano}, {"radadjust", "adjust p-values from radmeth output", main_radmeth_adjust}, {"radmerge", "merge significant CpGs in radmeth output", main_radmeth_merge}}}}, From bd30df2bf8c6e15b77ca6bd70ce7818e0745bfba Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 15 Oct 2025 16:15:56 -0700 Subject: [PATCH 04/33] src/radmeth/radmeth_design.hcpp: adding files for the Design class separated from the Regression classes --- src/radmeth/radmeth_design.cpp | 262 +++++++++++++++++++++++++++++++++ src/radmeth/radmeth_design.hpp | 79 ++++++++++ 2 files changed, 341 insertions(+) create mode 100644 src/radmeth/radmeth_design.cpp create mode 100644 src/radmeth/radmeth_design.hpp diff --git a/src/radmeth/radmeth_design.cpp b/src/radmeth/radmeth_design.cpp new file mode 100644 index 00000000..f08dd9ea --- /dev/null +++ b/src/radmeth/radmeth_design.cpp @@ -0,0 +1,262 @@ +/* Copyright (C) 2025 Andrew D Smith + * + * Authors: Andrew D. Smith + * + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. + * + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + */ + +#include "radmeth_design.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +static void +make_groups(Design &design) { + auto s = design.matrix; + std::sort(std::begin(s), std::end(s)); + s.erase(std::unique(std::begin(s), std::end(s)), std::end(s)); + design.groups = std::move(s); +} + +static void +assign_groups(Design &design) { + const auto &matrix = design.matrix; + const auto &s = design.groups; + const auto n_samples = design.n_samples(); + auto &group_id = design.group_id; + group_id.resize(n_samples); + for (auto i = 0u; i < n_samples; ++i) { + const auto x = std::find(std::cbegin(s), std::cend(s), matrix[i]); + group_id[i] = std::distance(std::cbegin(s), x); + } +} + +template +static void +transpose(const std::vector> &mat, + std::vector> &tmat) { + const auto n_row = std::size(mat); + const auto n_col = std::size(mat.front()); + tmat.resize(n_col, std::vector(n_row, 0.0)); + for (auto row_idx = 0u; row_idx < n_row; ++row_idx) + for (auto col_idx = 0u; col_idx < n_col; ++col_idx) + tmat[col_idx][row_idx] = mat[row_idx][col_idx]; +} + +std::istream & +operator>>(std::istream &is, Design &design) { + std::string header_encoding; + std::getline(is, header_encoding); + + std::istringstream header_is(header_encoding); + std::string header_name; + while (header_is >> header_name) + design.factor_names.push_back(header_name); + + std::string row; + while (std::getline(is, row)) { + if (row.empty()) + continue; + + std::istringstream row_is(row); + std::string token; + row_is >> token; + design.sample_names.push_back(token); + + std::vector matrix_row; + while (row_is >> token) { + if (std::size(token) != 1 || (token != "0" && token != "1")) + throw std::runtime_error("Must use binary factor levels:\n" + row); + matrix_row.push_back(token == "1"); + } + + if (std::size(matrix_row) != design.n_factors()) + throw std::runtime_error( + "each row must have as many columns as factors:\n" + row); + + design.matrix.push_back(matrix_row); + } + + transpose(design.matrix, design.tmatrix); + make_groups(design); + assign_groups(design); + + return is; +} + +[[nodiscard]] Design +Design::drop_factor(const std::size_t factor_idx) { + // clang-format off + Design design{ + factor_names, + sample_names, + matrix, + tmatrix, + groups, + group_id, + }; + // clang-format on + design.factor_names.erase(std::begin(design.factor_names) + factor_idx); + for (auto i = 0u; i < n_samples(); ++i) + design.matrix[i].erase(std::begin(design.matrix[i]) + factor_idx); + transpose(design.matrix, design.tmatrix); + make_groups(design); + assign_groups(design); + 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; + const auto n_samples = design.n_samples(); + const auto n_factors = design.n_factors(); + if (n_samples <= max_samples_to_report) { + for (std::size_t factor = 0; factor < n_factors; ++factor) { + if (factor != 0) + out << '\t'; + out << design.factor_names[factor]; + } + out << '\n'; + + for (std::size_t i = 0; i < n_samples; ++i) { + out << design.sample_names[i]; + for (std::size_t j = 0; j < n_factors; ++j) + out << '\t' << static_cast(design.matrix[i][j]); + out << '\n'; + } + } + + // compute the number of samples per group + const auto n_groups = design.n_groups(); + std::vector n_samples_per_group(n_groups, 0); + for (auto s_idx = 0u; s_idx < n_samples; ++s_idx) + ++n_samples_per_group[design.group_id[s_idx]]; + + out << "group_id | factor_levels (" << n_factors << " factors) | n_samples\n"; + for (std::size_t g_idx = 0; g_idx < n_groups; ++g_idx) { + out << g_idx; + for (std::size_t f_idx = 0; f_idx < n_factors; ++f_idx) + out << '\t' << static_cast(design.groups[g_idx][f_idx]); + out << '\t' << n_samples_per_group[g_idx] << '\n'; + } + + return out; +} + +void +ensure_sample_order(const std::string &table_filename, Design &design) { + 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); + design.order_samples(sample_names); +} + +[[nodiscard]] 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; +} + +[[nodiscard]] bool +consistent_sample_names(const Design &design, const std::string &header) { + std::istringstream iss(header); + auto nm_itr(std::begin(design.sample_names)); + const auto nm_end(std::end(design.sample_names)); + std::string token; + while (iss >> token && nm_itr != nm_end) + if (token != *nm_itr++) + return false; + return true; +} + +[[nodiscard]] Design +Design::read_design(const std::string &design_filename) { + std::ifstream design_file(design_filename); + if (!design_file) + throw std::runtime_error("could not open file: " + design_filename); + Design design; + design_file >> design; + return design; +} + +[[nodiscard]] std::uint32_t +Design::get_test_factor_idx(const std::string &test_factor) const { + const auto itr = + std::find(std::cbegin(factor_names), std::cend(factor_names), test_factor); + if (itr == std::cend(factor_names)) + throw std::runtime_error("factor not part of design: " + test_factor); + return std::distance(std::cbegin(factor_names), itr); +} + +[[nodiscard]] bool +Design::has_two_values(const std::size_t test_factor) const { + const auto &tcol = tmatrix[test_factor]; + for (const auto x : tcol) + if (x != tcol[0]) + return true; + return false; +} diff --git a/src/radmeth/radmeth_design.hpp b/src/radmeth/radmeth_design.hpp new file mode 100644 index 00000000..38d37cc6 --- /dev/null +++ b/src/radmeth/radmeth_design.hpp @@ -0,0 +1,79 @@ +/* Copyright (C) 2025 Andrew D Smith + * + * Author: Andrew D Smith + * + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. + * + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + */ + +#ifndef RADMETH_DESIGN_HPP +#define RADMETH_DESIGN_HPP + +#include +#include +#include +#include +#include + +struct Design { + std::vector factor_names; + std::vector sample_names; + 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]] static Design + read_design(const std::string &design_filename); + + [[nodiscard]] std::size_t + n_factors() const { + return std::size(factor_names); + } + + [[nodiscard]] std::size_t + n_groups() const { + return std::size(groups); + } + + [[nodiscard]] std::size_t + n_samples() const { + return std::size(sample_names); + } + + [[nodiscard]] Design + drop_factor(const std::size_t factor_idx); + + void + order_samples(const std::vector &ordered_names); + + [[nodiscard]] std::uint32_t + get_test_factor_idx(const std::string &test_factor) const; + + [[nodiscard]] bool + has_two_values(const std::size_t test_factor) const; +}; + +std::istream & +operator>>(std::istream &is, Design &design); + +std::ostream & +operator<<(std::ostream &os, const Design &design); + +void +ensure_sample_order(const std::string &table_filename, Design &design); + +[[nodiscard]] std::vector +get_sample_names_from_header(const std::string &header); + +[[nodiscard]] bool +consistent_sample_names(const Design &design, const std::string &header); + +#endif // RADMETH_DESIGN_HPP From 3b3b79bf1fc8c27212913618297bff28a3f2be75 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 15 Oct 2025 16:16:24 -0700 Subject: [PATCH 05/33] ../src/radmeth/radmeth_model_nano.hcpp: adding files for the RegressionNano class for radmeth-nano --- src/radmeth/radmeth_model_nano.cpp | 80 ++++++++++++++++++++++++++ src/radmeth/radmeth_model_nano.hpp | 90 ++++++++++++++++++++++++++++++ 2 files changed, 170 insertions(+) create mode 100644 src/radmeth/radmeth_model_nano.cpp create mode 100644 src/radmeth/radmeth_model_nano.hpp diff --git a/src/radmeth/radmeth_model_nano.cpp b/src/radmeth/radmeth_model_nano.cpp new file mode 100644 index 00000000..061afe80 --- /dev/null +++ b/src/radmeth/radmeth_model_nano.cpp @@ -0,0 +1,80 @@ +/* Copyright (C) 2025 Andrew D Smith + * + * Authors: Andrew D. Smith + * + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. + * + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + */ + +#include "radmeth_model_nano.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +double RegressionNano::tolerance = 1e-4; +double RegressionNano::stepsize = 0.01; +std::uint32_t RegressionNano::max_iter = 250; + +void +RegressionNano::parse(const std::string &line) { + const auto first_ws = line.find_first_of(" \t"); + + // Parse the row name (must be like: "chr:position:strand:context") + bool failed = (first_ws == std::string::npos); + + auto field_s = line.data(); + auto field_e = line.data() + first_ws; + rowname = std::string{field_s, field_e}; + std::replace(std::begin(rowname), std::end(rowname), ':', '\t'); + if (failed) + throw std::runtime_error("failed to parse label from:\n" + line); + + // Parse the counts of total reads and methylated reads + const auto is_sep = [](const auto x) { return std::isspace(x); }; + const auto not_sep = [](const auto x) { return std::isdigit(x); }; + + const auto line_end = line.data() + std::size(line); + mcountsf mc1{}; + mc.clear(); + while (field_e != line_end) { + // get the total count + field_s = std::find_if(field_e + 1, line_end, not_sep); + field_e = std::find_if(field_s + 1, line_end, is_sep); + { + const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_reads); + failed = failed || (ec != std::errc{}); + } + + // get the methylated count + field_s = std::find_if(field_e + 1, line_end, not_sep); + field_e = std::find_if(field_s + 1, line_end, is_sep); + { +#ifdef __APPLE__ + const int ret = std::sscanf(field_s, "%lf", &mc1.n_meth); + failed = failed || (ret < 1); +#else + const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_meth); + failed = failed || (ec != std::errc{}); +#endif + } + + mc.push_back(mc1); + } + + if (failed) + throw std::runtime_error("failed to parse counts from:\n" + line); +} diff --git a/src/radmeth/radmeth_model_nano.hpp b/src/radmeth/radmeth_model_nano.hpp new file mode 100644 index 00000000..9e3c71ea --- /dev/null +++ b/src/radmeth/radmeth_model_nano.hpp @@ -0,0 +1,90 @@ +/* Copyright (C) 2025 Andrew D Smith + * + * Author: Andrew D Smith + * + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. + * + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + */ + +#ifndef RADMETH_MODEL_NANO_HPP +#define RADMETH_MODEL_NANO_HPP + +#include "radmeth_design.hpp" + +#include +#include +#include +#include +#include + +struct mcountsf { + std::uint32_t n_reads{}; + double n_meth{}; +}; + +[[nodiscard]] inline std::istream & +operator>>(std::istream &in, mcountsf &rm) { + return in >> rm.n_reads >> rm.n_meth; +} + +struct vars_cache { + double p{}; + double a{}; + double b{}; + double lgamma_a{}; + double lgamma_b{}; + double lgamma_a_b{}; + double digamma_a{}; + double digamma_b{}; + double digamma_a_b{}; +}; + +struct RegressionNano { + static double tolerance; // 1e-3; + static double stepsize; // 0.001; + static std::uint32_t max_iter; // 250; + + Design design; + std::string rowname; + std::vector mc; + double max_loglik{}; + + std::vector cache; // scratch space + + void + parse(const std::string &line); + + [[nodiscard]] std::size_t + n_factors() const { + return design.n_factors(); + } + + [[nodiscard]] std::size_t + n_groups() const { + return design.n_groups(); + } + + [[nodiscard]] std::size_t + n_params() const { + return n_factors() + 1; + } + + [[nodiscard]] std::size_t + props_size() const { + return std::size(mc); + } + + [[nodiscard]] std::size_t + n_samples() const { + return design.n_samples(); + } +}; + +#endif // RADMETH_MODEL_NANO_HPP From da61e8da09fa9cee5373fb132f9dab6ef1c4b812 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 15 Oct 2025 16:17:05 -0700 Subject: [PATCH 06/33] ../src/radmeth/radmeth_optimize_nano.hcpp: adding optimization code in a separate file for optimizing paramters of the model for nano with fractional observations --- src/radmeth/radmeth_optimize_nano.cpp | 330 ++++++++++++++++++++++++++ src/radmeth/radmeth_optimize_nano.hpp | 27 +++ 2 files changed, 357 insertions(+) create mode 100644 src/radmeth/radmeth_optimize_nano.cpp create mode 100644 src/radmeth/radmeth_optimize_nano.hpp diff --git a/src/radmeth/radmeth_optimize_nano.cpp b/src/radmeth/radmeth_optimize_nano.cpp new file mode 100644 index 00000000..25024337 --- /dev/null +++ b/src/radmeth/radmeth_optimize_nano.cpp @@ -0,0 +1,330 @@ +/* Copyright (C) 2013-2025 Andrew D Smith + * + * Author: Andrew D. Smith + * + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. + * + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + */ + +#include "radmeth_optimize_nano.hpp" +#include "radmeth_model_nano.hpp" + +#include +#include // for gsl_sf_psi (digamma) +#include + +#include +#include +#include +#include +#include + +/* Coefficients for the Chebyschev polynomial for the digamma function in the + range 0-1 */ +// clang-format off +static constexpr std::array psi1_cs { + -0.038057080835217922, // == -0.019028540417608961*2 + 0.491415393029387130, + -0.056815747821244730, + 0.008357821225914313, + -0.001333232857994342, + 0.000220313287069308, + -0.000037040238178456, + 0.000006283793654854, + -0.000001071263908506, + 0.000000183128394654, + -0.000000031353509361, + 0.000000005372808776, + -0.000000000921168141, + 0.000000000157981265, + -0.000000000027098646, + 0.000000000004648722, + -0.000000000000797527, + 0.000000000000136827, + -0.000000000000023475, + 0.000000000000004027, + -0.000000000000000691, + 0.000000000000000118, + -0.000000000000000020 +}; +// clang-format on + +/* Alternate set of coefficients for the Chebyschev polynomial for the digamma + function */ +// clang-format off +static constexpr std::array psi2_cs = { + -0.0204749044678185, // == -0.01023745223390925*2 + -0.0101801271534859, + 0.0000559718725387, + -0.0000012917176570, + 0.0000000572858606, + -0.0000000038213539, + 0.0000000003397434, + -0.0000000000374838, + 0.0000000000048990, + -0.0000000000007344, + 0.0000000000001233, + -0.0000000000000228, + 0.0000000000000045, + -0.0000000000000009, + 0.0000000000000002, + -0.0000000000000000 , +}; +// clang-format on + +template +[[nodiscard]] static inline double +chebyschev(const std::array &coeffs, std::uint32_t order, + const double y) { + const auto y2 = 2.0 * y; + double d = 0.0; + double dd = 0.0; + for (auto j = order; j >= 1; j--) { + const auto temp = d; + d = y2 * d - dd + coeffs[j]; + dd = temp; + } + return y * d - dd + 0.5 * coeffs[0]; +} + +[[nodiscard]] static inline double +logistic(const double x) { + return 1.0 / (1.0 / std::exp(x) + 1.0); +} + +template +[[nodiscard]] static double +get_p(const std::vector &v, const gsl_vector *params) { + const auto a = v.data(); + return logistic(std::inner_product(a, a + std::size(v), params->data, 0.0)); +} + +static void +get_cache_lgamma(const std::vector &group, + const gsl_vector *params, const double phi, + vars_cache &cache) { + const auto p = get_p(group, params); + const auto a = p * phi; + const auto b = (1.0 - p) * phi; + cache.p = p; + cache.a = a; + cache.b = b; + cache.lgamma_a = std::lgamma(a); + cache.lgamma_b = std::lgamma(b); + cache.lgamma_a_b = std::lgamma(a + b); +} + +[[nodiscard]] static double +log_likelihood_nano(const gsl_vector *params, RegressionNano ®) { + const auto phi = 1.0 / std::exp(gsl_vector_get(params, reg.n_factors())); + + const auto n_groups = reg.n_groups(); + const auto &groups = reg.design.groups; + auto &cache = reg.cache; + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + get_cache_lgamma(groups[g_idx], params, phi, cache[g_idx]); + + const auto &group_id = reg.design.group_id; + const auto &mc = reg.mc; + double ll = 0.0; + const auto n_samples = reg.n_samples(); + for (auto i = 0u; i < n_samples; ++i) { + const auto y = mc[i].n_meth; + const auto n = mc[i].n_reads; + + const auto &c = cache[group_id[i]]; + const auto a = c.a; + const auto b = c.b; + + // clang-format off + ll += ((std::lgamma(y + a) - c.lgamma_a) + + (std::lgamma(n - y + b) - c.lgamma_b) + + (c.lgamma_a_b - std::lgamma(n + a + b))); + // clang-format on + } + + return ll; +} + +// digamma for x non-negative +static double +digamma(const double y) { + static constexpr auto psi_order = 7; // max=22; + static constexpr auto apsi_order = 7; // max=15; + if (y >= 2.0) { + const auto t = 8.0 / (y * y) - 1.0; + return std::log(y) - 0.5 / y + chebyschev(psi2_cs, apsi_order, t); + } + if (y < 1.0) + return -1.0 / y + chebyschev(psi1_cs, psi_order, 2.0 * y - 1.0); + return chebyschev(psi1_cs, psi_order, 2.0 * (y - 1.0) - 1.0); +} + +static void +get_cache_digamma(const std::vector &group, + const gsl_vector *params, const double phi, + vars_cache &cache) { + const auto p = get_p(group, params); + const auto a = p * phi; + const auto b = (1.0 - p) * phi; + cache.p = p; + cache.a = a; + cache.b = b; + cache.digamma_a = digamma(a); + cache.digamma_b = digamma(b); + cache.digamma_a_b = digamma(a + b); +} + +static void +gradient_nano(const gsl_vector *params, RegressionNano ®, + gsl_vector *output) { + const auto n_factors = reg.n_factors(); + const auto phi = 1.0 / std::exp(gsl_vector_get(params, n_factors)); + + const auto &mc = reg.mc; + const auto &matrix = reg.design.matrix; + + const auto n_groups = reg.n_groups(); + const auto &groups = reg.design.groups; + auto &cache = reg.cache; // ADS: reusing scratch space + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + get_cache_digamma(groups[g_idx], params, phi, cache[g_idx]); + + gsl_vector_set_all(output, 0.0); // init output to zero for all factors + auto &data = output->data; + + const auto &group_id = reg.design.group_id; + double grad_phi = 0.0; + const auto n_samples = reg.n_samples(); + for (auto i = 0u; i < n_samples; ++i) { + const auto y = mc[i].n_meth; + const auto n = mc[i].n_reads; + + const auto &c = cache[group_id[i]]; + const auto p = c.p; + const auto a = c.a; + const auto b = c.b; + + const auto digamma_a = c.digamma_a; + const auto digamma_b = c.digamma_b; + const auto digamma_y_a = digamma(y + a); + const auto digamma_n_y_b = digamma(n - y + b); + + // grad wrt p + const auto dlogl_dp = + phi * (digamma_y_a - digamma_a - digamma_n_y_b + digamma_b); + const auto dp_delta = p * (1.0 - p); // chain rule: d p / d param + const auto dlogl_delta = dlogl_dp * dp_delta; // grad wrt params + + auto matrix_itr = std::cbegin(matrix[i]); + const auto data_end = data + n_factors; + for (auto data_itr = data; data_itr != data_end; ++data_itr) + *data_itr += (*matrix_itr++) * dlogl_delta; + + const auto digamma_delta = c.digamma_a_b - digamma(n + a + b); + const auto dphi_term1 = p * (digamma_y_a - digamma_a + digamma_delta); + const auto dphi_term2 = + (1.0 - p) * (digamma_n_y_b - digamma_b + digamma_delta); + + grad_phi += dphi_term1 + dphi_term2; + } + + const auto grad_theta = -grad_phi * phi; + gsl_vector_set(output, n_factors, grad_theta); +} + +[[nodiscard]] static double +neg_loglik_nano(const gsl_vector *params, void *object) { + auto reg = static_cast(object); + return -log_likelihood_nano(params, *reg); +} + +static void +neg_gradient_nano(const gsl_vector *params, void *object, gsl_vector *output) { + auto reg = static_cast(object); + gradient_nano(params, *reg, output); + gsl_vector_scale(output, -1.0); +} + +static void +neg_loglik_and_grad_nano(const gsl_vector *params, void *object, + double *loglik_val, gsl_vector *d_loglik_val) { + *loglik_val = neg_loglik_nano(params, object); + neg_gradient_nano(params, object, d_loglik_val); +} + +void +fit_regression_model_nano(RegressionNano &r, std::vector &p_estimates, + double &dispersion_estimate) { + static constexpr auto init_dispersion_param = -2.5; + const auto stepsize = RegressionNano::stepsize; + const auto max_iter = RegressionNano::max_iter; + const auto n_groups = r.n_groups(); + r.cache.resize(n_groups); // make sure scratch space is allocated + + const std::size_t n_params = r.n_params(); + const auto tol = + std::sqrt(n_params) * r.n_samples() * RegressionNano::tolerance; + // clang-format off + auto loglik_bundle = gsl_multimin_function_fdf{ + &neg_loglik_nano, // objective function + &neg_gradient_nano, // gradient + &neg_loglik_and_grad_nano, // combined obj and grad + n_params, // number of model params + static_cast(&r) // parameters for objective and gradient functions + }; + // clang-format on + + // set the parameters: zero for "p" parameters and the final one for + // dispersion using the constant + auto params = gsl_vector_alloc(n_params); + gsl_vector_set_all(params, 0.0); + gsl_vector_set(params, n_params - 1, init_dispersion_param); + + // Alternatives: + // - gsl_multimin_fdfminimizer_conjugate_pr + // - gsl_multimin_fdfminimizer_conjugate_fr + // - gsl_multimin_fdfminimizer_vector_bfgs2 + // - gsl_multimin_fdfminimizer_steepest_descent + const auto minimizer = gsl_multimin_fdfminimizer_conjugate_pr; + auto s = gsl_multimin_fdfminimizer_alloc(minimizer, n_params); + + gsl_multimin_fdfminimizer_set(s, &loglik_bundle, params, stepsize, tol); + + int status = 0; + std::size_t iter = 0; + do { + status = gsl_multimin_fdfminimizer_iterate(s); // one iter and get status + if (status) + break; + // check status from gradient + status = gsl_multimin_test_gradient(s->gradient, tol); + } while (status == GSL_CONTINUE && ++iter < max_iter); + + /// ADS: the condition below might not always pass even if we are doing + /// ok. It's not clear how to check for failure. + + // if (status != GSL_SUCCESS) + // throw std::runtime_error("failed to fit model parameters"); + + const auto param_estimates = gsl_multimin_fdfminimizer_x(s); + + const auto &groups = r.design.groups; + p_estimates.clear(); + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + p_estimates.push_back(get_p(groups[g_idx], param_estimates)); + const auto disp_param = gsl_vector_get(param_estimates, n_params - 1); + dispersion_estimate = 1.0 / std::exp(disp_param); + + r.max_loglik = log_likelihood_nano(param_estimates, r); + + gsl_multimin_fdfminimizer_free(s); + gsl_vector_free(params); +} diff --git a/src/radmeth/radmeth_optimize_nano.hpp b/src/radmeth/radmeth_optimize_nano.hpp new file mode 100644 index 00000000..3487cbea --- /dev/null +++ b/src/radmeth/radmeth_optimize_nano.hpp @@ -0,0 +1,27 @@ +/* Copyright (C) 2025 Andrew D Smith + * + * Author: Andrew D. Smith + * + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. + * + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + */ + +#ifndef RADMETH_OPTIMIZE_NANO_HPP +#define RADMETH_OPTIMIZE_NANO_HPP + +#include + +struct RegressionNano; + +void +fit_regression_model_nano(RegressionNano &r, std::vector &p_estimates, + double &dispersion_estimate); + +#endif // RADMETH_OPTIMIZE_NANO_HPP From 45af7a8584628ae3dc5b54889aa56bc06f4ed494 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 15 Oct 2025 16:17:38 -0700 Subject: [PATCH 07/33] src/radmeth/radmeth_utils.hcpp: adding files for common functions shared by radmeth and radmeth-nano --- src/radmeth/radmeth_utils.cpp | 144 ++++++++++++++++++++++++++++++++++ src/radmeth/radmeth_utils.hpp | 43 ++++++++++ 2 files changed, 187 insertions(+) create mode 100644 src/radmeth/radmeth_utils.cpp create mode 100644 src/radmeth/radmeth_utils.hpp diff --git a/src/radmeth/radmeth_utils.cpp b/src/radmeth/radmeth_utils.cpp new file mode 100644 index 00000000..54a39c1c --- /dev/null +++ b/src/radmeth/radmeth_utils.cpp @@ -0,0 +1,144 @@ +/* Copyright (C) 2025 Andrew D Smith + * + * Author: Andrew D. Smith + * + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. + * + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + */ + +#include "radmeth_utils.hpp" + +#include +#include +#include +#include +#include +#include +#include + +[[nodiscard]] std::string +format_duration(const std::chrono::duration elapsed) { + static constexpr auto s_per_h = 3600; + static constexpr auto s_per_m = 60; + const double tot_s = elapsed.count(); + + // break down into hours, minutes, seconds + const std::uint32_t hours = tot_s / 3600; + const std::uint32_t minutes = (static_cast(tot_s) % s_per_h) / s_per_m; + const double seconds = tot_s - (hours * s_per_h) - (minutes * s_per_m); + + std::ostringstream oss; + oss << std::setfill('0') << std::setw(2) << hours << ":" << std::setfill('0') + << std::setw(2) << minutes << ":" << std::fixed << std::setprecision(2) + << std::setw(5) << seconds; + return oss.str(); +} + +file_progress::file_progress(const std::string &filename) : + one_thousand_over_filesize{1000.0 / std::filesystem::file_size(filename)} {} + +void +file_progress::operator()(std::ifstream &in) { + const std::size_t curr_offset = + in.eof() ? 1000 : in.tellg() * one_thousand_over_filesize; + if (curr_offset <= prev_offset) + return; + std::ios old_state(nullptr); + old_state.copyfmt(std::cerr); + std::cerr << "\r[progress: " << std::setw(5) << std::fixed + << std::setprecision(1) << (curr_offset / 10.0) + << (curr_offset == 1000 ? "%]\n" : "%]"); + std::cerr.copyfmt(old_state); + prev_offset = (curr_offset == 1000) ? std::numeric_limits::max() + : curr_offset; +} + +// 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. +[[nodiscard]] 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); + + // It is assumed that null model has one fewer factor than the full + // model. Hence the number of degrees of freedom is 1. + const std::size_t degrees_of_freedom = 1; + + // Log-likelihood ratio statistic has a chi-sqare distribution. + const double chisq_p = chi_square_cdf(log_lik_stat, degrees_of_freedom); + + const double p_value = 1.0 - chisq_p; + + return p_value; +} diff --git a/src/radmeth/radmeth_utils.hpp b/src/radmeth/radmeth_utils.hpp new file mode 100644 index 00000000..cae962c1 --- /dev/null +++ b/src/radmeth/radmeth_utils.hpp @@ -0,0 +1,43 @@ +/* Copyright (C) 2025 Andrew D Smith + * + * Author: Andrew D. Smith + * + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. + * + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + */ + +#ifndef RADMETH_UTILS_HPP +#define RADMETH_UTILS_HPP + +#include +#include +#include +#include + +[[nodiscard]] std::string +format_duration(const std::chrono::duration elapsed); + +struct file_progress { + double one_thousand_over_filesize{}; + std::size_t prev_offset{}; + explicit file_progress(const std::string &filename); + void + operator()(std::ifstream &in); +}; + +[[nodiscard]] double +llr_test(const double null_loglik, const double full_loglik); + +[[nodiscard]] inline double +overdispersion_factor(const std::uint32_t n_samples, const double dispersion) { + return (n_samples - 1) / (dispersion + 1); +} + +#endif // RADMETH_UTILS_HPP From 392ab7f38ce4b52bcafe921352b49276fe2f1d89 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 15 Oct 2025 16:18:11 -0700 Subject: [PATCH 08/33] src/radmeth/radmeth.cpp: factoring out code shared by radmeth and radmeth-nano --- src/radmeth/radmeth.cpp | 257 ++++++---------------------------------- 1 file changed, 34 insertions(+), 223 deletions(-) diff --git a/src/radmeth/radmeth.cpp b/src/radmeth/radmeth.cpp index fcf9bae1..2a8f215d 100644 --- a/src/radmeth/radmeth.cpp +++ b/src/radmeth/radmeth.cpp @@ -14,8 +14,10 @@ * more details. */ +#include "radmeth_design.hpp" #include "radmeth_model.hpp" #include "radmeth_optimize.hpp" +#include "radmeth_utils.hpp" // smithlab headers #include "GenomicRegion.hpp" @@ -26,7 +28,6 @@ #include #include #include -#include #include #include #include @@ -36,169 +37,25 @@ #include #include -[[nodiscard]] static std::string -format_duration(const std::chrono::duration elapsed) { - static constexpr auto s_per_h = 3600; - static constexpr auto s_per_m = 60; - const double tot_s = elapsed.count(); - - // break down into hours, minutes, seconds - const std::uint32_t hours = tot_s / 3600; - const std::uint32_t minutes = (static_cast(tot_s) % s_per_h) / s_per_m; - const double seconds = tot_s - (hours * s_per_h) - (minutes * s_per_m); - - std::ostringstream oss; - oss << std::setfill('0') << std::setw(2) << hours << ":" << std::setfill('0') - << std::setw(2) << minutes << ":" << std::fixed << std::setprecision(2) - << std::setw(5) << seconds; - return oss.str(); -} - -struct file_progress { - double one_thousand_over_filesize{}; - std::size_t prev_offset{}; - explicit file_progress(const std::string &filename) : - one_thousand_over_filesize{1000.0 / std::filesystem::file_size(filename)} {} - void - operator()(std::ifstream &in) { - const std::size_t curr_offset = - in.eof() ? 1000 : in.tellg() * one_thousand_over_filesize; - if (curr_offset <= prev_offset) - return; - std::ios old_state(nullptr); - old_state.copyfmt(std::cerr); - std::cerr << "\r[progress: " << std::setw(5) << std::fixed - << std::setprecision(1) << (curr_offset / 10.0) - << (curr_offset == 1000 ? "%]\n" : "%]"); - std::cerr.copyfmt(old_state); - prev_offset = (curr_offset == 1000) - ? std::numeric_limits::max() - : curr_offset; - } -}; - +template [[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)); - const auto nm_end(std::end(reg.design.sample_names)); - std::string token; - while (iss >> token && nm_itr != nm_end) - if (token != *nm_itr++) - return false; - 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. -[[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); - - // It is assumed that null model has one fewer factor than the full - // model. Hence the number of degrees of freedom is 1. - const std::size_t degrees_of_freedom = 1; - - // Log-likelihood ratio statistic has a chi-sqare distribution. - const double chisq_p = chi_square_cdf(log_lik_stat, degrees_of_freedom); - - const double p_value = 1.0 - chisq_p; - - return p_value; -} - -static bool -has_low_coverage(const Regression ®, const std::size_t test_factor) { +has_low_coverage(const RegressionType ®, const std::size_t test_factor) { bool cvrd_in_test_fact_smpls = false; const auto &tcol = reg.design.tmatrix[test_factor]; for (std::size_t i = 0; i < reg.n_samples() && !cvrd_in_test_fact_smpls; ++i) - cvrd_in_test_fact_smpls = (tcol[i] == 1 && reg.props.mc[i].n_reads != 0); + cvrd_in_test_fact_smpls = (tcol[i] == 1 && reg.mc[i].n_reads != 0); bool cvrd_in_other_smpls = false; for (std::size_t i = 0; i < reg.n_samples() && !cvrd_in_other_smpls; ++i) - cvrd_in_other_smpls = (tcol[i] != 1 && reg.props.mc[i].n_reads != 0); + cvrd_in_other_smpls = (tcol[i] != 1 && reg.mc[i].n_reads != 0); return !cvrd_in_test_fact_smpls || !cvrd_in_other_smpls; } +template [[nodiscard]] static bool -has_extreme_counts(const Regression ®) { - const auto &mc = reg.props.mc; +has_extreme_counts(const RegressionType ®) { + const auto &mc = reg.mc; bool full_meth = true; for (auto i = std::cbegin(mc); i != std::cend(mc) && full_meth; ++i) @@ -211,42 +68,6 @@ has_extreme_counts(const Regression ®) { return full_meth || no_meth; } -[[nodiscard]] static bool -has_two_values(const Regression ®, const std::size_t test_factor) { - const auto &tcol = reg.design.tmatrix[test_factor]; - for (const auto x : tcol) - if (x != tcol[0]) - return true; - return false; -} - -[[nodiscard]] static std::uint32_t -get_test_factor_idx(const Regression &model, const std::string &test_factor) { - const auto &factors = model.design.factor_names; - const auto itr = - std::find(std::cbegin(factors), std::cend(factors), test_factor); - - if (itr == std::cend(factors)) - throw std::runtime_error("factor not part of design: " + test_factor); - - return std::distance(std::cbegin(factors), itr); -} - -[[nodiscard]] static Design -read_design(const std::string &design_filename) { - std::ifstream design_file(design_filename); - if (!design_file) - throw std::runtime_error("could not open file: " + design_filename); - Design design; - design_file >> design; - return design; -} - -[[nodiscard]] static inline double -overdispersion_factor(const std::uint32_t n_samples, const double dispersion) { - return (n_samples - 1) / (dispersion + 1); -} - enum class row_status : std::uint8_t { ok, na, @@ -254,11 +75,12 @@ enum class row_status : std::uint8_t { na_extreme_cnt, }; +template static void radmeth(const bool show_progress, const bool more_na_info, const std::uint32_t n_threads, const std::string &table_filename, - const std::string &outfile, const Regression &alt_model, - const Regression &null_model, const std::uint32_t test_factor_idx) { + const std::string &outfile, const RegressionType &alt_model, + const RegressionType &null_model, const std::uint32_t test_factor_idx) { static constexpr auto buf_size = 1024; static constexpr auto max_lines = 16384; @@ -274,7 +96,7 @@ radmeth(const bool show_progress, const bool more_na_info, 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)) { + if (!consistent_sample_names(alt_model.design, sample_names_header)) { // clang-format off const auto message = R"( @@ -298,8 +120,8 @@ that the design matrix and the proportion table are correctly formatted. std::vector n_bytes(max_lines, 0); std::vector lines(max_lines); - std::vector alt_models(n_threads, alt_model); - std::vector null_models(n_threads, null_model); + std::vector alt_models(n_threads, alt_model); + std::vector null_models(n_threads, null_model); const auto n_groups = alt_model.n_groups(); @@ -328,7 +150,7 @@ that the design matrix and the proportion table are correctly formatted. // balances work better. if (b % n_threads != thread_id) continue; - t_alt_model.props.parse(lines[b]); + t_alt_model.parse(lines[b]); if (t_alt_model.props_size() != n_samples) throw std::runtime_error("found row with wrong number of columns"); @@ -344,7 +166,8 @@ that the design matrix and the proportion table are correctly formatted. fit_regression_model(t_alt_model, p_estim_alt, phi_estim_alt); - t_null_model.props = t_alt_model.props; + t_null_model.mc = t_alt_model.mc; + t_null_model.rowname = t_alt_model.rowname; fit_regression_model(t_null_model, p_estim_null, phi_estim_null); @@ -361,7 +184,7 @@ that the design matrix and the proportion table are correctly formatted. n_bytes[b] = [&] { // clang-format off const int n_prefix_bytes = - std::sprintf(bufs[b].data(), "%s\t", t_alt_model.rowname().data()); + std::sprintf(bufs[b].data(), "%s\t", t_alt_model.rowname.data()); // clang-format on if (n_prefix_bytes < 0) return n_prefix_bytes; @@ -432,19 +255,6 @@ 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 { @@ -480,7 +290,6 @@ main_radmeth(int argc, char *argv[]) { opt_parse.add_opt("max-iter", '\0', "max iterations when estimating parameters", false, Regression::max_iter); - std::vector leftover_args; opt_parse.parse(argc, argv, leftover_args); if (argc == 1 || opt_parse.help_requested()) { @@ -507,30 +316,26 @@ main_radmeth(int argc, char *argv[]) { if (verbose) std::cerr << "design table filename: " << design_filename << "\n\n"; - // initialize full design matrix from file - Regression alt_model; - alt_model.design = read_design(design_filename); - + Design design = Design::read_design(design_filename); if (verbose) - std::cerr << "Alternate model:\n" << alt_model.design << '\n'; + std::cerr << "Alternate model:\n" << design << '\n'; // Check that provided test factor name exists and find its index. - const auto test_factor_idx = get_test_factor_idx(alt_model, test_factor); + const auto test_factor_idx = design.get_test_factor_idx(test_factor); + + ensure_sample_order(table_filename, design); // verify that the design includes more than one level for the // test factor - if (!has_two_values(alt_model, test_factor_idx)) { - const auto first_level = alt_model.design.matrix[0][test_factor_idx]; + if (!design.has_two_values(test_factor_idx)) { + const auto first_level = design.matrix[0][test_factor_idx]; throw std::runtime_error("test factor only one level: " + test_factor + ", " + std::to_string(first_level)); } - Regression null_model = alt_model; - null_model.design = alt_model.design.drop_factor(test_factor_idx); + const Design null_design = design.drop_factor(test_factor_idx); if (verbose) - std::cerr << "Null model:\n" << null_model.design << '\n'; - - ensure_sample_order(table_filename, alt_model, null_model); + std::cerr << "Null model:\n" << null_design << '\n'; // clang-format off if (verbose) @@ -546,6 +351,12 @@ main_radmeth(int argc, char *argv[]) { << "estimated methylation is for test factor value (0 or 1)\n" << '\n'; // clang-format on + + Regression alt_model; + alt_model.design = design; + Regression null_model; + null_model.design = null_design; + const auto start_time = std::chrono::steady_clock::now(); radmeth(show_progress, more_na_info, n_threads, table_filename, outfile, alt_model, null_model, test_factor_idx); From 99e696ac9a29408bc41ead6c5bfb5ae1ad90d4bb Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 15 Oct 2025 16:18:45 -0700 Subject: [PATCH 09/33] src/radmeth/radmeth_model.hcpp: moving Design-related code to a separate file --- src/radmeth/radmeth_model.cpp | 188 +--------------------------------- src/radmeth/radmeth_model.hpp | 111 ++++++-------------- 2 files changed, 32 insertions(+), 267 deletions(-) diff --git a/src/radmeth/radmeth_model.cpp b/src/radmeth/radmeth_model.cpp index 1aec01fe..db39b8c7 100644 --- a/src/radmeth/radmeth_model.cpp +++ b/src/radmeth/radmeth_model.cpp @@ -20,9 +20,7 @@ #include #include #include -#include #include -#include #include double Regression::tolerance = 1e-4; @@ -30,7 +28,7 @@ double Regression::stepsize = 0.01; std::uint32_t Regression::max_iter = 250; void -SiteProp::parse(const std::string &line) { +Regression::parse(const std::string &line) { const auto first_ws = line.find_first_of(" \t"); // Parse the row name (must be like: "chr:position:strand:context") @@ -56,20 +54,14 @@ SiteProp::parse(const std::string &line) { field_e = std::find_if(field_s + 1, line_end, is_sep); { const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_reads); - failed = failed || (ec != std::errc{}); + failed = failed || (ptr != field_e); } - // get the methylated count field_s = std::find_if(field_e + 1, line_end, not_sep); field_e = std::find_if(field_s + 1, line_end, is_sep); { -#ifdef __APPLE__ - const int ret = std::sscanf(field_s, "%lf", &mc1.n_meth); - failed = failed || (ret < 1); -#else const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_meth); - failed = failed || (ec != std::errc{}); -#endif + failed = failed || (ptr != field_e); } mc.push_back(mc1); @@ -78,177 +70,3 @@ SiteProp::parse(const std::string &line) { if (failed) throw std::runtime_error("failed to parse counts from:\n" + line); } - -static void -make_groups(Design &design) { - auto s = design.matrix; - std::sort(std::begin(s), std::end(s)); - s.erase(std::unique(std::begin(s), std::end(s)), std::end(s)); - design.groups = std::move(s); -} - -static void -assign_groups(Design &design) { - const auto &matrix = design.matrix; - const auto &s = design.groups; - const auto n_samples = design.n_samples(); - auto &group_id = design.group_id; - group_id.resize(n_samples); - for (auto i = 0u; i < n_samples; ++i) { - const auto x = std::find(std::cbegin(s), std::cend(s), matrix[i]); - group_id[i] = std::distance(std::cbegin(s), x); - } -} - -template -static void -transpose(const std::vector> &mat, - std::vector> &tmat) { - const auto n_row = std::size(mat); - const auto n_col = std::size(mat.front()); - tmat.resize(n_col, std::vector(n_row, 0.0)); - for (auto row_idx = 0u; row_idx < n_row; ++row_idx) - for (auto col_idx = 0u; col_idx < n_col; ++col_idx) - tmat[col_idx][row_idx] = mat[row_idx][col_idx]; -} - -std::istream & -operator>>(std::istream &is, Design &design) { - std::string header_encoding; - std::getline(is, header_encoding); - - std::istringstream header_is(header_encoding); - std::string header_name; - while (header_is >> header_name) - design.factor_names.push_back(header_name); - - std::string row; - while (std::getline(is, row)) { - if (row.empty()) - continue; - - std::istringstream row_is(row); - std::string token; - row_is >> token; - design.sample_names.push_back(token); - - std::vector matrix_row; - while (row_is >> token) { - if (std::size(token) != 1 || (token != "0" && token != "1")) - throw std::runtime_error("Must use binary factor levels:\n" + row); - matrix_row.push_back(token == "1"); - } - - if (std::size(matrix_row) != design.n_factors()) - throw std::runtime_error( - "each row must have as many columns as factors:\n" + row); - - design.matrix.push_back(matrix_row); - } - - transpose(design.matrix, design.tmatrix); - make_groups(design); - assign_groups(design); - - return is; -} - -[[nodiscard]] Design -Design::drop_factor(const std::size_t factor_idx) { - // clang-format off - Design design{ - factor_names, - sample_names, - matrix, - tmatrix, - groups, - group_id, - }; - // clang-format on - design.factor_names.erase(std::begin(design.factor_names) + factor_idx); - for (auto i = 0u; i < n_samples(); ++i) - design.matrix[i].erase(std::begin(design.matrix[i]) + factor_idx); - transpose(design.matrix, design.tmatrix); - make_groups(design); - assign_groups(design); - 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; - const auto n_samples = design.n_samples(); - const auto n_factors = design.n_factors(); - if (n_samples <= max_samples_to_report) { - for (std::size_t factor = 0; factor < n_factors; ++factor) { - if (factor != 0) - out << '\t'; - out << design.factor_names[factor]; - } - out << '\n'; - - for (std::size_t i = 0; i < n_samples; ++i) { - out << design.sample_names[i]; - for (std::size_t j = 0; j < n_factors; ++j) - out << '\t' << static_cast(design.matrix[i][j]); - out << '\n'; - } - } - - // compute the number of samples per group - const auto n_groups = design.n_groups(); - std::vector n_samples_per_group(n_groups, 0); - for (auto s_idx = 0u; s_idx < n_samples; ++s_idx) - ++n_samples_per_group[design.group_id[s_idx]]; - - out << "group_id | factor_levels (" << n_factors << " factors) | n_samples\n"; - for (std::size_t g_idx = 0; g_idx < n_groups; ++g_idx) { - out << g_idx; - for (std::size_t f_idx = 0; f_idx < n_factors; ++f_idx) - out << '\t' << static_cast(design.groups[g_idx][f_idx]); - out << '\t' << n_samples_per_group[g_idx] << '\n'; - } - - return out; -} diff --git a/src/radmeth/radmeth_model.hpp b/src/radmeth/radmeth_model.hpp index 013e51f4..0be0099d 100644 --- a/src/radmeth/radmeth_model.hpp +++ b/src/radmeth/radmeth_model.hpp @@ -1,69 +1,32 @@ -/* Copyright (C) 2013-2023 University of Southern California and - * Egor Dolzhenko - * Andrew D Smith - * Guilherme Sena +/* Copyright (C) 2025 Andrew D Smith * - * Authors: Andrew D. Smith and Egor Dolzhenko and Guilherme Sena + * Author: Andrew D Smith * - * This program is free software: you can redistribute it and/or - * modify it under the terms of the GNU General Public License as - * published by the Free Software Foundation, either version 3 of the - * License, or (at your option) any later version. + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. * - * This program is distributed in the hope that it will be useful, but - * WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * General Public License for more details. + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. */ #ifndef RADMETH_MODEL_HPP #define RADMETH_MODEL_HPP +#include "radmeth_design.hpp" + #include #include #include #include #include -struct Design { - std::vector factor_names; - std::vector sample_names; - 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 { - return std::size(factor_names); - } - - [[nodiscard]] std::size_t - n_groups() const { - return std::size(groups); - } - - [[nodiscard]] std::size_t - n_samples() const { - return std::size(sample_names); - } - - [[nodiscard]] Design - drop_factor(const std::size_t factor_idx); - - void - order_samples(const std::vector &ordered_names); -}; - -std::istream & -operator>>(std::istream &is, Design &design); - -std::ostream & -operator<<(std::ostream &os, const Design &design); - struct mcounts { std::uint32_t n_reads{}; - double n_meth{}; + std::uint32_t n_meth{}; }; [[nodiscard]] inline std::istream & @@ -71,24 +34,10 @@ operator>>(std::istream &in, mcounts &rm) { return in >> rm.n_reads >> rm.n_meth; } -struct SiteProp { - std::string rowname; - std::vector mc; - - void - parse(const std::string &line); -}; - -struct vars_cache { - double p{}; - double a{}; - double b{}; - double lgamma_a{}; - double lgamma_b{}; - double lgamma_a_b{}; - double digamma_a{}; - double digamma_b{}; - double digamma_a_b{}; +struct cumul_counts { + std::vector m_counts; + std::vector r_counts; + std::vector d_counts; }; struct Regression { @@ -97,10 +46,18 @@ struct Regression { static std::uint32_t max_iter; // 250; Design design; - SiteProp props; + std::string rowname; + std::vector mc; double max_loglik{}; - std::vector cache; // scratch space + // scratch space + std::vector cumul; + std::vector p_v; + std::vector cache; + std::uint32_t max_r_count{}; + + void + parse(const std::string &line); [[nodiscard]] std::size_t n_factors() const { @@ -119,23 +76,13 @@ struct Regression { [[nodiscard]] std::size_t props_size() const { - return std::size(props.mc); + return std::size(mc); } [[nodiscard]] std::size_t n_samples() const { return design.n_samples(); } - - [[nodiscard]] std::string - rowname() const { - return props.rowname; - } - - void - order_samples(const std::vector &ordered_names) { - design.order_samples(ordered_names); - } }; -#endif +#endif // RADMETH_MODEL_HPP From cde6b33f1a48bf2d78a93d09e6a85e43543a0066 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 15 Oct 2025 16:22:15 -0700 Subject: [PATCH 10/33] src/radmeth/radmeth_optimize.cpp: adding back the code for the discrete optimization, with the newer code moved to the new file radmeth_optimize_nano.hcpp for fractional --- src/radmeth/radmeth_optimize.cpp | 356 ++++++++++++++----------------- 1 file changed, 165 insertions(+), 191 deletions(-) diff --git a/src/radmeth/radmeth_optimize.cpp b/src/radmeth/radmeth_optimize.cpp index d2b90bdf..77f6aaba 100644 --- a/src/radmeth/radmeth_optimize.cpp +++ b/src/radmeth/radmeth_optimize.cpp @@ -1,99 +1,29 @@ -/* Copyright (C) 2013-2025 Andrew D Smith +/* Copyright (C) 2025 Andrew D. Smith * * Author: Andrew D. Smith * - * This program is free software: you can redistribute it and/or - * modify it under the terms of the GNU General Public License as - * published by the Free Software Foundation, either version 3 of the - * License, or (at your option) any later version. + * This program is free software: you can redistribute it and/or modify it + * under the terms the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. * - * This program is distributed in the hope that it will be useful, but - * WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * General Public License for more details. + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. */ #include "radmeth_optimize.hpp" #include "radmeth_model.hpp" #include -#include // for gsl_sf_psi (digamma) #include #include -#include #include #include #include -/* Coefficients for the Chebyschev polynomial for the digamma function in the - range 0-1 */ -// clang-format off -static constexpr std::array psi1_cs { - -0.038057080835217922, // == -0.019028540417608961*2 - 0.491415393029387130, - -0.056815747821244730, - 0.008357821225914313, - -0.001333232857994342, - 0.000220313287069308, - -0.000037040238178456, - 0.000006283793654854, - -0.000001071263908506, - 0.000000183128394654, - -0.000000031353509361, - 0.000000005372808776, - -0.000000000921168141, - 0.000000000157981265, - -0.000000000027098646, - 0.000000000004648722, - -0.000000000000797527, - 0.000000000000136827, - -0.000000000000023475, - 0.000000000000004027, - -0.000000000000000691, - 0.000000000000000118, - -0.000000000000000020 -}; -// clang-format on - -/* Alternate set of coefficients for the Chebyschev polynomial for the digamma - function */ -// clang-format off -static constexpr std::array psi2_cs = { - -0.0204749044678185, // == -0.01023745223390925*2 - -0.0101801271534859, - 0.0000559718725387, - -0.0000012917176570, - 0.0000000572858606, - -0.0000000038213539, - 0.0000000003397434, - -0.0000000000374838, - 0.0000000000048990, - -0.0000000000007344, - 0.0000000000001233, - -0.0000000000000228, - 0.0000000000000045, - -0.0000000000000009, - 0.0000000000000002, - -0.0000000000000000 , -}; -// clang-format on - -template -[[nodiscard]] static inline double -chebyschev(const std::array &coeffs, std::uint32_t order, - const double y) { - const auto y2 = 2.0 * y; - double d = 0.0; - double dd = 0.0; - for (auto j = order; j >= 1; j--) { - const auto temp = d; - d = y2 * d - dd + coeffs[j]; - dd = temp; - } - return y * d - dd + 0.5 * coeffs[0]; -} - [[nodiscard]] static inline double logistic(const double x) { return 1.0 / (1.0 / std::exp(x) + 1.0); @@ -106,137 +36,132 @@ get_p(const std::vector &v, const gsl_vector *params) { return logistic(std::inner_product(a, a + std::size(v), params->data, 0.0)); } -static void -get_cache_lgamma(const std::vector &group, - const gsl_vector *params, const double phi, - vars_cache &cache) { - const auto p = get_p(group, params); - const auto a = p * phi; - const auto b = (1.0 - p) * phi; - cache.p = p; - cache.a = a; - cache.b = b; - cache.lgamma_a = std::lgamma(a); - cache.lgamma_b = std::lgamma(b); - cache.lgamma_a_b = std::lgamma(a + b); +static inline auto +set_max_r_count(Regression ®) { + const auto &cumul = reg.cumul; + const auto max_itr = std::max_element( + std::cbegin(cumul), std::cend(cumul), [](const auto &a, const auto &b) { + return std::size(a.r_counts) < std::size(b.r_counts); + }); + reg.max_r_count = std::size(max_itr->r_counts); + // ADS: avoid the realloc that can happen even for resize(smaller_size) + if (reg.max_r_count > std::size(reg.cache)) + reg.cache.resize(reg.max_r_count); +} + +static inline auto +cache_log1p_factors(Regression ®, const double phi) { + const std::size_t max_k = reg.max_r_count; + auto &cache = reg.cache; + for (std::size_t k = 0; k < max_k; ++k) + cache[k] = std::log1p(phi * (k - 1.0)); +} + +static inline auto +cache_dispersion_effect(Regression ®, const double phi) { + const std::size_t max_k = reg.max_r_count; + auto &cache = reg.cache; + for (std::size_t k = 0; k < max_k; ++k) + cache[k] = (k - 1.0) / (1.0 + phi * (k - 1.0)); } [[nodiscard]] static double log_likelihood(const gsl_vector *params, Regression ®) { - const auto phi = 1.0 / std::exp(gsl_vector_get(params, reg.n_factors())); + const auto phi = logistic(gsl_vector_get(params, reg.design.n_factors())); + const auto one_minus_phi = 1.0 - phi; const auto n_groups = reg.n_groups(); const auto &groups = reg.design.groups; - auto &cache = reg.cache; - for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) - get_cache_lgamma(groups[g_idx], params, phi, cache[g_idx]); - - const auto &group_id = reg.design.group_id; - const auto &mc = reg.props.mc; - double ll = 0.0; - const auto n_samples = reg.n_samples(); - for (auto i = 0u; i < n_samples; ++i) { - const auto y = mc[i].n_meth; - const auto n = mc[i].n_reads; - - const auto &c = cache[group_id[i]]; - const auto a = c.a; - const auto b = c.b; - - // clang-format off - ll += ((std::lgamma(y + a) - c.lgamma_a) + - (std::lgamma(n - y + b) - c.lgamma_b) + - (c.lgamma_a_b - std::lgamma(n + a + b))); - // clang-format on - } - - return ll; -} - -// digamma for x non-negative -static double -digamma(const double y) { - static constexpr auto psi_order = 7; // max=22; - static constexpr auto apsi_order = 7; // max=15; - if (y >= 2.0) { - const auto t = 8.0 / (y * y) - 1.0; - return std::log(y) - 0.5 / y + chebyschev(psi2_cs, apsi_order, t); + const auto &cumul = reg.cumul; + + // ADS: precompute the log1p(phi * (k - 1.0)) values, which are reused for + // each group. + cache_log1p_factors(reg, phi); + const auto &log1p_fact_v = reg.cache; + + double log_lik = 0.0; + for (std::size_t g_idx = 0; g_idx < n_groups; ++g_idx) { + const auto p = get_p(groups[g_idx], params); + const auto one_minus_p = 1.0 - p; + + const auto term1 = one_minus_phi * p; + const auto &cumul_y = cumul[g_idx].m_counts; + for (std::size_t k = 0; k < std::size(cumul_y); ++k) + log_lik += cumul_y[k] * std::log(term1 + phi * k); + + const auto term2 = one_minus_phi * one_minus_p; + const auto &cumul_d = cumul[g_idx].d_counts; + for (std::size_t k = 0; k < std::size(cumul_d); ++k) + log_lik += cumul_d[k] * std::log(term2 + phi * k); + + const auto &cumul_n = cumul[g_idx].r_counts; + for (std::size_t k = 0; k < std::size(cumul_n); ++k) + log_lik -= cumul_n[k] * log1p_fact_v[k]; } - if (y < 1.0) - return -1.0 / y + chebyschev(psi1_cs, psi_order, 2.0 * y - 1.0); - return chebyschev(psi1_cs, psi_order, 2.0 * (y - 1.0) - 1.0); -} - -static void -get_cache_digamma(const std::vector &group, - const gsl_vector *params, const double phi, - vars_cache &cache) { - const auto p = get_p(group, params); - const auto a = p * phi; - const auto b = (1.0 - p) * phi; - cache.p = p; - cache.a = a; - cache.b = b; - cache.digamma_a = digamma(a); - cache.digamma_b = digamma(b); - cache.digamma_a_b = digamma(a + b); + return log_lik; } static void gradient(const gsl_vector *params, Regression ®, gsl_vector *output) { - const auto n_factors = reg.n_factors(); - const auto phi = 1.0 / std::exp(gsl_vector_get(params, n_factors)); - - const auto &mc = reg.props.mc; - const auto &matrix = reg.design.matrix; + const auto n_factors = reg.design.n_factors(); + const auto phi = logistic(gsl_vector_get(params, n_factors)); + const auto one_minus_phi = 1.0 - phi; const auto n_groups = reg.n_groups(); const auto &groups = reg.design.groups; - auto &cache = reg.cache; // ADS: reusing scratch space + const auto &cumul = reg.cumul; + + auto &p_v = reg.p_v; // ADS: reusing scratch space for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) - get_cache_digamma(groups[g_idx], params, phi, cache[g_idx]); + p_v[g_idx] = get_p(groups[g_idx], params); + + cache_dispersion_effect(reg, phi); + const auto &dispersion_effect = reg.cache; // (k-1)/(1 + phi(k-1)) - gsl_vector_set_all(output, 0.0); // init output to zero for all factors + // init output to zero for all factors + gsl_vector_set_all(output, 0.0); auto &data = output->data; - const auto &group_id = reg.design.group_id; - double grad_phi = 0.0; - const auto n_samples = reg.n_samples(); - for (auto i = 0u; i < n_samples; ++i) { - const auto y = mc[i].n_meth; - const auto n = mc[i].n_reads; - - const auto &c = cache[group_id[i]]; - const auto p = c.p; - const auto a = c.a; - const auto b = c.b; - - const auto digamma_a = c.digamma_a; - const auto digamma_b = c.digamma_b; - const auto digamma_y_a = digamma(y + a); - const auto digamma_n_y_b = digamma(n - y + b); - - // grad wrt p - const auto dlogl_dp = - phi * (digamma_y_a - digamma_a - digamma_n_y_b + digamma_b); - const auto dp_delta = p * (1.0 - p); // chain rule: d p / d param - const auto dlogl_delta = dlogl_dp * dp_delta; // grad wrt params - - auto matrix_itr = std::cbegin(matrix[i]); - const auto data_end = data + n_factors; - for (auto data_itr = data; data_itr != data_end; ++data_itr) - *data_itr += (*matrix_itr++) * dlogl_delta; - - const auto digamma_delta = c.digamma_a_b - digamma(n + a + b); - const auto dphi_term1 = p * (digamma_y_a - digamma_a + digamma_delta); - const auto dphi_term2 = - (1.0 - p) * (digamma_n_y_b - digamma_b + digamma_delta); - - grad_phi += dphi_term1 + dphi_term2; + double disp_deriv = 0.0; + for (std::size_t g_idx = 0; g_idx < n_groups; ++g_idx) { + const auto p = p_v[g_idx]; + const auto one_minus_p = 1.0 - p; + + double deriv = 0.0; + + const auto denom_term1 = one_minus_phi * p; + const auto &cumul_y = cumul[g_idx].m_counts; + const auto y_lim = std::size(cumul_y); + for (auto k = 0u; k < y_lim; ++k) { + const auto common_factor = cumul_y[k] / (denom_term1 + phi * k); + deriv += common_factor; + disp_deriv += (k - p) * common_factor; + } + + const auto denom_term2 = one_minus_phi * one_minus_p; + const auto &cumul_d = cumul[g_idx].d_counts; + const auto d_lim = std::size(cumul_d); + for (auto k = 0u; k < d_lim; ++k) { + const auto common_factor = cumul_d[k] / (denom_term2 + phi * k); + deriv -= common_factor; + disp_deriv += (k - one_minus_p) * common_factor; + } + + const auto &cumul_n = cumul[g_idx].r_counts; + const auto n_lim = std::size(cumul_n); + for (auto k = 0u; k < n_lim; ++k) + disp_deriv -= cumul_n[k] * dispersion_effect[k]; + + const auto &g = groups[g_idx]; + const auto denom_term1_one_minus_p = denom_term1 * one_minus_p; + for (auto fact_idx = 0u; fact_idx < n_factors; ++fact_idx) { + const auto level = g[fact_idx]; + if (level == 0) + continue; + data[fact_idx] += deriv * (denom_term1_one_minus_p * level); + } } - - const auto grad_theta = -grad_phi * phi; - gsl_vector_set(output, n_factors, grad_theta); + gsl_vector_set(output, n_factors, disp_deriv * (phi * one_minus_phi)); } [[nodiscard]] static double @@ -259,14 +184,63 @@ neg_loglik_and_grad(const gsl_vector *params, void *object, double *loglik_val, neg_gradient(params, object, d_loglik_val); } +static void +get_cumulative(const std::vector &group_id, + const std::uint32_t n_groups, const std::vector &mc, + std::vector &cumul) { + const auto n_cols = std::size(mc); + cumul.clear(); + cumul.resize(n_groups); + + const auto comp_cumul = [&](auto get_value, auto get_vector) { + // phase 1: determine max value for each group + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) { + std::uint32_t max_v{}; + for (auto c_idx = 0u; c_idx < n_cols; ++c_idx) { + if (group_id[c_idx] == g_idx) { + const auto val = get_value(mc[c_idx]); + if (val > max_v) + max_v = val; + } + } + get_vector(cumul[g_idx]).resize(max_v, 0); + } + + // phase 2: fill cumulative counts + for (auto c_idx = 0u; c_idx < n_cols; ++c_idx) { + const auto g_idx = group_id[c_idx]; + const auto val = get_value(mc[c_idx]); + auto &vec = get_vector(cumul[g_idx]); + for (auto i = 0u; i < val; ++i) + ++vec[i]; + } + }; + // call the lambda 3 times for m_counts, r_counts, d_counts + comp_cumul( + [](const mcounts &m) { return m.n_meth; }, + [](cumul_counts &c) -> std::vector & { return c.m_counts; }); + + comp_cumul( + [](const mcounts &m) { return m.n_reads; }, + [](cumul_counts &c) -> std::vector & { return c.r_counts; }); + + comp_cumul( + [](const mcounts &m) { return m.n_reads - m.n_meth; }, + [](cumul_counts &c) -> std::vector & { return c.d_counts; }); +} + void fit_regression_model(Regression &r, std::vector &p_estimates, double &dispersion_estimate) { static constexpr auto init_dispersion_param = -2.5; const auto stepsize = Regression::stepsize; const auto max_iter = Regression::max_iter; + const auto n_groups = r.n_groups(); - r.cache.resize(n_groups); // make sure scratch space is allocated + get_cumulative(r.design.group_id, n_groups, r.mc, r.cumul); + set_max_r_count(r); + + r.p_v.resize(n_groups); const std::size_t n_params = r.n_params(); const auto tol = std::sqrt(n_params) * r.n_samples() * Regression::tolerance; From 639aa0c7592eb5dcdd48f161f396300a8fa5dd0f Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Thu, 16 Oct 2025 11:15:43 -0700 Subject: [PATCH 11/33] src/radmeth/radmeth_nano.cpp: adding this file, forgotten previously, and almost identical to regular radmeth.cpp --- src/radmeth/radmeth_nano.cpp | 375 +++++++++++++++++++++++++++++++++++ 1 file changed, 375 insertions(+) create mode 100644 src/radmeth/radmeth_nano.cpp diff --git a/src/radmeth/radmeth_nano.cpp b/src/radmeth/radmeth_nano.cpp new file mode 100644 index 00000000..bca070b9 --- /dev/null +++ b/src/radmeth/radmeth_nano.cpp @@ -0,0 +1,375 @@ +/* Copyright (C) 2013-2025 Andrew D Smith + * + * Author: Andrew D. Smith + * Contributors: Egor Dolzhenko and Guilherme Sena + * + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. + * + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + */ + +#include "radmeth_design.hpp" +#include "radmeth_model_nano.hpp" +#include "radmeth_optimize_nano.hpp" +#include "radmeth_utils.hpp" + +// smithlab headers +#include "GenomicRegion.hpp" +#include "OptionParser.hpp" +#include "smithlab_os.hpp" +#include "smithlab_utils.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +template +[[nodiscard]] static bool +has_low_coverage(const RegressionType ®, const std::size_t test_factor) { + bool cvrd_in_test_fact_smpls = false; + const auto &tcol = reg.design.tmatrix[test_factor]; + for (std::size_t i = 0; i < reg.n_samples() && !cvrd_in_test_fact_smpls; ++i) + cvrd_in_test_fact_smpls = (tcol[i] == 1 && reg.mc[i].n_reads != 0); + + bool cvrd_in_other_smpls = false; + for (std::size_t i = 0; i < reg.n_samples() && !cvrd_in_other_smpls; ++i) + cvrd_in_other_smpls = (tcol[i] != 1 && reg.mc[i].n_reads != 0); + + return !cvrd_in_test_fact_smpls || !cvrd_in_other_smpls; +} + +template +[[nodiscard]] static bool +has_extreme_counts(const RegressionType ®) { + const auto &mc = reg.mc; + + bool full_meth = true; + for (auto i = std::cbegin(mc); i != std::cend(mc) && full_meth; ++i) + full_meth = (i->n_reads == i->n_meth); + + bool no_meth = true; + for (auto i = std::cbegin(mc); i != std::cend(mc) && no_meth; ++i) + no_meth = (i->n_meth == 0.0); + + return full_meth || no_meth; +} + +enum class row_status : std::uint8_t { + ok, + na, + na_low_cov, + na_extreme_cnt, +}; + +template +static void +radmeth(const bool show_progress, const bool more_na_info, + const std::uint32_t n_threads, const std::string &table_filename, + const std::string &outfile, const RegressionType &alt_model, + const RegressionType &null_model, const std::uint32_t test_factor_idx) { + static constexpr auto buf_size = 1024; + static constexpr auto max_lines = 16384; + + // ADS: open the data table file + std::ifstream table_file(table_filename); + if (!table_file) + throw std::runtime_error("could not open file: " + table_filename); + + // Make sure that the first line of the proportion table file contains + // names of the samples. Throw an exception if the names or their order + // in the proportion table does not match those in the full design matrix. + 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.design, sample_names_header)) { + // clang-format off + const auto message = + R"( +does not match factor names or their order in the design matrix. Check +that the design matrix and the proportion table are correctly formatted. +)"; + // clang-format on + throw std::runtime_error("header:\n" + sample_names_header + message); + } + + const std::size_t n_samples = alt_model.n_samples(); + + std::ofstream out(outfile); + if (!out) + throw std::runtime_error("failed to open output file: " + outfile); + + file_progress progress{table_filename}; + + std::vector> bufs(max_lines, + std::vector(buf_size, 0)); + std::vector n_bytes(max_lines, 0); + std::vector lines(max_lines); + + std::vector alt_models(n_threads, alt_model); + std::vector null_models(n_threads, null_model); + + const auto n_groups = alt_model.n_groups(); + + // Iterate over rows in the file and do llr test on proportions from + // each. Do this in sets of rows to avoid having to spawn too many threads. + while (true) { + std::uint32_t n_lines = 0; + while (n_lines < max_lines && std::getline(table_file, lines[n_lines])) + ++n_lines; + if (n_lines == 0) + break; + + std::vector threads; + for (auto thread_id = 0u; thread_id < n_threads; ++thread_id) { + threads.emplace_back([&, thread_id] { + std::vector p_estim_alt; + std::vector p_estim_null; + double phi_estim_alt{}; + double phi_estim_null{}; + + auto &t_alt_model = alt_models[thread_id]; + auto &t_null_model = null_models[thread_id]; + for (auto b = 0u; b < n_lines; ++b) { + // ADS: rows done by different threads are interleaved because the + // difficult (e.g., high-coverage) rows can be consecutive and this + // balances work better. + if (b % n_threads != thread_id) + continue; + t_alt_model.parse(lines[b]); + if (t_alt_model.props_size() != n_samples) + throw std::runtime_error("found row with wrong number of columns"); + + const auto p_val_status = [&]() -> std::tuple { + // Skip the test if (1) no coverage in all cases or in all + // controls, or (2) the site is completely methylated or + // completely unmethylated across all samples. + if (has_low_coverage(t_alt_model, test_factor_idx)) + return std::tuple{1.0, row_status::na_low_cov}; + + if (has_extreme_counts(t_alt_model)) + return std::tuple{1.0, row_status::na_extreme_cnt}; + + fit_regression_model_nano(t_alt_model, p_estim_alt, phi_estim_alt); + + t_null_model.mc = t_alt_model.mc; + t_null_model.rowname = t_alt_model.rowname; + + fit_regression_model_nano(t_null_model, p_estim_null, + phi_estim_null); + + const double p_value = + llr_test(t_null_model.max_loglik, t_alt_model.max_loglik); + + return (p_value != p_value) ? std::tuple{1.0, row_status::na} + : std::tuple{p_value, row_status::ok}; + }(); + // ADS: avoid capture structured binding in C++17 + const auto p_val = std::get<0>(p_val_status); + const auto status = std::get<1>(p_val_status); + + n_bytes[b] = [&] { + // clang-format off + const int n_prefix_bytes = + std::sprintf(bufs[b].data(), "%s\t", t_alt_model.rowname.data()); + // clang-format on + if (n_prefix_bytes < 0) + return n_prefix_bytes; + + auto cursor = bufs[b].data() + n_prefix_bytes; + + const int n_pval_bytes = [&] { + if (status == row_status::ok) + return std::sprintf(cursor, "%.6g", p_val); + if (!more_na_info || status == row_status::na) + return std::sprintf(cursor, "NA"); + if (status == row_status::na_extreme_cnt) + return std::sprintf(cursor, "NA_EXTREME_CNT"); + // if (status == row_status::na_low_cov) + return std::sprintf(cursor, "NA_LOW_COV"); + }(); + + if (n_pval_bytes < 0) + return n_pval_bytes; + + cursor += n_pval_bytes; + + const int n_param_bytes = [&] { + std::int32_t n_bytes = 0; + if (p_estim_alt.empty()) + p_estim_alt.resize(n_groups, 0.0); + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) { + const int n = std::sprintf(cursor, "\t%f", p_estim_alt[g_idx]); + cursor += n; + if (n < 0) + return n; + n_bytes += n; + } + const auto od = overdispersion_factor(n_samples, phi_estim_alt); + const int n = std::sprintf(cursor, "\t%f\n", od); + if (n < 0) + return n; + n_bytes += n; + return n_bytes; + }(); + + if (n_param_bytes < 0) + return n_param_bytes; + + return n_prefix_bytes + n_pval_bytes + n_param_bytes; + }(); + } + }); + } + + for (auto &thread : threads) + thread.join(); + + for (auto i = 0u; i < n_lines; ++i) { + if (n_bytes[i] < 0) + throw std::runtime_error("failed to write to output buffer"); + out.write(bufs[i].data(), n_bytes[i]); + } + + if (show_progress) + progress(table_file); + + if (n_lines < n_threads) + break; + } + + if (show_progress) + progress(table_file); +} + +int +main_radmeth_nano(int argc, char *argv[]) { + try { + static const std::string description = + "calculate differential methylation scores (optimized for nanopore)"; + + std::string outfile; + std::string test_factor; + bool verbose{false}; + bool show_progress{false}; + bool more_na_info{false}; + std::uint32_t n_threads{1}; + + /****************** COMMAND LINE OPTIONS ********************/ + OptionParser opt_parse(strip_path(argv[0]), description, + " "); + opt_parse.set_show_defaults(); + opt_parse.add_opt("out", 'o', "output file", true, outfile); + opt_parse.add_opt("threads", 't', "number of threads to use", false, + n_threads); + opt_parse.add_opt("verbose", 'v', "print more run info", false, verbose); + opt_parse.add_opt("progress", '\0', "show progress", false, show_progress); + opt_parse.add_opt( + "na-info", 'n', + "if a p-value is not calculated, print NAs in more detail: " + "low count (NA_LOW_COV) extreme values (NA_EXTREME_CNT) or " + "numerical errors in likelihood ratios (NA)", + false, more_na_info); + opt_parse.add_opt("factor", 'f', "a factor to test", true, test_factor); + opt_parse.add_opt("tolerance", '\0', + "numerical tolerance to test convergence", false, + RegressionNano::tolerance); + opt_parse.add_opt("max-iter", '\0', + "max iterations when estimating parameters", false, + RegressionNano::max_iter); + std::vector leftover_args; + opt_parse.parse(argc, argv, leftover_args); + if (argc == 1 || opt_parse.help_requested()) { + std::cerr << opt_parse.help_message() << std::endl + << opt_parse.about_message() << std::endl; + return EXIT_SUCCESS; + } + if (opt_parse.about_requested()) { + std::cerr << opt_parse.about_message() << std::endl; + return EXIT_SUCCESS; + } + if (opt_parse.option_missing()) { + std::cerr << opt_parse.option_missing_message() << std::endl; + return EXIT_SUCCESS; + } + if (leftover_args.size() != 2) { + std::cerr << opt_parse.help_message() << std::endl; + return EXIT_SUCCESS; + } + const std::string design_filename(leftover_args.front()); + const std::string table_filename(leftover_args.back()); + /****************** END COMMAND LINE OPTIONS *****************/ + + if (verbose) + std::cerr << "design table filename: " << design_filename << "\n\n"; + + Design design = Design::read_design(design_filename); + if (verbose) + std::cerr << "Alternate model:\n" << design << '\n'; + + // Check that provided test factor name exists and find its index. + const auto test_factor_idx = design.get_test_factor_idx(test_factor); + + ensure_sample_order(table_filename, design); + + // verify that the design includes more than one level for the + // test factor + if (!design.has_two_values(test_factor_idx)) { + const auto first_level = design.matrix[0][test_factor_idx]; + throw std::runtime_error("test factor only one level: " + test_factor + + ", " + std::to_string(first_level)); + } + + const Design null_design = design.drop_factor(test_factor_idx); + if (verbose) + std::cerr << "Null model:\n" << null_design << '\n'; + + // clang-format off + if (verbose) + std::cerr << "Output columns:\n" + << "(1) chrom\n" + << "(2) position\n" + << "(3) strand\n" + << "(4) cytosine context\n" + << "(5) p-value\n" + << "(6) estimated methylation 0\n" + << "(7) estimated methylation 1\n" + << "(8) overdispersion factor (variance above binomial)\n" + << "estimated methylation is for test factor value (0 or 1)\n" + << '\n'; + // clang-format on + + RegressionNano alt_model; + alt_model.design = design; + RegressionNano null_model; + null_model.design = null_design; + + const auto start_time = std::chrono::steady_clock::now(); + radmeth(show_progress, more_na_info, n_threads, table_filename, outfile, + alt_model, null_model, test_factor_idx); + const auto stop_time = std::chrono::steady_clock::now(); + + if (verbose) + std::cerr << "[total time: " << format_duration(stop_time - start_time) + << "]\n"; + } + catch (const std::exception &e) { + std::cerr << e.what() << '\n'; + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +} From be62144238a1232670710481c019b94d25d0aa25 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Thu, 16 Oct 2025 17:17:23 -0700 Subject: [PATCH 12/33] src/radmeth/radmeth_optimize.cpp and _nano.cpp: adding a const --- src/radmeth/radmeth_optimize.cpp | 2 +- src/radmeth/radmeth_optimize_nano.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/radmeth/radmeth_optimize.cpp b/src/radmeth/radmeth_optimize.cpp index 77f6aaba..9612f902 100644 --- a/src/radmeth/radmeth_optimize.cpp +++ b/src/radmeth/radmeth_optimize.cpp @@ -120,7 +120,7 @@ gradient(const gsl_vector *params, Regression ®, gsl_vector *output) { // init output to zero for all factors gsl_vector_set_all(output, 0.0); - auto &data = output->data; + const auto &data = output->data; double disp_deriv = 0.0; for (std::size_t g_idx = 0; g_idx < n_groups; ++g_idx) { diff --git a/src/radmeth/radmeth_optimize_nano.cpp b/src/radmeth/radmeth_optimize_nano.cpp index 25024337..c3a2f4ac 100644 --- a/src/radmeth/radmeth_optimize_nano.cpp +++ b/src/radmeth/radmeth_optimize_nano.cpp @@ -198,7 +198,7 @@ gradient_nano(const gsl_vector *params, RegressionNano ®, get_cache_digamma(groups[g_idx], params, phi, cache[g_idx]); gsl_vector_set_all(output, 0.0); // init output to zero for all factors - auto &data = output->data; + const auto &data = output->data; const auto &group_id = reg.design.group_id; double grad_phi = 0.0; From 5830de33f4b50e7e680af89d62e65eab88a5aa40 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Thu, 16 Oct 2025 17:19:46 -0700 Subject: [PATCH 13/33] src/radmeth/radmeth_nano.cpp: minor cleanup --- src/radmeth/radmeth_nano.cpp | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/src/radmeth/radmeth_nano.cpp b/src/radmeth/radmeth_nano.cpp index bca070b9..87b1fd72 100644 --- a/src/radmeth/radmeth_nano.cpp +++ b/src/radmeth/radmeth_nano.cpp @@ -1,7 +1,6 @@ -/* Copyright (C) 2013-2025 Andrew D Smith +/* Copyright (C) 2025 Andrew D Smith * * Author: Andrew D. Smith - * Contributors: Egor Dolzhenko and Guilherme Sena * * This program is free software: you can redistribute it and/or modify it * under the terms of the GNU General Public License as published by the Free @@ -77,10 +76,11 @@ enum class row_status : std::uint8_t { template static void -radmeth(const bool show_progress, const bool more_na_info, - const std::uint32_t n_threads, const std::string &table_filename, - const std::string &outfile, const RegressionType &alt_model, - const RegressionType &null_model, const std::uint32_t test_factor_idx) { +radmeth_nano(const bool show_progress, const bool more_na_info, + const std::uint32_t n_threads, const std::string &table_filename, + const std::string &outfile, const RegressionType &alt_model, + const RegressionType &null_model, + const std::uint32_t test_factor_idx) { static constexpr auto buf_size = 1024; static constexpr auto max_lines = 16384; @@ -260,7 +260,7 @@ int main_radmeth_nano(int argc, char *argv[]) { try { static const std::string description = - "calculate differential methylation scores (optimized for nanopore)"; + "calculate differential methylation scores for nanopore data"; std::string outfile; std::string test_factor; @@ -335,12 +335,11 @@ main_radmeth_nano(int argc, char *argv[]) { } const Design null_design = design.drop_factor(test_factor_idx); - if (verbose) - std::cerr << "Null model:\n" << null_design << '\n'; // clang-format off if (verbose) - std::cerr << "Output columns:\n" + std::cerr << "Null model:\n" << null_design << '\n' + << "Output columns:\n" << "(1) chrom\n" << "(2) position\n" << "(3) strand\n" @@ -359,8 +358,8 @@ main_radmeth_nano(int argc, char *argv[]) { null_model.design = null_design; const auto start_time = std::chrono::steady_clock::now(); - radmeth(show_progress, more_na_info, n_threads, table_filename, outfile, - alt_model, null_model, test_factor_idx); + radmeth_nano(show_progress, more_na_info, n_threads, table_filename, + outfile, alt_model, null_model, test_factor_idx); const auto stop_time = std::chrono::steady_clock::now(); if (verbose) From d9ab6699accd8d2f52d793fcffb40c676d782a1c Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Thu, 16 Oct 2025 19:10:21 -0700 Subject: [PATCH 14/33] src/radmeth/radmeth.cpp: ensuring the reported model design for the alt model has absent samples removed --- src/radmeth/radmeth.cpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/radmeth/radmeth.cpp b/src/radmeth/radmeth.cpp index 2a8f215d..c5642742 100644 --- a/src/radmeth/radmeth.cpp +++ b/src/radmeth/radmeth.cpp @@ -293,20 +293,20 @@ main_radmeth(int argc, char *argv[]) { std::vector leftover_args; opt_parse.parse(argc, argv, leftover_args); if (argc == 1 || opt_parse.help_requested()) { - std::cerr << opt_parse.help_message() << std::endl - << opt_parse.about_message() << std::endl; + std::cerr << opt_parse.help_message() << '\n' + << opt_parse.about_message() << '\n'; return EXIT_SUCCESS; } if (opt_parse.about_requested()) { - std::cerr << opt_parse.about_message() << std::endl; + std::cerr << opt_parse.about_message() << '\n'; return EXIT_SUCCESS; } if (opt_parse.option_missing()) { - std::cerr << opt_parse.option_missing_message() << std::endl; + std::cerr << opt_parse.option_missing_message() << '\n'; return EXIT_SUCCESS; } if (leftover_args.size() != 2) { - std::cerr << opt_parse.help_message() << std::endl; + std::cerr << opt_parse.help_message() << '\n'; return EXIT_SUCCESS; } const std::string design_filename(leftover_args.front()); @@ -317,14 +317,15 @@ main_radmeth(int argc, char *argv[]) { std::cerr << "design table filename: " << design_filename << "\n\n"; Design design = Design::read_design(design_filename); - if (verbose) - std::cerr << "Alternate model:\n" << design << '\n'; // Check that provided test factor name exists and find its index. const auto test_factor_idx = design.get_test_factor_idx(test_factor); ensure_sample_order(table_filename, design); + if (verbose) + std::cerr << "Alternate model:\n" << design << '\n'; + // verify that the design includes more than one level for the // test factor if (!design.has_two_values(test_factor_idx)) { From 07dc1fdb58b36f244a3c4184dca9fa10331bccf1 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Thu, 16 Oct 2025 19:10:54 -0700 Subject: [PATCH 15/33] src/radmeth/radmeth_nano.cpp: ensuring the reported model design for the alt model has absent samples removed --- src/radmeth/radmeth_nano.cpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/radmeth/radmeth_nano.cpp b/src/radmeth/radmeth_nano.cpp index 87b1fd72..7451e2dd 100644 --- a/src/radmeth/radmeth_nano.cpp +++ b/src/radmeth/radmeth_nano.cpp @@ -294,20 +294,20 @@ main_radmeth_nano(int argc, char *argv[]) { std::vector leftover_args; opt_parse.parse(argc, argv, leftover_args); if (argc == 1 || opt_parse.help_requested()) { - std::cerr << opt_parse.help_message() << std::endl - << opt_parse.about_message() << std::endl; + std::cerr << opt_parse.help_message() << '\n' + << opt_parse.about_message() << '\n'; return EXIT_SUCCESS; } if (opt_parse.about_requested()) { - std::cerr << opt_parse.about_message() << std::endl; + std::cerr << opt_parse.about_message() << '\n'; return EXIT_SUCCESS; } if (opt_parse.option_missing()) { - std::cerr << opt_parse.option_missing_message() << std::endl; + std::cerr << opt_parse.option_missing_message() << '\n'; return EXIT_SUCCESS; } if (leftover_args.size() != 2) { - std::cerr << opt_parse.help_message() << std::endl; + std::cerr << opt_parse.help_message() << '\n'; return EXIT_SUCCESS; } const std::string design_filename(leftover_args.front()); @@ -318,14 +318,15 @@ main_radmeth_nano(int argc, char *argv[]) { std::cerr << "design table filename: " << design_filename << "\n\n"; Design design = Design::read_design(design_filename); - if (verbose) - std::cerr << "Alternate model:\n" << design << '\n'; // Check that provided test factor name exists and find its index. const auto test_factor_idx = design.get_test_factor_idx(test_factor); ensure_sample_order(table_filename, design); + if (verbose) + std::cerr << "Alternate model:\n" << design << '\n'; + // verify that the design includes more than one level for the // test factor if (!design.has_two_values(test_factor_idx)) { From b53c7ba1828c70a38e490916bd4b5a980c639c6a Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 17 Oct 2025 11:56:09 -0700 Subject: [PATCH 16/33] src/radmeth/radmeth_optimize_nano.cpp: making the chebyschev approx full order for now --- src/radmeth/radmeth_optimize_nano.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/radmeth/radmeth_optimize_nano.cpp b/src/radmeth/radmeth_optimize_nano.cpp index c3a2f4ac..e9f5a48c 100644 --- a/src/radmeth/radmeth_optimize_nano.cpp +++ b/src/radmeth/radmeth_optimize_nano.cpp @@ -75,7 +75,7 @@ static constexpr std::array psi2_cs = { 0.0000000000000045, -0.0000000000000009, 0.0000000000000002, - -0.0000000000000000 , + -0.0000000000000000 }; // clang-format on @@ -156,15 +156,15 @@ log_likelihood_nano(const gsl_vector *params, RegressionNano ®) { // digamma for x non-negative static double digamma(const double y) { - static constexpr auto psi_order = 7; // max=22; - static constexpr auto apsi_order = 7; // max=15; + static constexpr auto psi1_order = 22; // max=22; + static constexpr auto psi2_order = 15; // max=15; if (y >= 2.0) { const auto t = 8.0 / (y * y) - 1.0; - return std::log(y) - 0.5 / y + chebyschev(psi2_cs, apsi_order, t); + return std::log(y) - 0.5 / y + chebyschev(psi2_cs, psi2_order, t); } if (y < 1.0) - return -1.0 / y + chebyschev(psi1_cs, psi_order, 2.0 * y - 1.0); - return chebyschev(psi1_cs, psi_order, 2.0 * (y - 1.0) - 1.0); + return -1.0 / y + chebyschev(psi1_cs, psi1_order, 2.0 * y - 1.0); + return chebyschev(psi1_cs, psi1_order, 2.0 * (y - 1.0) - 1.0); } static void From 780ebc2c9bf8a3488875e529b00ff1c50e31cabf Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Fri, 17 Oct 2025 16:31:25 -0700 Subject: [PATCH 17/33] src/radmeth/radmeth_optimize.cpp: minor tweaks for speed and clarity --- src/radmeth/radmeth_optimize.cpp | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/radmeth/radmeth_optimize.cpp b/src/radmeth/radmeth_optimize.cpp index 9612f902..0863c34c 100644 --- a/src/radmeth/radmeth_optimize.cpp +++ b/src/radmeth/radmeth_optimize.cpp @@ -61,8 +61,10 @@ static inline auto cache_dispersion_effect(Regression ®, const double phi) { const std::size_t max_k = reg.max_r_count; auto &cache = reg.cache; - for (std::size_t k = 0; k < max_k; ++k) - cache[k] = (k - 1.0) / (1.0 + phi * (k - 1.0)); + double j = -1.0; + const auto lim = std::cbegin(cache) + max_k; + for (auto it = std::begin(cache); it != lim; ++it, ++j) + *it = j / (1.0 + phi * j); } [[nodiscard]] static double @@ -86,16 +88,19 @@ log_likelihood(const gsl_vector *params, Regression ®) { const auto term1 = one_minus_phi * p; const auto &cumul_y = cumul[g_idx].m_counts; - for (std::size_t k = 0; k < std::size(cumul_y); ++k) + std::size_t sz = std::size(cumul_y); + for (std::size_t k = 0; k < sz; ++k) log_lik += cumul_y[k] * std::log(term1 + phi * k); const auto term2 = one_minus_phi * one_minus_p; const auto &cumul_d = cumul[g_idx].d_counts; - for (std::size_t k = 0; k < std::size(cumul_d); ++k) + sz = std::size(cumul_d); + for (std::size_t k = 0; k < sz; ++k) log_lik += cumul_d[k] * std::log(term2 + phi * k); const auto &cumul_n = cumul[g_idx].r_counts; - for (std::size_t k = 0; k < std::size(cumul_n); ++k) + sz = std::size(cumul_n); + for (std::size_t k = 0; k < sz; ++k) log_lik -= cumul_n[k] * log1p_fact_v[k]; } return log_lik; From bc71c535016e3104eb685a1bd2ed5a211ed4116a Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 19 Oct 2025 17:49:17 -0700 Subject: [PATCH 18/33] src/radmeth/radmeth_model.cpp: removed during renaming --- src/radmeth/radmeth_model.cpp | 72 ----------------------------------- 1 file changed, 72 deletions(-) delete mode 100644 src/radmeth/radmeth_model.cpp diff --git a/src/radmeth/radmeth_model.cpp b/src/radmeth/radmeth_model.cpp deleted file mode 100644 index db39b8c7..00000000 --- a/src/radmeth/radmeth_model.cpp +++ /dev/null @@ -1,72 +0,0 @@ -/* Copyright (C) 2025 Andrew D Smith - * - * Authors: Andrew D. Smith - * - * This program is free software: you can redistribute it and/or - * modify it under the terms of the GNU General Public License as - * published by the Free Software Foundation, either version 3 of the - * License, or (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, but - * WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * General Public License for more details. - */ - -#include "radmeth_model.hpp" - -#include -#include -#include -#include -#include -#include -#include - -double Regression::tolerance = 1e-4; -double Regression::stepsize = 0.01; -std::uint32_t Regression::max_iter = 250; - -void -Regression::parse(const std::string &line) { - const auto first_ws = line.find_first_of(" \t"); - - // Parse the row name (must be like: "chr:position:strand:context") - bool failed = (first_ws == std::string::npos); - - auto field_s = line.data(); - auto field_e = line.data() + first_ws; - rowname = std::string{field_s, field_e}; - std::replace(std::begin(rowname), std::end(rowname), ':', '\t'); - if (failed) - throw std::runtime_error("failed to parse label from:\n" + line); - - // Parse the counts of total reads and methylated reads - const auto is_sep = [](const auto x) { return std::isspace(x); }; - const auto not_sep = [](const auto x) { return std::isdigit(x); }; - - const auto line_end = line.data() + std::size(line); - mcounts mc1{}; - mc.clear(); - while (field_e != line_end) { - // get the total count - field_s = std::find_if(field_e + 1, line_end, not_sep); - field_e = std::find_if(field_s + 1, line_end, is_sep); - { - const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_reads); - failed = failed || (ptr != field_e); - } - - field_s = std::find_if(field_e + 1, line_end, not_sep); - field_e = std::find_if(field_s + 1, line_end, is_sep); - { - const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_meth); - failed = failed || (ptr != field_e); - } - - mc.push_back(mc1); - } - - if (failed) - throw std::runtime_error("failed to parse counts from:\n" + line); -} From 001adb22111b9832faf1d3211d90557d54867dac Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 19 Oct 2025 17:51:39 -0700 Subject: [PATCH 19/33] src/radmeth/radmeth_optimize.hcpp: renamed --- src/radmeth/radmeth_optimize.cpp | 307 ------------------------------- src/radmeth/radmeth_optimize.hpp | 30 --- 2 files changed, 337 deletions(-) delete mode 100644 src/radmeth/radmeth_optimize.cpp delete mode 100644 src/radmeth/radmeth_optimize.hpp diff --git a/src/radmeth/radmeth_optimize.cpp b/src/radmeth/radmeth_optimize.cpp deleted file mode 100644 index 0863c34c..00000000 --- a/src/radmeth/radmeth_optimize.cpp +++ /dev/null @@ -1,307 +0,0 @@ -/* Copyright (C) 2025 Andrew D. Smith - * - * Author: Andrew D. Smith - * - * This program is free software: you can redistribute it and/or modify it - * under the terms the GNU General Public License as published by the Free - * Software Foundation, either version 3 of the License, or (at your option) - * any later version. - * - * This program is distributed in the hope that it will be useful, but WITHOUT - * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for - * more details. - */ - -#include "radmeth_optimize.hpp" -#include "radmeth_model.hpp" - -#include -#include - -#include -#include -#include -#include - -[[nodiscard]] static inline double -logistic(const double x) { - return 1.0 / (1.0 / std::exp(x) + 1.0); -} - -template -[[nodiscard]] static double -get_p(const std::vector &v, const gsl_vector *params) { - const auto a = v.data(); - return logistic(std::inner_product(a, a + std::size(v), params->data, 0.0)); -} - -static inline auto -set_max_r_count(Regression ®) { - const auto &cumul = reg.cumul; - const auto max_itr = std::max_element( - std::cbegin(cumul), std::cend(cumul), [](const auto &a, const auto &b) { - return std::size(a.r_counts) < std::size(b.r_counts); - }); - reg.max_r_count = std::size(max_itr->r_counts); - // ADS: avoid the realloc that can happen even for resize(smaller_size) - if (reg.max_r_count > std::size(reg.cache)) - reg.cache.resize(reg.max_r_count); -} - -static inline auto -cache_log1p_factors(Regression ®, const double phi) { - const std::size_t max_k = reg.max_r_count; - auto &cache = reg.cache; - for (std::size_t k = 0; k < max_k; ++k) - cache[k] = std::log1p(phi * (k - 1.0)); -} - -static inline auto -cache_dispersion_effect(Regression ®, const double phi) { - const std::size_t max_k = reg.max_r_count; - auto &cache = reg.cache; - double j = -1.0; - const auto lim = std::cbegin(cache) + max_k; - for (auto it = std::begin(cache); it != lim; ++it, ++j) - *it = j / (1.0 + phi * j); -} - -[[nodiscard]] static double -log_likelihood(const gsl_vector *params, Regression ®) { - const auto phi = logistic(gsl_vector_get(params, reg.design.n_factors())); - const auto one_minus_phi = 1.0 - phi; - - const auto n_groups = reg.n_groups(); - const auto &groups = reg.design.groups; - const auto &cumul = reg.cumul; - - // ADS: precompute the log1p(phi * (k - 1.0)) values, which are reused for - // each group. - cache_log1p_factors(reg, phi); - const auto &log1p_fact_v = reg.cache; - - double log_lik = 0.0; - for (std::size_t g_idx = 0; g_idx < n_groups; ++g_idx) { - const auto p = get_p(groups[g_idx], params); - const auto one_minus_p = 1.0 - p; - - const auto term1 = one_minus_phi * p; - const auto &cumul_y = cumul[g_idx].m_counts; - std::size_t sz = std::size(cumul_y); - for (std::size_t k = 0; k < sz; ++k) - log_lik += cumul_y[k] * std::log(term1 + phi * k); - - const auto term2 = one_minus_phi * one_minus_p; - const auto &cumul_d = cumul[g_idx].d_counts; - sz = std::size(cumul_d); - for (std::size_t k = 0; k < sz; ++k) - log_lik += cumul_d[k] * std::log(term2 + phi * k); - - const auto &cumul_n = cumul[g_idx].r_counts; - sz = std::size(cumul_n); - for (std::size_t k = 0; k < sz; ++k) - log_lik -= cumul_n[k] * log1p_fact_v[k]; - } - return log_lik; -} - -static void -gradient(const gsl_vector *params, Regression ®, gsl_vector *output) { - const auto n_factors = reg.design.n_factors(); - const auto phi = logistic(gsl_vector_get(params, n_factors)); - const auto one_minus_phi = 1.0 - phi; - - const auto n_groups = reg.n_groups(); - const auto &groups = reg.design.groups; - const auto &cumul = reg.cumul; - - auto &p_v = reg.p_v; // ADS: reusing scratch space - for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) - p_v[g_idx] = get_p(groups[g_idx], params); - - cache_dispersion_effect(reg, phi); - const auto &dispersion_effect = reg.cache; // (k-1)/(1 + phi(k-1)) - - // init output to zero for all factors - gsl_vector_set_all(output, 0.0); - const auto &data = output->data; - - double disp_deriv = 0.0; - for (std::size_t g_idx = 0; g_idx < n_groups; ++g_idx) { - const auto p = p_v[g_idx]; - const auto one_minus_p = 1.0 - p; - - double deriv = 0.0; - - const auto denom_term1 = one_minus_phi * p; - const auto &cumul_y = cumul[g_idx].m_counts; - const auto y_lim = std::size(cumul_y); - for (auto k = 0u; k < y_lim; ++k) { - const auto common_factor = cumul_y[k] / (denom_term1 + phi * k); - deriv += common_factor; - disp_deriv += (k - p) * common_factor; - } - - const auto denom_term2 = one_minus_phi * one_minus_p; - const auto &cumul_d = cumul[g_idx].d_counts; - const auto d_lim = std::size(cumul_d); - for (auto k = 0u; k < d_lim; ++k) { - const auto common_factor = cumul_d[k] / (denom_term2 + phi * k); - deriv -= common_factor; - disp_deriv += (k - one_minus_p) * common_factor; - } - - const auto &cumul_n = cumul[g_idx].r_counts; - const auto n_lim = std::size(cumul_n); - for (auto k = 0u; k < n_lim; ++k) - disp_deriv -= cumul_n[k] * dispersion_effect[k]; - - const auto &g = groups[g_idx]; - const auto denom_term1_one_minus_p = denom_term1 * one_minus_p; - for (auto fact_idx = 0u; fact_idx < n_factors; ++fact_idx) { - const auto level = g[fact_idx]; - if (level == 0) - continue; - data[fact_idx] += deriv * (denom_term1_one_minus_p * level); - } - } - gsl_vector_set(output, n_factors, disp_deriv * (phi * one_minus_phi)); -} - -[[nodiscard]] static double -neg_loglik(const gsl_vector *params, void *object) { - auto reg = static_cast(object); - return -log_likelihood(params, *reg); -} - -static void -neg_gradient(const gsl_vector *params, void *object, gsl_vector *output) { - auto reg = static_cast(object); - gradient(params, *reg, output); - gsl_vector_scale(output, -1.0); -} - -static void -neg_loglik_and_grad(const gsl_vector *params, void *object, double *loglik_val, - gsl_vector *d_loglik_val) { - *loglik_val = neg_loglik(params, object); - neg_gradient(params, object, d_loglik_val); -} - -static void -get_cumulative(const std::vector &group_id, - const std::uint32_t n_groups, const std::vector &mc, - std::vector &cumul) { - const auto n_cols = std::size(mc); - cumul.clear(); - cumul.resize(n_groups); - - const auto comp_cumul = [&](auto get_value, auto get_vector) { - // phase 1: determine max value for each group - for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) { - std::uint32_t max_v{}; - for (auto c_idx = 0u; c_idx < n_cols; ++c_idx) { - if (group_id[c_idx] == g_idx) { - const auto val = get_value(mc[c_idx]); - if (val > max_v) - max_v = val; - } - } - get_vector(cumul[g_idx]).resize(max_v, 0); - } - - // phase 2: fill cumulative counts - for (auto c_idx = 0u; c_idx < n_cols; ++c_idx) { - const auto g_idx = group_id[c_idx]; - const auto val = get_value(mc[c_idx]); - auto &vec = get_vector(cumul[g_idx]); - for (auto i = 0u; i < val; ++i) - ++vec[i]; - } - }; - // call the lambda 3 times for m_counts, r_counts, d_counts - comp_cumul( - [](const mcounts &m) { return m.n_meth; }, - [](cumul_counts &c) -> std::vector & { return c.m_counts; }); - - comp_cumul( - [](const mcounts &m) { return m.n_reads; }, - [](cumul_counts &c) -> std::vector & { return c.r_counts; }); - - comp_cumul( - [](const mcounts &m) { return m.n_reads - m.n_meth; }, - [](cumul_counts &c) -> std::vector & { return c.d_counts; }); -} - -void -fit_regression_model(Regression &r, std::vector &p_estimates, - double &dispersion_estimate) { - static constexpr auto init_dispersion_param = -2.5; - const auto stepsize = Regression::stepsize; - const auto max_iter = Regression::max_iter; - - const auto n_groups = r.n_groups(); - get_cumulative(r.design.group_id, n_groups, r.mc, r.cumul); - set_max_r_count(r); - - r.p_v.resize(n_groups); - - const std::size_t n_params = r.n_params(); - const auto tol = std::sqrt(n_params) * r.n_samples() * Regression::tolerance; - // clang-format off - auto loglik_bundle = gsl_multimin_function_fdf{ - &neg_loglik, // objective function - &neg_gradient, // gradient - &neg_loglik_and_grad, // combined obj and grad - n_params, // number of model params - static_cast(&r) // parameters for objective and gradient functions - }; - // clang-format on - - // set the parameters: zero for "p" parameters and the final one for - // dispersion using the constant - auto params = gsl_vector_alloc(n_params); - gsl_vector_set_all(params, 0.0); - gsl_vector_set(params, n_params - 1, init_dispersion_param); - - // Alternatives: - // - gsl_multimin_fdfminimizer_conjugate_pr - // - gsl_multimin_fdfminimizer_conjugate_fr - // - gsl_multimin_fdfminimizer_vector_bfgs2 - // - gsl_multimin_fdfminimizer_steepest_descent - const auto minimizer = gsl_multimin_fdfminimizer_conjugate_pr; - auto s = gsl_multimin_fdfminimizer_alloc(minimizer, n_params); - - gsl_multimin_fdfminimizer_set(s, &loglik_bundle, params, stepsize, tol); - - int status = 0; - std::size_t iter = 0; - do { - status = gsl_multimin_fdfminimizer_iterate(s); // one iter and get status - if (status) - break; - // check status from gradient - status = gsl_multimin_test_gradient(s->gradient, tol); - } while (status == GSL_CONTINUE && ++iter < max_iter); - - /// ADS: the condition below might not always pass even if we are doing - /// ok. It's not clear how to check for failure. - - // if (status != GSL_SUCCESS) - // throw std::runtime_error("failed to fit model parameters"); - - const auto param_estimates = gsl_multimin_fdfminimizer_x(s); - - const auto &groups = r.design.groups; - p_estimates.clear(); - for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) - p_estimates.push_back(get_p(groups[g_idx], param_estimates)); - const auto disp_param = gsl_vector_get(param_estimates, n_params - 1); - dispersion_estimate = 1.0 / std::exp(disp_param); - - r.max_loglik = log_likelihood(param_estimates, r); - - gsl_multimin_fdfminimizer_free(s); - gsl_vector_free(params); -} diff --git a/src/radmeth/radmeth_optimize.hpp b/src/radmeth/radmeth_optimize.hpp deleted file mode 100644 index d0209fb2..00000000 --- a/src/radmeth/radmeth_optimize.hpp +++ /dev/null @@ -1,30 +0,0 @@ -/* Copyright (C) 2013-2023 University of Southern California and - * Egor Dolzhenko - * Andrew D Smith - * Guilherme Sena - * - * Authors: Andrew D. Smith and Egor Dolzhenko and Guilherme Sena - * - * This program is free software: you can redistribute it and/or - * modify it under the terms of the GNU General Public License as - * published by the Free Software Foundation, either version 3 of the - * License, or (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, but - * WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * General Public License for more details. - */ - -#ifndef RADMETH_OPTIMIZE_HPP -#define RADMETH_OPTIMIZE_HPP - -#include - -struct Regression; - -void -fit_regression_model(Regression &r, std::vector &p_estimates, - double &dispersion_estimate); - -#endif From 6c7a2f90ba407dd50765ca3fd13a8766e73cd7d7 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 19 Oct 2025 17:52:02 -0700 Subject: [PATCH 20/33] src/radmeth/radmeth_optimize_series.hcpp: renamed --- src/radmeth/radmeth_optimize_series.cpp | 619 ++++++++++++++++++++++++ src/radmeth/radmeth_optimize_series.hpp | 29 ++ 2 files changed, 648 insertions(+) create mode 100644 src/radmeth/radmeth_optimize_series.cpp create mode 100644 src/radmeth/radmeth_optimize_series.hpp diff --git a/src/radmeth/radmeth_optimize_series.cpp b/src/radmeth/radmeth_optimize_series.cpp new file mode 100644 index 00000000..e23bd367 --- /dev/null +++ b/src/radmeth/radmeth_optimize_series.cpp @@ -0,0 +1,619 @@ +/* Copyright (C) 2025 Andrew D. Smith + * + * Author: Andrew D. Smith + * + * This program is free software: you can redistribute it and/or modify it + * under the terms the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. + * + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + */ + +#include "radmeth_optimize_series.hpp" + +#include "radmeth_model.hpp" +#include "radmeth_optimize_params.hpp" + +#include +#include + +#include +#include +#include +#include +#include + +#include ///////////////////// +#include ///////////////////// + +[[nodiscard]] static inline auto +logistic(const double x) -> double { + return 1.0 / (1.0 / std::exp(x) + 1.0); +} + +template +[[nodiscard]] static inline auto +get_p(const std::vector &v, const gsl_vector *params) -> double { + const auto a = v.data(); + return logistic(std::inner_product(a, a + std::size(v), params->data, 0.0)); +} + +template +[[nodiscard]] static inline auto +get_p(const std::vector &v, const std::vector ¶ms) -> double { + const auto a = v.data(); + return logistic(std::inner_product(a, a + std::size(v), params.data(), 0.0)); +} + +static inline auto +set_max_r_count(Regression ®) { + const auto &cumul = reg.cumul; + const auto max_itr = std::max_element( + std::cbegin(cumul), std::cend(cumul), [](const auto &a, const auto &b) { + return std::size(a.r_counts) < std::size(b.r_counts); + }); + reg.max_r_count = std::size(max_itr->r_counts); + // ADS: avoid the realloc that can happen even for resize(smaller_size) + if (reg.max_r_count > std::size(reg.cache_log1p_factors)) { + reg.cache_log1p_factors.resize(reg.max_r_count); + reg.cache_dispersion_effect.resize(reg.max_r_count); + } +} + +static inline auto +get_cached_log1p_factors(Regression ®, const double phi) { + const std::size_t max_k = reg.max_r_count; + auto &cache_log1p_factors = reg.cache_log1p_factors; + for (std::size_t k = 0; k < max_k; ++k) + cache_log1p_factors[k] = std::log1p(phi * (k - 1.0)); +} + +static inline auto +get_cached_dispersion_effect(Regression ®, const double phi) { + const std::size_t max_k = reg.max_r_count; + auto &cache_dispersion_effect = reg.cache_dispersion_effect; + double j = -1.0; + const auto lim = std::cbegin(cache_dispersion_effect) + max_k; + for (auto it = std::begin(cache_dispersion_effect); it != lim; ++it, ++j) + *it = j / (1.0 + phi * j); +} + +[[nodiscard]] static double +log_likelihood(const gsl_vector *params, Regression ®) { + const auto phi = logistic(gsl_vector_get(params, reg.design.n_factors())); + const auto one_minus_phi = 1.0 - phi; + + const auto n_groups = reg.n_groups(); + const auto &groups = reg.design.groups; + const auto &cumul = reg.cumul; + + // ADS: precompute log1p(phi * (k - 1.0)), reused for each group + get_cached_log1p_factors(reg, phi); + const auto &log1p_fact_v = reg.cache_log1p_factors; + + double log_lik = 0.0; + for (std::size_t g_idx = 0; g_idx < n_groups; ++g_idx) { + const auto p = get_p(groups[g_idx], params); + const auto one_minus_p = 1.0 - p; + + const auto term1 = one_minus_phi * p; + const auto &cumul_y = cumul[g_idx].m_counts; + std::size_t sz = std::size(cumul_y); + for (std::size_t k = 0; k < sz; ++k) + log_lik += cumul_y[k] * std::log(term1 + phi * k); + + const auto term2 = one_minus_phi * one_minus_p; + const auto &cumul_d = cumul[g_idx].d_counts; + sz = std::size(cumul_d); + for (std::size_t k = 0; k < sz; ++k) + log_lik += cumul_d[k] * std::log(term2 + phi * k); + + const auto &cumul_n = cumul[g_idx].r_counts; + sz = std::size(cumul_n); + for (std::size_t k = 0; k < sz; ++k) + log_lik -= cumul_n[k] * log1p_fact_v[k]; + } + return log_lik; +} + +[[nodiscard]] static double +log_likelihood(const std::vector ¶ms, + Regression ®) { + const auto phi = logistic(params[reg.design.n_factors()]); + const auto one_minus_phi = 1.0 - phi; + + const auto n_groups = reg.n_groups(); + const auto &groups = reg.design.groups; + const auto &cumul = reg.cumul; + + // ADS: precompute log1p(phi * (k - 1.0)), reused for each group + get_cached_log1p_factors(reg, phi); + const auto &log1p_fact_v = reg.cache_log1p_factors; + + double log_lik = 0.0; + for (std::size_t g_idx = 0; g_idx < n_groups; ++g_idx) { + const auto p = get_p(groups[g_idx], params); + const auto one_minus_p = 1.0 - p; + + const auto term1 = one_minus_phi * p; + const auto &cumul_y = cumul[g_idx].m_counts; + std::size_t sz = std::size(cumul_y); + for (std::size_t k = 0; k < sz; ++k) + log_lik += cumul_y[k] * std::log(term1 + phi * k); + + const auto term2 = one_minus_phi * one_minus_p; + const auto &cumul_d = cumul[g_idx].d_counts; + sz = std::size(cumul_d); + for (std::size_t k = 0; k < sz; ++k) + log_lik += cumul_d[k] * std::log(term2 + phi * k); + + const auto &cumul_n = cumul[g_idx].r_counts; + sz = std::size(cumul_n); + for (std::size_t k = 0; k < sz; ++k) + log_lik -= cumul_n[k] * log1p_fact_v[k]; + } + return -log_lik; +} + +static void +gradient(const gsl_vector *params, Regression ®, + gsl_vector *output) { + const auto n_factors = reg.design.n_factors(); + const auto phi = logistic(gsl_vector_get(params, n_factors)); + const auto one_minus_phi = 1.0 - phi; + + const auto n_groups = reg.n_groups(); + const auto &groups = reg.design.groups; + const auto &cumul = reg.cumul; + + auto &p_v = reg.p_v; + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + p_v[g_idx] = get_p(groups[g_idx], params); + + get_cached_dispersion_effect(reg, phi); // (k-1)/(1 + phi(k-1)) + const auto &dispersion_effect = reg.cache_dispersion_effect; + + // init output to zero for all factors + gsl_vector_set_all(output, 0.0); + const auto &data = output->data; + + double disp_deriv = 0.0; + for (std::size_t g_idx = 0; g_idx < n_groups; ++g_idx) { + const auto p = p_v[g_idx]; + const auto one_minus_p = 1.0 - p; + + double deriv = 0.0; + + const auto denom_term1 = one_minus_phi * p; + const auto &cumul_y = cumul[g_idx].m_counts; + const auto y_lim = std::size(cumul_y); + for (auto k = 0u; k < y_lim; ++k) { + const auto common_factor = cumul_y[k] / (denom_term1 + phi * k); + deriv += common_factor; + disp_deriv += (k - p) * common_factor; + } + + const auto denom_term2 = one_minus_phi * one_minus_p; + const auto &cumul_d = cumul[g_idx].d_counts; + const auto d_lim = std::size(cumul_d); + for (auto k = 0u; k < d_lim; ++k) { + const auto common_factor = cumul_d[k] / (denom_term2 + phi * k); + deriv -= common_factor; + disp_deriv += (k - one_minus_p) * common_factor; + } + + const auto &cumul_n = cumul[g_idx].r_counts; + const auto n_lim = std::size(cumul_n); + for (auto k = 0u; k < n_lim; ++k) + disp_deriv -= cumul_n[k] * dispersion_effect[k]; + + const auto &g = groups[g_idx]; + const auto denom_term1_one_minus_p = denom_term1 * one_minus_p; + for (auto fact_idx = 0u; fact_idx < n_factors; ++fact_idx) { + const auto level = g[fact_idx]; + if (level == 0) + continue; + data[fact_idx] += deriv * (denom_term1_one_minus_p * level); + } + } + gsl_vector_set(output, n_factors, disp_deriv * (phi * one_minus_phi)); +} + +static void +gradient(const std::vector ¶ms, Regression ®, + std::vector &output) { + const auto n_factors = reg.design.n_factors(); + const auto phi = logistic(params[n_factors]); + const auto one_minus_phi = 1.0 - phi; + + const auto n_groups = reg.n_groups(); + const auto &groups = reg.design.groups; + const auto &cumul = reg.cumul; + + auto &p_v = reg.p_v; + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + p_v[g_idx] = get_p(groups[g_idx], params); + + get_cached_dispersion_effect(reg, phi); // (k-1)/(1 + phi(k-1)) + const auto &dispersion_effect = reg.cache_dispersion_effect; + + // init output to zero for all factors + output.clear(); + output.resize(std::size(params), 0.0); + auto &data = output; + + double disp_deriv = 0.0; + for (std::size_t g_idx = 0; g_idx < n_groups; ++g_idx) { + const auto p = p_v[g_idx]; + const auto one_minus_p = 1.0 - p; + + // if (p == 0.0 || one_minus_p == 0.0 || phi == 0.0 || one_minus_phi == 0.0) + // { + // std::cout << p << '\t' << one_minus_p << std::endl; + // std::cout << phi << '\t' << one_minus_phi << std::endl; + // } + + double deriv = 0.0; + + const auto denom_term1 = one_minus_phi * p; + const auto &cumul_y = cumul[g_idx].m_counts; + const auto y_lim = std::size(cumul_y); + for (auto k = 0u; k < y_lim; ++k) { + const auto common_factor = cumul_y[k] / (denom_term1 + phi * k); + deriv += common_factor; + disp_deriv += (k - p) * common_factor; + // if (!std::isfinite(disp_deriv)) { + // // clang-format off + // std::cout << "common_factor: " << common_factor << '\n' + // << "cumul_y[k]: " << cumul_y[k] << '\n' + // << "denom_term1: " << denom_term1 << '\n' + // << "phi: " << phi << '\n' + // << "k: " << k << '\n' + // << "p: " << p << '\n' + // << "one_minus_phi: " << one_minus_phi << '\n' + // << "phi: " << phi << '\n' + // << "params[n_factors]: " << params[n_factors] << '\n' + // << "(denom_term1 + phi * k): " << (denom_term1 + phi * k) + // << '\n' + // << std::endl; + // throw std::runtime_error("bad disp_deriv 1"); + // // clang-format on + // } + } + + const auto denom_term2 = one_minus_phi * one_minus_p; + const auto &cumul_d = cumul[g_idx].d_counts; + const auto d_lim = std::size(cumul_d); + for (auto k = 0u; k < d_lim; ++k) { + const auto common_factor = cumul_d[k] / (denom_term2 + phi * k); + deriv -= common_factor; + disp_deriv += (k - one_minus_p) * common_factor; + // if (!std::isfinite(disp_deriv)) + // throw std::runtime_error("bad disp_deriv 2"); + } + + const auto &cumul_n = cumul[g_idx].r_counts; + const auto n_lim = std::size(cumul_n); + for (auto k = 0u; k < n_lim; ++k) { + disp_deriv -= cumul_n[k] * dispersion_effect[k]; + // if (!std::isfinite(disp_deriv)) + // throw std::runtime_error("bad disp_deriv 3"); + } + + const auto &g = groups[g_idx]; + const auto denom_term1_one_minus_p = denom_term1 * one_minus_p; + for (auto fact_idx = 0u; fact_idx < n_factors; ++fact_idx) { + const auto level = g[fact_idx]; + if (level == 0) + continue; + data[fact_idx] += deriv * (denom_term1_one_minus_p * level); + } + } + // if (!std::isfinite(disp_deriv)) + // throw std::runtime_error("bad disp_deriv"); + output[n_factors] = disp_deriv * (phi * one_minus_phi); + // for (auto i = 0u; i <= n_factors; ++i) { + // if (!std::isfinite(output[i])) + // throw std::runtime_error("bad2"); + // } + for (auto i = 0u; i <= n_factors; ++i) { + output[i] *= -1.0; + // if (!std::isfinite(output[i])) + // throw std::runtime_error("bad"); + } + // std::cout << "grad=" << disp_deriv * (phi * one_minus_phi) << std::endl; +} + +[[nodiscard]] static double +neg_loglik(const gsl_vector *params, void *object) { + auto reg = static_cast *>(object); + return -log_likelihood(params, *reg); +} + +static void +neg_gradient(const gsl_vector *params, void *object, gsl_vector *output) { + auto reg = static_cast *>(object); + gradient(params, *reg, output); + gsl_vector_scale(output, -1.0); +} + +static void +neg_loglik_and_grad(const gsl_vector *params, void *object, double *loglik_val, + gsl_vector *d_loglik_val) { + *loglik_val = neg_loglik(params, object); + neg_gradient(params, object, d_loglik_val); +} + +static void +get_cumulative(const std::vector &group_id, + const std::uint32_t n_groups, + const std::vector> &mc, + std::vector &cumul) { + const auto n_cols = std::size(mc); + cumul.clear(); + cumul.resize(n_groups); + + const auto comp_cumul = [&](auto get_value, auto get_vector) { + // phase 1: determine max value for each group + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) { + std::uint32_t max_v{}; + for (auto c_idx = 0u; c_idx < n_cols; ++c_idx) { + if (group_id[c_idx] == g_idx) { + const auto val = get_value(mc[c_idx]); + if (val > max_v) + max_v = val; + } + } + get_vector(cumul[g_idx]).resize(max_v, 0); + } + + // phase 2: fill cumulative counts + for (auto c_idx = 0u; c_idx < n_cols; ++c_idx) { + const auto g_idx = group_id[c_idx]; + const auto val = get_value(mc[c_idx]); + auto &vec = get_vector(cumul[g_idx]); + for (auto i = 0u; i < val; ++i) + ++vec[i]; + } + }; + // call 3 times: m_counts, r_counts, d_counts + comp_cumul( + [](const mcounts &m) { return m.n_meth; }, + [](cumul_counts &c) -> std::vector & { return c.m_counts; }); + + comp_cumul( + [](const mcounts &m) { return m.n_reads; }, + [](cumul_counts &c) -> std::vector & { return c.r_counts; }); + + comp_cumul( + [](const mcounts &m) { return m.n_reads - m.n_meth; }, + [](cumul_counts &c) -> std::vector & { return c.d_counts; }); +} + +static inline auto +conjugate_pr(std::vector ¶ms, Regression ®, + const double alpha) -> int { + const auto n = std::size(params); + const auto max_iter = radmeth_optimize_params::max_iter; + const auto tol = radmeth_optimize_params::tolerance; + + std::vector g(n); + std::vector g_prev(n); + gradient(params, reg, g); + + auto d = g; + std::transform(std::begin(d), std::end(d), std::begin(d), + [](const auto x) { return -1.0 * x; }); + + // const auto ll = log_likelihood(params, reg); + // std::cout << "ll=" << ll << std::endl; + + double gdot{}; + for (const double gi : g) + gdot += gi * gi; + double denom = gdot; + double gnorm = std::sqrt(gdot); + + std::size_t iter{}; + for (iter = 0; iter < max_iter; ++iter) { + // std::memcpy(g_prev.data(), g.data(), sizeof(double) * n); + // update parameters based on directions and step size + for (std::size_t i = 0; i < n; ++i) + params[i] += alpha * d[i]; + + gradient(params, reg, g); + + // get numerator and for norm of gradient + gdot = 0.0; + for (const double gi : g) + gdot += gi * gi; + gnorm = std::sqrt(gdot); + + // stop if norm of gradient is small enough + if (gnorm < tol) + break; + + // // get denom + // double denom = 0.0; + // for (const double gpi : g_prev) + // denom += gpi * gpi; + + // update conjugate directions + const double beta = std::max(0.0, gdot / denom); + for (auto i = 0u; i < n; ++i) + d[i] = -g[i] + beta * d[i]; + + denom = gdot; + } + return iter; + // if (gnorm < tol) { + // std::cout << "Converged in " << iter << " iterations" << "\t" << gnorm + // << '\t' << tol << std::endl; + // } + // else { + // std::cout << "Did not converge in " << max_iter << " iterations\n"; + // } +} + +void +fit_regression_model(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate) { + static constexpr auto init_dispersion_param = -2.5; + const auto stepsize = radmeth_optimize_params::stepsize; + const auto max_iter = radmeth_optimize_params::max_iter; + + const auto n_groups = r.n_groups(); + get_cumulative(r.design.group_id, n_groups, r.mc, r.cumul); + set_max_r_count(r); + + r.p_v.resize(n_groups); + + const std::size_t n_params = r.n_params(); + const auto tol = + std::sqrt(n_params) * r.n_samples() * radmeth_optimize_params::tolerance; + // clang-format off + auto loglik_bundle = gsl_multimin_function_fdf{ + &neg_loglik, // objective function + &neg_gradient, // gradient + &neg_loglik_and_grad, // combined obj and grad + n_params, // number of model params + static_cast(&r) // params for objective and gradients + }; + // clang-format on + + // std::vector x(n_params, 0.0); + // // for (auto i = 0u; i < n_params; ++i) + // // x[i] = gsl_vector_get(param_estimates, i); + // x.back() = -2.5; + // bool all_finite = false; + // auto start_time_a = std::chrono::steady_clock::now(); + // auto ii = 0; + // for (ii = 0; ii < 4 && !all_finite; ++ii) { + // double alpha = std::pow(10.0, -1.0 * (ii + 2)); + // std::cout << "alpha=" << alpha << std::endl; + // std::vector x_out(x); + // conjugate_pr(x_out, r, alpha); + // all_finite = true; + // for (auto i = 0u; i < std::size(x_out) && all_finite; ++i) + // all_finite = std::isfinite(x_out[i]); + // if (all_finite) + // x = x_out; + // } + // // if (!all_finite) { + // // alpha = 1e-03; + // // x_out = x; + // // ii = conjugate_pr(x, r, alpha); + // // all_finite = true; + // // for (auto i = 0u; i < std::size(x) && all_finite; ++i) + // // all_finite = std::isfinite(x[i]); + // // if (!all_finite) { + // // alpha = 1e-04; + // // x_out = x; + // // ii = conjugate_pr(x, r, alpha); + // // all_finite = true; + // // for (auto i = 0u; i < std::size(x) && all_finite; ++i) + // // all_finite = std::isfinite(x[i]); + // // if (!all_finite) { + // // alpha = 1e-05; + // // x_out = x; + // // ii = conjugate_pr(x, r, alpha); + // // all_finite = true; + // // for (auto i = 0u; i < std::size(x) && all_finite; ++i) + // // all_finite = std::isfinite(x[i]); + // // if (!all_finite) { + // // alpha = 1e-06; + // // x_out = x; + // // ii = conjugate_pr(x, r, alpha); + // // all_finite = true; + // // for (auto i = 0u; i < std::size(x) && all_finite; ++i) + // // all_finite = std::isfinite(x[i]); + // // } + // // } + // // } + // // } + // auto stop_time_a = std::chrono::steady_clock::now(); + // // std::cout << "[total time: " << (stop_time_a - start_time_a).count() << + // // "]\t" + // // << ii << std::endl; + + // parameters: 0 for "p" params and final one is a constant for dispersion + auto params = gsl_vector_alloc(n_params); + gsl_vector_set_all(params, 0.0); + gsl_vector_set(params, n_params - 1, init_dispersion_param); + + // Alternatives: + // - gsl_multimin_fdfminimizer_conjugate_pr + // - gsl_multimin_fdfminimizer_conjugate_fr + // - gsl_multimin_fdfminimizer_vector_bfgs2 + // - gsl_multimin_fdfminimizer_steepest_descent + + const auto minimizer = gsl_multimin_fdfminimizer_conjugate_fr; + auto s = gsl_multimin_fdfminimizer_alloc(minimizer, n_params); + gsl_multimin_fdfminimizer_set(s, &loglik_bundle, params, stepsize, tol); + + // gsl_multimin_function minex_func; + + // minex_func.n = n_params; + // minex_func.f = &neg_loglik; + // minex_func.params = static_cast(&r); + + // /* Set initial step sizes to 1 */ + // auto ss = gsl_vector_alloc(n_params); + // gsl_vector_set_all(ss, 1.0); + + // const auto minimizer = gsl_multimin_fminimizer_nmsimplex2; + // auto s = gsl_multimin_fminimizer_alloc(minimizer, n_params); + // gsl_multimin_fminimizer_set(s, &minex_func, params, ss); + + auto start_time_b = std::chrono::steady_clock::now(); + int status = 0; + std::size_t iter = 0; + + double size{}; + do { + status = gsl_multimin_fdfminimizer_iterate(s); // one iter and get + // status + // status = gsl_multimin_fminimizer_iterate(s); // one iter and get status + if (status) + break; + // check status from gradient + + status = gsl_multimin_test_gradient(s->gradient, tol); + + // size = gsl_multimin_fminimizer_size(s); + // status = gsl_multimin_test_size(size, 1e-2); + + } while (status == GSL_CONTINUE && ++iter < max_iter); + // ADS: can't use (status != GSL_SUCCESS) + auto stop_time_b = std::chrono::steady_clock::now(); + // std::cout << "[total time: " + // << (stop_time_b - start_time_b).count() + // // - + // // (stop_time_a - start_time_a).count() + // << "]\t" << iter << std::endl; + + const auto param_estimates = gsl_multimin_fdfminimizer_x(s); + // const auto param_estimates = gsl_multimin_fminimizer_x(s); + // for (auto i = 0u; i < n_params; ++i) + // std::cout << "x[i]=" << gsl_vector_get(param_estimates, i) << std::endl; + + // gsl_vector_fprintf(stdout, param_estimates, "%g"); + + const auto &groups = r.design.groups; + p_estimates.clear(); + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + p_estimates.push_back(get_p(groups[g_idx], param_estimates)); + const auto disp_param = gsl_vector_get(param_estimates, n_params - 1); + dispersion_estimate = 1.0 / std::exp(disp_param); + + r.max_loglik = log_likelihood(param_estimates, r); + + gsl_multimin_fdfminimizer_free(s); + // gsl_multimin_fminimizer_free(s); + gsl_vector_free(params); +} diff --git a/src/radmeth/radmeth_optimize_series.hpp b/src/radmeth/radmeth_optimize_series.hpp new file mode 100644 index 00000000..12b4b312 --- /dev/null +++ b/src/radmeth/radmeth_optimize_series.hpp @@ -0,0 +1,29 @@ +/* Copyright (C) 2025 Andrew D. + * + * Author: Andrew D. Smith + * + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. + * + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + */ + +#ifndef RADMETH_OPTIMIZE_SERIES_HPP +#define RADMETH_OPTIMIZE_SERIES_HPP + +#include +#include + +template struct Regression; + +void +fit_regression_model(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate); + +#endif // RADMETH_OPTIMIZE_SERIES_HPP From 7b69ce5658a5d8612417905eca474acbeef7c889 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 19 Oct 2025 17:52:35 -0700 Subject: [PATCH 21/33] src/radmeth/radmeth_optimize_gamma.hcpp: renamed and adding non-gradient optimization choice for experiments --- src/radmeth/radmeth_optimize_gamma.cpp | 422 +++++++++++++++++++++++++ src/radmeth/radmeth_optimize_gamma.hpp | 44 +++ 2 files changed, 466 insertions(+) create mode 100644 src/radmeth/radmeth_optimize_gamma.cpp create mode 100644 src/radmeth/radmeth_optimize_gamma.hpp diff --git a/src/radmeth/radmeth_optimize_gamma.cpp b/src/radmeth/radmeth_optimize_gamma.cpp new file mode 100644 index 00000000..592f5fb0 --- /dev/null +++ b/src/radmeth/radmeth_optimize_gamma.cpp @@ -0,0 +1,422 @@ +/* Copyright (C) 2013-2025 Andrew D Smith + * + * Author: Andrew D. Smith + * + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. + * + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + */ + +#include "radmeth_model.hpp" +#include "radmeth_optimize_gamma.hpp" +#include "radmeth_optimize_params.hpp" + +#include +#include + +#include +#include +#include +#include +#include + +/* Coefficients for the Chebyschev polynomial for the digamma function in the + range 0-1 */ +// clang-format off +static constexpr std::array psi1_cs { + -0.038057080835217922, // == -0.019028540417608961*2 + 0.491415393029387130, + -0.056815747821244730, + 0.008357821225914313, + -0.001333232857994342, + 0.000220313287069308, + -0.000037040238178456, + 0.000006283793654854, + -0.000001071263908506, + 0.000000183128394654, + -0.000000031353509361, + 0.000000005372808776, + -0.000000000921168141, + 0.000000000157981265, + -0.000000000027098646, + 0.000000000004648722, + -0.000000000000797527, + 0.000000000000136827, + -0.000000000000023475, + 0.000000000000004027, + -0.000000000000000691, + 0.000000000000000118, + -0.000000000000000020 +}; +// clang-format on + +/* Alternate set of coefficients for the Chebyschev polynomial for the digamma + function */ +// clang-format off +static constexpr std::array psi2_cs = { + -0.0204749044678185, // == -0.01023745223390925*2 + -0.0101801271534859, + 0.0000559718725387, + -0.0000012917176570, + 0.0000000572858606, + -0.0000000038213539, + 0.0000000003397434, + -0.0000000000374838, + 0.0000000000048990, + -0.0000000000007344, + 0.0000000000001233, + -0.0000000000000228, + 0.0000000000000045, + -0.0000000000000009, + 0.0000000000000002, + -0.0000000000000000 +}; +// clang-format on + +template +[[nodiscard]] static inline double +chebyschev(const std::array &coeffs, std::uint32_t order, + const double y) { + const auto y2 = 2.0 * y; + double d = 0.0; + double dd = 0.0; + for (auto j = order; j >= 1; j--) { + const auto temp = d; + d = y2 * d - dd + coeffs[j]; + dd = temp; + } + return y * d - dd + 0.5 * coeffs[0]; +} + +// digamma for x non-negative +static inline auto +digamma(const double y) -> double { + static constexpr auto psi1_order = 22; // max=22; + static constexpr auto psi2_order = 15; // max=15; + if (y >= 2.0) { + const auto t = 8.0 / (y * y) - 1.0; + return std::log(y) - 0.5 / y + chebyschev(psi2_cs, psi2_order, t); + } + if (y < 1.0) + return -1.0 / y + chebyschev(psi1_cs, psi1_order, 2.0 * y - 1.0); + return chebyschev(psi1_cs, psi1_order, 2.0 * (y - 1.0) - 1.0); +} + +[[nodiscard]] static inline auto +logistic(const double x) -> double { + return 1.0 / (1.0 / std::exp(x) + 1.0); +} + +template +[[nodiscard]] static auto +get_p(const std::vector &v, const gsl_vector *params) -> double { + const auto a = v.data(); + return logistic(std::inner_product(a, a + std::size(v), params->data, 0.0)); +} + +static auto +get_cache_lgamma(const std::vector &group, + const gsl_vector *params, const double phi, + vars_cache &cache) { + const auto p = get_p(group, params); + const auto a = p * phi; + const auto b = (1.0 - p) * phi; + cache.p = p; + cache.a = a; + cache.b = b; + cache.lgamma_a = std::lgamma(a); + cache.lgamma_b = std::lgamma(b); + cache.lgamma_a_b = std::lgamma(a + b); +} + +template +[[nodiscard]] static auto +log_likelihood(const gsl_vector *params, Regression ®) -> double { + const auto phi = 1.0 / std::exp(gsl_vector_get(params, reg.n_factors())); + + const auto n_groups = reg.n_groups(); + const auto &groups = reg.design.groups; + auto &cache = reg.cache; + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + get_cache_lgamma(groups[g_idx], params, phi, cache[g_idx]); + + const auto &group_id = reg.design.group_id; + const auto &mc = reg.mc; + double ll = 0.0; + const auto n_samples = reg.n_samples(); + for (auto i = 0u; i < n_samples; ++i) { + const auto y = mc[i].n_meth; + const auto n = mc[i].n_reads; + + const auto &c = cache[group_id[i]]; + const auto a = c.a; + const auto b = c.b; + + // clang-format off + ll += ((std::lgamma(y + a) - c.lgamma_a) + + (std::lgamma(n - y + b) - c.lgamma_b) + + (c.lgamma_a_b - std::lgamma(n + a + b))); + // clang-format on + } + + return ll; +} + +static void +get_cache_digamma(const std::vector &group, + const gsl_vector *params, const double phi, + vars_cache &cache) { + const auto p = get_p(group, params); + const auto a = p * phi; + const auto b = (1.0 - p) * phi; + cache.p = p; + cache.a = a; + cache.b = b; + cache.digamma_a = digamma(a); + cache.digamma_b = digamma(b); + cache.digamma_a_b = digamma(a + b); +} + +template +static void +gradient(const gsl_vector *params, Regression ®, gsl_vector *output) { + const auto n_factors = reg.n_factors(); + const auto phi = 1.0 / std::exp(gsl_vector_get(params, n_factors)); + + const auto &mc = reg.mc; + const auto &matrix = reg.design.matrix; + + const auto n_groups = reg.n_groups(); + const auto &groups = reg.design.groups; + auto &cache = reg.cache; // ADS: reusing scratch space + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + get_cache_digamma(groups[g_idx], params, phi, cache[g_idx]); + + gsl_vector_set_all(output, 0.0); // init output to zero for all factors + const auto &data = output->data; + + const auto &group_id = reg.design.group_id; + double grad_phi = 0.0; + const auto n_samples = reg.n_samples(); + for (auto i = 0u; i < n_samples; ++i) { + const double y = mc[i].n_meth; // ADS: n_meth maybe int or float + const auto n = mc[i].n_reads; + + const auto &c = cache[group_id[i]]; + const auto p = c.p; + const auto a = c.a; + const auto b = c.b; + + const auto digamma_a = c.digamma_a; + const auto digamma_b = c.digamma_b; + const auto digamma_y_a = digamma(y + a); + const auto digamma_n_y_b = digamma(n - y + b); + + // grad wrt p + const auto dlogl_dp = + phi * (digamma_y_a - digamma_a - digamma_n_y_b + digamma_b); + const auto dp_delta = p * (1.0 - p); // chain rule: d p / d param + const auto dlogl_delta = dlogl_dp * dp_delta; // grad wrt params + + auto matrix_itr = std::cbegin(matrix[i]); + const auto data_end = data + n_factors; + for (auto data_itr = data; data_itr != data_end; ++data_itr) + *data_itr += (*matrix_itr++) * dlogl_delta; + + const auto digamma_delta = c.digamma_a_b - digamma(n + a + b); + const auto dphi_term1 = p * (digamma_y_a - digamma_a + digamma_delta); + const auto dphi_term2 = + (1.0 - p) * (digamma_n_y_b - digamma_b + digamma_delta); + + grad_phi += dphi_term1 + dphi_term2; + } + + const auto grad_theta = -grad_phi * phi; + gsl_vector_set(output, n_factors, grad_theta); +} + +template +[[nodiscard]] static double +neg_loglik(const gsl_vector *params, void *object) { + auto reg = static_cast *>(object); + return -log_likelihood(params, *reg); +} + +template +static void +neg_gradient(const gsl_vector *params, void *object, gsl_vector *output) { + auto reg = static_cast *>(object); + gradient(params, *reg, output); + gsl_vector_scale(output, -1.0); +} + +template +static void +neg_loglik_and_grad(const gsl_vector *params, void *object, double *loglik_val, + gsl_vector *d_loglik_val) { + *loglik_val = neg_loglik(params, object); + neg_gradient(params, object, d_loglik_val); +} + +template +static void +fit_regression_model(Regression &r, std::vector &p_estimates, + double &dispersion_estimate) { + static constexpr auto init_dispersion_param = -1.0; + const auto stepsize = radmeth_optimize_params::stepsize; + const auto max_iter = radmeth_optimize_params::max_iter; + const auto n_groups = r.n_groups(); + r.cache.resize(n_groups); // make sure scratch space is allocated + + const std::size_t n_params = r.n_params(); + const auto tol = + std::sqrt(n_params) * r.n_samples() * radmeth_optimize_params::tolerance; + + // set the parameters: zero for "p" parameters and the final one for + // dispersion using the constant + auto params = gsl_vector_alloc(n_params); + gsl_vector_set_all(params, 0.0); + gsl_vector_set(params, n_params - 1, init_dispersion_param); + + auto minex_func = gsl_multimin_function{ + &neg_loglik, + n_params, + static_cast(&r), + }; + + // set initial step size for all dims + auto step_sizes = gsl_vector_alloc(n_params); + gsl_vector_set_all(step_sizes, stepsize); + + const auto minimizer = gsl_multimin_fminimizer_nmsimplex2; + auto s = gsl_multimin_fminimizer_alloc(minimizer, n_params); + gsl_multimin_fminimizer_set(s, &minex_func, params, step_sizes); + + int status = 0; + std::size_t iter = 0; + double size{}; + + do { + status = gsl_multimin_fminimizer_iterate(s); // one iter and get status + if (status) + break; + + size = gsl_multimin_fminimizer_size(s); + status = gsl_multimin_test_size(size, tol); + } while (status == GSL_CONTINUE && ++iter < max_iter); + // ADS: can't use (status != GSL_SUCCESS) + + const auto param_estimates = gsl_multimin_fminimizer_x(s); + + const auto &groups = r.design.groups; + p_estimates.clear(); + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + p_estimates.push_back(get_p(groups[g_idx], param_estimates)); + const auto disp_param = gsl_vector_get(param_estimates, n_params - 1); + dispersion_estimate = 1.0 / std::exp(disp_param); + + r.max_loglik = log_likelihood(param_estimates, r); + + gsl_multimin_fminimizer_free(s); + gsl_vector_free(step_sizes); + gsl_vector_free(params); +} + +template +static void +fit_regression_model_fdf(Regression &r, std::vector &p_estimates, + double &dispersion_estimate) { + static constexpr auto init_dispersion_param = -1.0; + const auto stepsize = radmeth_optimize_params::stepsize; + const auto max_iter = radmeth_optimize_params::max_iter; + const auto n_groups = r.n_groups(); + r.cache.resize(n_groups); // make sure scratch space is allocated + + const std::size_t n_params = r.n_params(); + const auto tol = + std::sqrt(n_params) * r.n_samples() * radmeth_optimize_params::tolerance; + // clang-format off + auto loglik_bundle = gsl_multimin_function_fdf{ + &neg_loglik, // objective function + &neg_gradient, // gradient + &neg_loglik_and_grad, // combined obj and grad + n_params, // number of model params + static_cast(&r) // parameters for objective and gradient functions + }; + // clang-format on + + // set the parameters: zero for "p" parameters and the final one for + // dispersion using the constant + auto params = gsl_vector_alloc(n_params); + gsl_vector_set_all(params, 0.0); + gsl_vector_set(params, n_params - 1, init_dispersion_param); + + // Alternatives: + // - gsl_multimin_fdfminimizer_conjugate_pr + // - gsl_multimin_fdfminimizer_conjugate_fr + // - gsl_multimin_fdfminimizer_vector_bfgs2 + // - gsl_multimin_fdfminimizer_steepest_descent + const auto minimizer = gsl_multimin_fdfminimizer_conjugate_pr; + auto s = gsl_multimin_fdfminimizer_alloc(minimizer, n_params); + + gsl_multimin_fdfminimizer_set(s, &loglik_bundle, params, stepsize, tol); + + int status = 0; + std::size_t iter = 0; + + do { + status = gsl_multimin_fdfminimizer_iterate(s); // one iter and get status + if (status) + break; + // check status from gradient + status = gsl_multimin_test_gradient(s->gradient, tol); + } while (status == GSL_CONTINUE && ++iter < max_iter); + // ADS: can't use (status != GSL_SUCCESS) + + const auto param_estimates = gsl_multimin_fdfminimizer_x(s); + const auto &groups = r.design.groups; + p_estimates.clear(); + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + p_estimates.push_back(get_p(groups[g_idx], param_estimates)); + const auto disp_param = gsl_vector_get(param_estimates, n_params - 1); + dispersion_estimate = 1.0 / std::exp(disp_param); + + r.max_loglik = log_likelihood(param_estimates, r); + + gsl_multimin_fdfminimizer_free(s); + gsl_vector_free(params); +} + +void +fit_regression_model_gamma(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate) { + fit_regression_model(r, p_estimates, dispersion_estimate); +} + +void +fit_regression_model_gamma(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate) { + fit_regression_model(r, p_estimates, dispersion_estimate); +} + +void +fit_regression_model_gamma_fdf(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate) { + fit_regression_model_fdf(r, p_estimates, dispersion_estimate); +} + +void +fit_regression_model_gamma_fdf(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate) { + fit_regression_model_fdf(r, p_estimates, dispersion_estimate); +} diff --git a/src/radmeth/radmeth_optimize_gamma.hpp b/src/radmeth/radmeth_optimize_gamma.hpp new file mode 100644 index 00000000..1949b6fb --- /dev/null +++ b/src/radmeth/radmeth_optimize_gamma.hpp @@ -0,0 +1,44 @@ +/* Copyright (C) 2025 Andrew D Smith + * + * Author: Andrew D. Smith + * + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. + * + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + */ + +#ifndef RADMETH_OPTIMIZE_GAMMA_HPP +#define RADMETH_OPTIMIZE_GAMMA_HPP + +#include +#include + +template struct Regression; + +void +fit_regression_model_gamma(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate); + +void +fit_regression_model_gamma(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate); + +void +fit_regression_model_gamma_fdf(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate); + +void +fit_regression_model_gamma_fdf(Regression &r, + std::vector &p_estimates, + double &dispersion_estimate); + +#endif // RADMETH_OPTIMIZE_GAMMA_HPP From 3be1c88451dba7c70766abd931056f9dd33f83e4 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 19 Oct 2025 17:52:51 -0700 Subject: [PATCH 22/33] src/radmeth/radmeth_optimize_nano.hcpp: removed during renaming --- src/radmeth/radmeth_optimize_nano.cpp | 330 -------------------------- src/radmeth/radmeth_optimize_nano.hpp | 27 --- 2 files changed, 357 deletions(-) delete mode 100644 src/radmeth/radmeth_optimize_nano.cpp delete mode 100644 src/radmeth/radmeth_optimize_nano.hpp diff --git a/src/radmeth/radmeth_optimize_nano.cpp b/src/radmeth/radmeth_optimize_nano.cpp deleted file mode 100644 index e9f5a48c..00000000 --- a/src/radmeth/radmeth_optimize_nano.cpp +++ /dev/null @@ -1,330 +0,0 @@ -/* Copyright (C) 2013-2025 Andrew D Smith - * - * Author: Andrew D. Smith - * - * This program is free software: you can redistribute it and/or modify it - * under the terms of the GNU General Public License as published by the Free - * Software Foundation, either version 3 of the License, or (at your option) - * any later version. - * - * This program is distributed in the hope that it will be useful, but WITHOUT - * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for - * more details. - */ - -#include "radmeth_optimize_nano.hpp" -#include "radmeth_model_nano.hpp" - -#include -#include // for gsl_sf_psi (digamma) -#include - -#include -#include -#include -#include -#include - -/* Coefficients for the Chebyschev polynomial for the digamma function in the - range 0-1 */ -// clang-format off -static constexpr std::array psi1_cs { - -0.038057080835217922, // == -0.019028540417608961*2 - 0.491415393029387130, - -0.056815747821244730, - 0.008357821225914313, - -0.001333232857994342, - 0.000220313287069308, - -0.000037040238178456, - 0.000006283793654854, - -0.000001071263908506, - 0.000000183128394654, - -0.000000031353509361, - 0.000000005372808776, - -0.000000000921168141, - 0.000000000157981265, - -0.000000000027098646, - 0.000000000004648722, - -0.000000000000797527, - 0.000000000000136827, - -0.000000000000023475, - 0.000000000000004027, - -0.000000000000000691, - 0.000000000000000118, - -0.000000000000000020 -}; -// clang-format on - -/* Alternate set of coefficients for the Chebyschev polynomial for the digamma - function */ -// clang-format off -static constexpr std::array psi2_cs = { - -0.0204749044678185, // == -0.01023745223390925*2 - -0.0101801271534859, - 0.0000559718725387, - -0.0000012917176570, - 0.0000000572858606, - -0.0000000038213539, - 0.0000000003397434, - -0.0000000000374838, - 0.0000000000048990, - -0.0000000000007344, - 0.0000000000001233, - -0.0000000000000228, - 0.0000000000000045, - -0.0000000000000009, - 0.0000000000000002, - -0.0000000000000000 -}; -// clang-format on - -template -[[nodiscard]] static inline double -chebyschev(const std::array &coeffs, std::uint32_t order, - const double y) { - const auto y2 = 2.0 * y; - double d = 0.0; - double dd = 0.0; - for (auto j = order; j >= 1; j--) { - const auto temp = d; - d = y2 * d - dd + coeffs[j]; - dd = temp; - } - return y * d - dd + 0.5 * coeffs[0]; -} - -[[nodiscard]] static inline double -logistic(const double x) { - return 1.0 / (1.0 / std::exp(x) + 1.0); -} - -template -[[nodiscard]] static double -get_p(const std::vector &v, const gsl_vector *params) { - const auto a = v.data(); - return logistic(std::inner_product(a, a + std::size(v), params->data, 0.0)); -} - -static void -get_cache_lgamma(const std::vector &group, - const gsl_vector *params, const double phi, - vars_cache &cache) { - const auto p = get_p(group, params); - const auto a = p * phi; - const auto b = (1.0 - p) * phi; - cache.p = p; - cache.a = a; - cache.b = b; - cache.lgamma_a = std::lgamma(a); - cache.lgamma_b = std::lgamma(b); - cache.lgamma_a_b = std::lgamma(a + b); -} - -[[nodiscard]] static double -log_likelihood_nano(const gsl_vector *params, RegressionNano ®) { - const auto phi = 1.0 / std::exp(gsl_vector_get(params, reg.n_factors())); - - const auto n_groups = reg.n_groups(); - const auto &groups = reg.design.groups; - auto &cache = reg.cache; - for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) - get_cache_lgamma(groups[g_idx], params, phi, cache[g_idx]); - - const auto &group_id = reg.design.group_id; - const auto &mc = reg.mc; - double ll = 0.0; - const auto n_samples = reg.n_samples(); - for (auto i = 0u; i < n_samples; ++i) { - const auto y = mc[i].n_meth; - const auto n = mc[i].n_reads; - - const auto &c = cache[group_id[i]]; - const auto a = c.a; - const auto b = c.b; - - // clang-format off - ll += ((std::lgamma(y + a) - c.lgamma_a) + - (std::lgamma(n - y + b) - c.lgamma_b) + - (c.lgamma_a_b - std::lgamma(n + a + b))); - // clang-format on - } - - return ll; -} - -// digamma for x non-negative -static double -digamma(const double y) { - static constexpr auto psi1_order = 22; // max=22; - static constexpr auto psi2_order = 15; // max=15; - if (y >= 2.0) { - const auto t = 8.0 / (y * y) - 1.0; - return std::log(y) - 0.5 / y + chebyschev(psi2_cs, psi2_order, t); - } - if (y < 1.0) - return -1.0 / y + chebyschev(psi1_cs, psi1_order, 2.0 * y - 1.0); - return chebyschev(psi1_cs, psi1_order, 2.0 * (y - 1.0) - 1.0); -} - -static void -get_cache_digamma(const std::vector &group, - const gsl_vector *params, const double phi, - vars_cache &cache) { - const auto p = get_p(group, params); - const auto a = p * phi; - const auto b = (1.0 - p) * phi; - cache.p = p; - cache.a = a; - cache.b = b; - cache.digamma_a = digamma(a); - cache.digamma_b = digamma(b); - cache.digamma_a_b = digamma(a + b); -} - -static void -gradient_nano(const gsl_vector *params, RegressionNano ®, - gsl_vector *output) { - const auto n_factors = reg.n_factors(); - const auto phi = 1.0 / std::exp(gsl_vector_get(params, n_factors)); - - const auto &mc = reg.mc; - const auto &matrix = reg.design.matrix; - - const auto n_groups = reg.n_groups(); - const auto &groups = reg.design.groups; - auto &cache = reg.cache; // ADS: reusing scratch space - for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) - get_cache_digamma(groups[g_idx], params, phi, cache[g_idx]); - - gsl_vector_set_all(output, 0.0); // init output to zero for all factors - const auto &data = output->data; - - const auto &group_id = reg.design.group_id; - double grad_phi = 0.0; - const auto n_samples = reg.n_samples(); - for (auto i = 0u; i < n_samples; ++i) { - const auto y = mc[i].n_meth; - const auto n = mc[i].n_reads; - - const auto &c = cache[group_id[i]]; - const auto p = c.p; - const auto a = c.a; - const auto b = c.b; - - const auto digamma_a = c.digamma_a; - const auto digamma_b = c.digamma_b; - const auto digamma_y_a = digamma(y + a); - const auto digamma_n_y_b = digamma(n - y + b); - - // grad wrt p - const auto dlogl_dp = - phi * (digamma_y_a - digamma_a - digamma_n_y_b + digamma_b); - const auto dp_delta = p * (1.0 - p); // chain rule: d p / d param - const auto dlogl_delta = dlogl_dp * dp_delta; // grad wrt params - - auto matrix_itr = std::cbegin(matrix[i]); - const auto data_end = data + n_factors; - for (auto data_itr = data; data_itr != data_end; ++data_itr) - *data_itr += (*matrix_itr++) * dlogl_delta; - - const auto digamma_delta = c.digamma_a_b - digamma(n + a + b); - const auto dphi_term1 = p * (digamma_y_a - digamma_a + digamma_delta); - const auto dphi_term2 = - (1.0 - p) * (digamma_n_y_b - digamma_b + digamma_delta); - - grad_phi += dphi_term1 + dphi_term2; - } - - const auto grad_theta = -grad_phi * phi; - gsl_vector_set(output, n_factors, grad_theta); -} - -[[nodiscard]] static double -neg_loglik_nano(const gsl_vector *params, void *object) { - auto reg = static_cast(object); - return -log_likelihood_nano(params, *reg); -} - -static void -neg_gradient_nano(const gsl_vector *params, void *object, gsl_vector *output) { - auto reg = static_cast(object); - gradient_nano(params, *reg, output); - gsl_vector_scale(output, -1.0); -} - -static void -neg_loglik_and_grad_nano(const gsl_vector *params, void *object, - double *loglik_val, gsl_vector *d_loglik_val) { - *loglik_val = neg_loglik_nano(params, object); - neg_gradient_nano(params, object, d_loglik_val); -} - -void -fit_regression_model_nano(RegressionNano &r, std::vector &p_estimates, - double &dispersion_estimate) { - static constexpr auto init_dispersion_param = -2.5; - const auto stepsize = RegressionNano::stepsize; - const auto max_iter = RegressionNano::max_iter; - const auto n_groups = r.n_groups(); - r.cache.resize(n_groups); // make sure scratch space is allocated - - const std::size_t n_params = r.n_params(); - const auto tol = - std::sqrt(n_params) * r.n_samples() * RegressionNano::tolerance; - // clang-format off - auto loglik_bundle = gsl_multimin_function_fdf{ - &neg_loglik_nano, // objective function - &neg_gradient_nano, // gradient - &neg_loglik_and_grad_nano, // combined obj and grad - n_params, // number of model params - static_cast(&r) // parameters for objective and gradient functions - }; - // clang-format on - - // set the parameters: zero for "p" parameters and the final one for - // dispersion using the constant - auto params = gsl_vector_alloc(n_params); - gsl_vector_set_all(params, 0.0); - gsl_vector_set(params, n_params - 1, init_dispersion_param); - - // Alternatives: - // - gsl_multimin_fdfminimizer_conjugate_pr - // - gsl_multimin_fdfminimizer_conjugate_fr - // - gsl_multimin_fdfminimizer_vector_bfgs2 - // - gsl_multimin_fdfminimizer_steepest_descent - const auto minimizer = gsl_multimin_fdfminimizer_conjugate_pr; - auto s = gsl_multimin_fdfminimizer_alloc(minimizer, n_params); - - gsl_multimin_fdfminimizer_set(s, &loglik_bundle, params, stepsize, tol); - - int status = 0; - std::size_t iter = 0; - do { - status = gsl_multimin_fdfminimizer_iterate(s); // one iter and get status - if (status) - break; - // check status from gradient - status = gsl_multimin_test_gradient(s->gradient, tol); - } while (status == GSL_CONTINUE && ++iter < max_iter); - - /// ADS: the condition below might not always pass even if we are doing - /// ok. It's not clear how to check for failure. - - // if (status != GSL_SUCCESS) - // throw std::runtime_error("failed to fit model parameters"); - - const auto param_estimates = gsl_multimin_fdfminimizer_x(s); - - const auto &groups = r.design.groups; - p_estimates.clear(); - for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) - p_estimates.push_back(get_p(groups[g_idx], param_estimates)); - const auto disp_param = gsl_vector_get(param_estimates, n_params - 1); - dispersion_estimate = 1.0 / std::exp(disp_param); - - r.max_loglik = log_likelihood_nano(param_estimates, r); - - gsl_multimin_fdfminimizer_free(s); - gsl_vector_free(params); -} diff --git a/src/radmeth/radmeth_optimize_nano.hpp b/src/radmeth/radmeth_optimize_nano.hpp deleted file mode 100644 index 3487cbea..00000000 --- a/src/radmeth/radmeth_optimize_nano.hpp +++ /dev/null @@ -1,27 +0,0 @@ -/* Copyright (C) 2025 Andrew D Smith - * - * Author: Andrew D. Smith - * - * This program is free software: you can redistribute it and/or modify it - * under the terms of the GNU General Public License as published by the Free - * Software Foundation, either version 3 of the License, or (at your option) - * any later version. - * - * This program is distributed in the hope that it will be useful, but WITHOUT - * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for - * more details. - */ - -#ifndef RADMETH_OPTIMIZE_NANO_HPP -#define RADMETH_OPTIMIZE_NANO_HPP - -#include - -struct RegressionNano; - -void -fit_regression_model_nano(RegressionNano &r, std::vector &p_estimates, - double &dispersion_estimate); - -#endif // RADMETH_OPTIMIZE_NANO_HPP From 65a86ea9a701e0769b3463fffeedbdd34bf42b33 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 19 Oct 2025 17:53:34 -0700 Subject: [PATCH 23/33] src/radmeth/radmeth_optimize_params.hpp: moving params in namespace --- src/radmeth/radmeth_optimize_params.hpp | 27 +++++++++++++++++++++++++ 1 file changed, 27 insertions(+) create mode 100644 src/radmeth/radmeth_optimize_params.hpp diff --git a/src/radmeth/radmeth_optimize_params.hpp b/src/radmeth/radmeth_optimize_params.hpp new file mode 100644 index 00000000..a9cbd64a --- /dev/null +++ b/src/radmeth/radmeth_optimize_params.hpp @@ -0,0 +1,27 @@ +/* Copyright (C) 2025 Andrew D Smith + * + * Author: Andrew D. Smith + * + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. + * + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. + */ + +#ifndef RADMETH_OPTIMIZE_PARAMS_HPP +#define RADMETH_OPTIMIZE_PARAMS_HPP + +#include + +namespace radmeth_optimize_params { +inline double tolerance = 1e-4; +inline double stepsize = 0.01; +inline std::uint32_t max_iter = 250; +}; // namespace radmeth_optimize_params + +#endif // RADMETH_OPTIMIZE_PARAMS_HPP From b27869b52c626b12cdbe4953cf9190498dfd1160 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 19 Oct 2025 17:55:01 -0700 Subject: [PATCH 24/33] src/radmeth/radmeth_model.hpp: collapsing the Regression model classes and templating the counts for integer vs. fractional and putting all the cache structures into one class --- src/radmeth/radmeth_model.hpp | 74 ++++++++++++++++++++++++++++++++--- 1 file changed, 68 insertions(+), 6 deletions(-) diff --git a/src/radmeth/radmeth_model.hpp b/src/radmeth/radmeth_model.hpp index 0be0099d..a202664e 100644 --- a/src/radmeth/radmeth_model.hpp +++ b/src/radmeth/radmeth_model.hpp @@ -18,19 +18,22 @@ #include "radmeth_design.hpp" +#include +#include #include #include #include #include #include -struct mcounts { +template struct mcounts { std::uint32_t n_reads{}; - std::uint32_t n_meth{}; + T n_meth{}; }; +template [[nodiscard]] inline std::istream & -operator>>(std::istream &in, mcounts &rm) { +operator>>(std::istream &in, mcounts &rm) { return in >> rm.n_reads >> rm.n_meth; } @@ -40,20 +43,34 @@ struct cumul_counts { std::vector d_counts; }; -struct Regression { +struct vars_cache { + double p{}; + double a{}; + double b{}; + double lgamma_a{}; + double lgamma_b{}; + double lgamma_a_b{}; + double digamma_a{}; + double digamma_b{}; + double digamma_a_b{}; +}; + +template struct Regression { static double tolerance; // 1e-3; static double stepsize; // 0.001; static std::uint32_t max_iter; // 250; Design design; std::string rowname; - std::vector mc; + std::vector> mc; double max_loglik{}; // scratch space std::vector cumul; std::vector p_v; - std::vector cache; + std::vector cache_log1p_factors; + std::vector cache_dispersion_effect; + std::vector cache; // scratch space std::uint32_t max_r_count{}; void @@ -85,4 +102,49 @@ struct Regression { } }; +template +void +Regression::parse(const std::string &line) { + const auto first_ws = line.find_first_of(" \t"); + + // Parse the row name (must be like: "chr:position:strand:context") + bool failed = (first_ws == std::string::npos); + + auto field_s = line.data(); + auto field_e = line.data() + first_ws; + rowname = std::string{field_s, field_e}; + std::replace(std::begin(rowname), std::end(rowname), ':', '\t'); + if (failed) + throw std::runtime_error("failed to parse label from:\n" + line); + + // Parse the counts of total reads and methylated reads + const auto is_sep = [](const auto x) { return std::isspace(x); }; + const auto not_sep = [](const auto x) { return std::isdigit(x); }; + + const auto line_end = line.data() + std::size(line); + mcounts mc1{}; + mc.clear(); + while (field_e != line_end) { + // get the total count + field_s = std::find_if(field_e + 1, line_end, not_sep); + field_e = std::find_if(field_s + 1, line_end, is_sep); + { + const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_reads); + failed = failed || (ptr != field_e); + } + + field_s = std::find_if(field_e + 1, line_end, not_sep); + field_e = std::find_if(field_s + 1, line_end, is_sep); + { + const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_meth); + failed = failed || (ptr != field_e); + } + + mc.push_back(mc1); + } + + if (failed) + throw std::runtime_error("failed to parse counts from:\n" + line); +} + #endif // RADMETH_MODEL_HPP From 086f58a8bd3079a64132725401c4370ebb595219 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 19 Oct 2025 17:55:45 -0700 Subject: [PATCH 25/33] src/radmeth/radmeth.cpp: update to use template for Regression class --- src/radmeth/radmeth.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/radmeth/radmeth.cpp b/src/radmeth/radmeth.cpp index c5642742..0da6bfc7 100644 --- a/src/radmeth/radmeth.cpp +++ b/src/radmeth/radmeth.cpp @@ -16,7 +16,8 @@ #include "radmeth_design.hpp" #include "radmeth_model.hpp" -#include "radmeth_optimize.hpp" +#include "radmeth_optimize_params.hpp" +#include "radmeth_optimize_series.hpp" #include "radmeth_utils.hpp" // smithlab headers @@ -286,10 +287,10 @@ main_radmeth(int argc, char *argv[]) { opt_parse.add_opt("factor", 'f', "a factor to test", true, test_factor); opt_parse.add_opt("tolerance", '\0', "numerical tolerance to test convergence", false, - Regression::tolerance); + radmeth_optimize_params::tolerance); opt_parse.add_opt("max-iter", '\0', "max iterations when estimating parameters", false, - Regression::max_iter); + radmeth_optimize_params::max_iter); std::vector leftover_args; opt_parse.parse(argc, argv, leftover_args); if (argc == 1 || opt_parse.help_requested()) { @@ -353,9 +354,9 @@ main_radmeth(int argc, char *argv[]) { << '\n'; // clang-format on - Regression alt_model; + Regression alt_model; alt_model.design = design; - Regression null_model; + Regression null_model; null_model.design = null_design; const auto start_time = std::chrono::steady_clock::now(); From 8e5d4643bb5df6d9d13e558f1fdc192fa8a71090 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 19 Oct 2025 17:56:11 -0700 Subject: [PATCH 26/33] src/radmeth/radmeth_nano.cpp: update to use template on Regression class instead of RegressionNano which was removed --- src/radmeth/radmeth_nano.cpp | 45 ++++++++++++++++++------------------ 1 file changed, 22 insertions(+), 23 deletions(-) diff --git a/src/radmeth/radmeth_nano.cpp b/src/radmeth/radmeth_nano.cpp index 7451e2dd..4bb41a73 100644 --- a/src/radmeth/radmeth_nano.cpp +++ b/src/radmeth/radmeth_nano.cpp @@ -14,8 +14,9 @@ */ #include "radmeth_design.hpp" -#include "radmeth_model_nano.hpp" -#include "radmeth_optimize_nano.hpp" +#include "radmeth_model.hpp" +#include "radmeth_optimize_gamma.hpp" +#include "radmeth_optimize_params.hpp" #include "radmeth_utils.hpp" // smithlab headers @@ -36,9 +37,8 @@ #include #include -template [[nodiscard]] static bool -has_low_coverage(const RegressionType ®, const std::size_t test_factor) { +has_low_coverage(const Regression ®, const std::size_t test_factor) { bool cvrd_in_test_fact_smpls = false; const auto &tcol = reg.design.tmatrix[test_factor]; for (std::size_t i = 0; i < reg.n_samples() && !cvrd_in_test_fact_smpls; ++i) @@ -51,9 +51,8 @@ has_low_coverage(const RegressionType ®, const std::size_t test_factor) { return !cvrd_in_test_fact_smpls || !cvrd_in_other_smpls; } -template [[nodiscard]] static bool -has_extreme_counts(const RegressionType ®) { +has_extreme_counts(const Regression ®) { const auto &mc = reg.mc; bool full_meth = true; @@ -74,13 +73,12 @@ enum class row_status : std::uint8_t { na_extreme_cnt, }; -template static void -radmeth_nano(const bool show_progress, const bool more_na_info, - const std::uint32_t n_threads, const std::string &table_filename, - const std::string &outfile, const RegressionType &alt_model, - const RegressionType &null_model, - const std::uint32_t test_factor_idx) { +radmeth(const bool show_progress, const bool more_na_info, + const std::uint32_t n_threads, const std::string &table_filename, + const std::string &outfile, const Regression &alt_model, + const Regression &null_model, + const std::uint32_t test_factor_idx) { static constexpr auto buf_size = 1024; static constexpr auto max_lines = 16384; @@ -120,8 +118,8 @@ that the design matrix and the proportion table are correctly formatted. std::vector n_bytes(max_lines, 0); std::vector lines(max_lines); - std::vector alt_models(n_threads, alt_model); - std::vector null_models(n_threads, null_model); + std::vector> alt_models(n_threads, alt_model); + std::vector> null_models(n_threads, null_model); const auto n_groups = alt_model.n_groups(); @@ -164,13 +162,14 @@ that the design matrix and the proportion table are correctly formatted. if (has_extreme_counts(t_alt_model)) return std::tuple{1.0, row_status::na_extreme_cnt}; - fit_regression_model_nano(t_alt_model, p_estim_alt, phi_estim_alt); + fit_regression_model_gamma_fdf(t_alt_model, p_estim_alt, + phi_estim_alt); t_null_model.mc = t_alt_model.mc; t_null_model.rowname = t_alt_model.rowname; - fit_regression_model_nano(t_null_model, p_estim_null, - phi_estim_null); + fit_regression_model_gamma_fdf(t_null_model, p_estim_null, + phi_estim_null); const double p_value = llr_test(t_null_model.max_loglik, t_alt_model.max_loglik); @@ -287,10 +286,10 @@ main_radmeth_nano(int argc, char *argv[]) { opt_parse.add_opt("factor", 'f', "a factor to test", true, test_factor); opt_parse.add_opt("tolerance", '\0', "numerical tolerance to test convergence", false, - RegressionNano::tolerance); + radmeth_optimize_params::tolerance); opt_parse.add_opt("max-iter", '\0', "max iterations when estimating parameters", false, - RegressionNano::max_iter); + radmeth_optimize_params::max_iter); std::vector leftover_args; opt_parse.parse(argc, argv, leftover_args); if (argc == 1 || opt_parse.help_requested()) { @@ -353,14 +352,14 @@ main_radmeth_nano(int argc, char *argv[]) { << '\n'; // clang-format on - RegressionNano alt_model; + Regression alt_model; alt_model.design = design; - RegressionNano null_model; + Regression null_model; null_model.design = null_design; const auto start_time = std::chrono::steady_clock::now(); - radmeth_nano(show_progress, more_na_info, n_threads, table_filename, - outfile, alt_model, null_model, test_factor_idx); + radmeth(show_progress, more_na_info, n_threads, table_filename, outfile, + alt_model, null_model, test_factor_idx); const auto stop_time = std::chrono::steady_clock::now(); if (verbose) From b893e22e6f8917ed379ec8087473ee2949f6e89e Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 19 Oct 2025 17:58:07 -0700 Subject: [PATCH 27/33] src/radmeth/radmeth_model_nano.hcpp: removing due to collapse into templated class --- src/radmeth/radmeth_model_nano.cpp | 80 -------------------------- src/radmeth/radmeth_model_nano.hpp | 90 ------------------------------ 2 files changed, 170 deletions(-) delete mode 100644 src/radmeth/radmeth_model_nano.cpp delete mode 100644 src/radmeth/radmeth_model_nano.hpp diff --git a/src/radmeth/radmeth_model_nano.cpp b/src/radmeth/radmeth_model_nano.cpp deleted file mode 100644 index 061afe80..00000000 --- a/src/radmeth/radmeth_model_nano.cpp +++ /dev/null @@ -1,80 +0,0 @@ -/* Copyright (C) 2025 Andrew D Smith - * - * Authors: Andrew D. Smith - * - * This program is free software: you can redistribute it and/or modify it - * under the terms of the GNU General Public License as published by the Free - * Software Foundation, either version 3 of the License, or (at your option) - * any later version. - * - * This program is distributed in the hope that it will be useful, but WITHOUT - * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for - * more details. - */ - -#include "radmeth_model_nano.hpp" - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -double RegressionNano::tolerance = 1e-4; -double RegressionNano::stepsize = 0.01; -std::uint32_t RegressionNano::max_iter = 250; - -void -RegressionNano::parse(const std::string &line) { - const auto first_ws = line.find_first_of(" \t"); - - // Parse the row name (must be like: "chr:position:strand:context") - bool failed = (first_ws == std::string::npos); - - auto field_s = line.data(); - auto field_e = line.data() + first_ws; - rowname = std::string{field_s, field_e}; - std::replace(std::begin(rowname), std::end(rowname), ':', '\t'); - if (failed) - throw std::runtime_error("failed to parse label from:\n" + line); - - // Parse the counts of total reads and methylated reads - const auto is_sep = [](const auto x) { return std::isspace(x); }; - const auto not_sep = [](const auto x) { return std::isdigit(x); }; - - const auto line_end = line.data() + std::size(line); - mcountsf mc1{}; - mc.clear(); - while (field_e != line_end) { - // get the total count - field_s = std::find_if(field_e + 1, line_end, not_sep); - field_e = std::find_if(field_s + 1, line_end, is_sep); - { - const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_reads); - failed = failed || (ec != std::errc{}); - } - - // get the methylated count - field_s = std::find_if(field_e + 1, line_end, not_sep); - field_e = std::find_if(field_s + 1, line_end, is_sep); - { -#ifdef __APPLE__ - const int ret = std::sscanf(field_s, "%lf", &mc1.n_meth); - failed = failed || (ret < 1); -#else - const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_meth); - failed = failed || (ec != std::errc{}); -#endif - } - - mc.push_back(mc1); - } - - if (failed) - throw std::runtime_error("failed to parse counts from:\n" + line); -} diff --git a/src/radmeth/radmeth_model_nano.hpp b/src/radmeth/radmeth_model_nano.hpp deleted file mode 100644 index 9e3c71ea..00000000 --- a/src/radmeth/radmeth_model_nano.hpp +++ /dev/null @@ -1,90 +0,0 @@ -/* Copyright (C) 2025 Andrew D Smith - * - * Author: Andrew D Smith - * - * This program is free software: you can redistribute it and/or modify it - * under the terms of the GNU General Public License as published by the Free - * Software Foundation, either version 3 of the License, or (at your option) - * any later version. - * - * This program is distributed in the hope that it will be useful, but WITHOUT - * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for - * more details. - */ - -#ifndef RADMETH_MODEL_NANO_HPP -#define RADMETH_MODEL_NANO_HPP - -#include "radmeth_design.hpp" - -#include -#include -#include -#include -#include - -struct mcountsf { - std::uint32_t n_reads{}; - double n_meth{}; -}; - -[[nodiscard]] inline std::istream & -operator>>(std::istream &in, mcountsf &rm) { - return in >> rm.n_reads >> rm.n_meth; -} - -struct vars_cache { - double p{}; - double a{}; - double b{}; - double lgamma_a{}; - double lgamma_b{}; - double lgamma_a_b{}; - double digamma_a{}; - double digamma_b{}; - double digamma_a_b{}; -}; - -struct RegressionNano { - static double tolerance; // 1e-3; - static double stepsize; // 0.001; - static std::uint32_t max_iter; // 250; - - Design design; - std::string rowname; - std::vector mc; - double max_loglik{}; - - std::vector cache; // scratch space - - void - parse(const std::string &line); - - [[nodiscard]] std::size_t - n_factors() const { - return design.n_factors(); - } - - [[nodiscard]] std::size_t - n_groups() const { - return design.n_groups(); - } - - [[nodiscard]] std::size_t - n_params() const { - return n_factors() + 1; - } - - [[nodiscard]] std::size_t - props_size() const { - return std::size(mc); - } - - [[nodiscard]] std::size_t - n_samples() const { - return design.n_samples(); - } -}; - -#endif // RADMETH_MODEL_NANO_HPP From cba7966018c85f8f29371acc8f4a2fa4ad709b1a Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 19 Oct 2025 17:58:56 -0700 Subject: [PATCH 28/33] Makefile.am: updates for (ongoing) change to names of radmeth sources --- Makefile.am | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/Makefile.am b/Makefile.am index e24502cd..58a1c915 100644 --- a/Makefile.am +++ b/Makefile.am @@ -189,12 +189,12 @@ libdnmtools_a_SOURCES += \ # ADS: additional radmeth sources to help isolate the parts using GSL for # optimization libdnmtools_a_SOURCES += \ - src/radmeth/radmeth_optimize.hpp \ - src/radmeth/radmeth_optimize_nano.hpp \ + src/radmeth/radmeth_optimize_series.hpp \ + src/radmeth/radmeth_optimize_gamma.hpp \ + src/radmeth/radmeth_optimize_params.hpp \ src/radmeth/radmeth_model.hpp \ src/radmeth/radmeth_utils.hpp \ - src/radmeth/radmeth_design.hpp \ - src/radmeth/radmeth_model_nano.hpp + src/radmeth/radmeth_design.hpp LDADD = libdnmtools.a src/abismal/libabismal.a src/smithlab_cpp/libsmithlab_cpp.a @@ -241,11 +241,9 @@ dnmtools_SOURCES += src/amrfinder/amrtester.cpp dnmtools_SOURCES += src/radmeth/dmr.cpp dnmtools_SOURCES += src/radmeth/methdiff.cpp dnmtools_SOURCES += src/radmeth/radmeth_utils.cpp -dnmtools_SOURCES += src/radmeth/radmeth_optimize.cpp -dnmtools_SOURCES += src/radmeth/radmeth_optimize_nano.cpp +dnmtools_SOURCES += src/radmeth/radmeth_optimize_series.cpp +dnmtools_SOURCES += src/radmeth/radmeth_optimize_gamma.cpp dnmtools_SOURCES += src/radmeth/radmeth_design.cpp -dnmtools_SOURCES += src/radmeth/radmeth_model.cpp -dnmtools_SOURCES += src/radmeth/radmeth_model_nano.cpp dnmtools_SOURCES += src/radmeth/radmeth.cpp dnmtools_SOURCES += src/radmeth/radmeth_nano.cpp dnmtools_SOURCES += src/radmeth/radmeth-adjust.cpp From 2541eedaca0ba0db227186106b527718f29ba654 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 19 Oct 2025 19:28:48 -0700 Subject: [PATCH 29/33] src/radmeth/radmeth_optimize_series.cpp: cleaning up junk code left in from experiments --- src/radmeth/radmeth_optimize_series.cpp | 164 +----------------------- 1 file changed, 2 insertions(+), 162 deletions(-) diff --git a/src/radmeth/radmeth_optimize_series.cpp b/src/radmeth/radmeth_optimize_series.cpp index e23bd367..5f92ac11 100644 --- a/src/radmeth/radmeth_optimize_series.cpp +++ b/src/radmeth/radmeth_optimize_series.cpp @@ -27,9 +27,6 @@ #include #include -#include ///////////////////// -#include ///////////////////// - [[nodiscard]] static inline auto logistic(const double x) -> double { return 1.0 / (1.0 / std::exp(x) + 1.0); @@ -394,71 +391,6 @@ get_cumulative(const std::vector &group_id, [](cumul_counts &c) -> std::vector & { return c.d_counts; }); } -static inline auto -conjugate_pr(std::vector ¶ms, Regression ®, - const double alpha) -> int { - const auto n = std::size(params); - const auto max_iter = radmeth_optimize_params::max_iter; - const auto tol = radmeth_optimize_params::tolerance; - - std::vector g(n); - std::vector g_prev(n); - gradient(params, reg, g); - - auto d = g; - std::transform(std::begin(d), std::end(d), std::begin(d), - [](const auto x) { return -1.0 * x; }); - - // const auto ll = log_likelihood(params, reg); - // std::cout << "ll=" << ll << std::endl; - - double gdot{}; - for (const double gi : g) - gdot += gi * gi; - double denom = gdot; - double gnorm = std::sqrt(gdot); - - std::size_t iter{}; - for (iter = 0; iter < max_iter; ++iter) { - // std::memcpy(g_prev.data(), g.data(), sizeof(double) * n); - // update parameters based on directions and step size - for (std::size_t i = 0; i < n; ++i) - params[i] += alpha * d[i]; - - gradient(params, reg, g); - - // get numerator and for norm of gradient - gdot = 0.0; - for (const double gi : g) - gdot += gi * gi; - gnorm = std::sqrt(gdot); - - // stop if norm of gradient is small enough - if (gnorm < tol) - break; - - // // get denom - // double denom = 0.0; - // for (const double gpi : g_prev) - // denom += gpi * gpi; - - // update conjugate directions - const double beta = std::max(0.0, gdot / denom); - for (auto i = 0u; i < n; ++i) - d[i] = -g[i] + beta * d[i]; - - denom = gdot; - } - return iter; - // if (gnorm < tol) { - // std::cout << "Converged in " << iter << " iterations" << "\t" << gnorm - // << '\t' << tol << std::endl; - // } - // else { - // std::cout << "Did not converge in " << max_iter << " iterations\n"; - // } -} - void fit_regression_model(Regression &r, std::vector &p_estimates, @@ -486,61 +418,6 @@ fit_regression_model(Regression &r, }; // clang-format on - // std::vector x(n_params, 0.0); - // // for (auto i = 0u; i < n_params; ++i) - // // x[i] = gsl_vector_get(param_estimates, i); - // x.back() = -2.5; - // bool all_finite = false; - // auto start_time_a = std::chrono::steady_clock::now(); - // auto ii = 0; - // for (ii = 0; ii < 4 && !all_finite; ++ii) { - // double alpha = std::pow(10.0, -1.0 * (ii + 2)); - // std::cout << "alpha=" << alpha << std::endl; - // std::vector x_out(x); - // conjugate_pr(x_out, r, alpha); - // all_finite = true; - // for (auto i = 0u; i < std::size(x_out) && all_finite; ++i) - // all_finite = std::isfinite(x_out[i]); - // if (all_finite) - // x = x_out; - // } - // // if (!all_finite) { - // // alpha = 1e-03; - // // x_out = x; - // // ii = conjugate_pr(x, r, alpha); - // // all_finite = true; - // // for (auto i = 0u; i < std::size(x) && all_finite; ++i) - // // all_finite = std::isfinite(x[i]); - // // if (!all_finite) { - // // alpha = 1e-04; - // // x_out = x; - // // ii = conjugate_pr(x, r, alpha); - // // all_finite = true; - // // for (auto i = 0u; i < std::size(x) && all_finite; ++i) - // // all_finite = std::isfinite(x[i]); - // // if (!all_finite) { - // // alpha = 1e-05; - // // x_out = x; - // // ii = conjugate_pr(x, r, alpha); - // // all_finite = true; - // // for (auto i = 0u; i < std::size(x) && all_finite; ++i) - // // all_finite = std::isfinite(x[i]); - // // if (!all_finite) { - // // alpha = 1e-06; - // // x_out = x; - // // ii = conjugate_pr(x, r, alpha); - // // all_finite = true; - // // for (auto i = 0u; i < std::size(x) && all_finite; ++i) - // // all_finite = std::isfinite(x[i]); - // // } - // // } - // // } - // // } - // auto stop_time_a = std::chrono::steady_clock::now(); - // // std::cout << "[total time: " << (stop_time_a - start_time_a).count() << - // // "]\t" - // // << ii << std::endl; - // parameters: 0 for "p" params and final one is a constant for dispersion auto params = gsl_vector_alloc(n_params); gsl_vector_set_all(params, 0.0); @@ -551,59 +428,23 @@ fit_regression_model(Regression &r, // - gsl_multimin_fdfminimizer_conjugate_fr // - gsl_multimin_fdfminimizer_vector_bfgs2 // - gsl_multimin_fdfminimizer_steepest_descent - - const auto minimizer = gsl_multimin_fdfminimizer_conjugate_fr; + const auto minimizer = gsl_multimin_fdfminimizer_vector_bfgs2; auto s = gsl_multimin_fdfminimizer_alloc(minimizer, n_params); gsl_multimin_fdfminimizer_set(s, &loglik_bundle, params, stepsize, tol); - // gsl_multimin_function minex_func; - - // minex_func.n = n_params; - // minex_func.f = &neg_loglik; - // minex_func.params = static_cast(&r); - - // /* Set initial step sizes to 1 */ - // auto ss = gsl_vector_alloc(n_params); - // gsl_vector_set_all(ss, 1.0); - - // const auto minimizer = gsl_multimin_fminimizer_nmsimplex2; - // auto s = gsl_multimin_fminimizer_alloc(minimizer, n_params); - // gsl_multimin_fminimizer_set(s, &minex_func, params, ss); - - auto start_time_b = std::chrono::steady_clock::now(); int status = 0; std::size_t iter = 0; - double size{}; do { - status = gsl_multimin_fdfminimizer_iterate(s); // one iter and get - // status - // status = gsl_multimin_fminimizer_iterate(s); // one iter and get status + status = gsl_multimin_fdfminimizer_iterate(s); // one iter and get status if (status) break; // check status from gradient - status = gsl_multimin_test_gradient(s->gradient, tol); - - // size = gsl_multimin_fminimizer_size(s); - // status = gsl_multimin_test_size(size, 1e-2); - } while (status == GSL_CONTINUE && ++iter < max_iter); // ADS: can't use (status != GSL_SUCCESS) - auto stop_time_b = std::chrono::steady_clock::now(); - // std::cout << "[total time: " - // << (stop_time_b - start_time_b).count() - // // - - // // (stop_time_a - start_time_a).count() - // << "]\t" << iter << std::endl; const auto param_estimates = gsl_multimin_fdfminimizer_x(s); - // const auto param_estimates = gsl_multimin_fminimizer_x(s); - // for (auto i = 0u; i < n_params; ++i) - // std::cout << "x[i]=" << gsl_vector_get(param_estimates, i) << std::endl; - - // gsl_vector_fprintf(stdout, param_estimates, "%g"); - const auto &groups = r.design.groups; p_estimates.clear(); for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) @@ -614,6 +455,5 @@ fit_regression_model(Regression &r, r.max_loglik = log_likelihood(param_estimates, r); gsl_multimin_fdfminimizer_free(s); - // gsl_multimin_fminimizer_free(s); gsl_vector_free(params); } From 7a755741a7feb0e2c9ed89d6e375d0f721d4820b Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 19 Oct 2025 19:49:02 -0700 Subject: [PATCH 30/33] src/radmeth/radmeth_optimize_series.cpp: cleaning more junk code --- src/radmeth/radmeth_optimize_series.cpp | 151 ------------------------ 1 file changed, 151 deletions(-) diff --git a/src/radmeth/radmeth_optimize_series.cpp b/src/radmeth/radmeth_optimize_series.cpp index 5f92ac11..bd54f087 100644 --- a/src/radmeth/radmeth_optimize_series.cpp +++ b/src/radmeth/radmeth_optimize_series.cpp @@ -39,13 +39,6 @@ get_p(const std::vector &v, const gsl_vector *params) -> double { return logistic(std::inner_product(a, a + std::size(v), params->data, 0.0)); } -template -[[nodiscard]] static inline auto -get_p(const std::vector &v, const std::vector ¶ms) -> double { - const auto a = v.data(); - return logistic(std::inner_product(a, a + std::size(v), params.data(), 0.0)); -} - static inline auto set_max_r_count(Regression ®) { const auto &cumul = reg.cumul; @@ -117,45 +110,6 @@ log_likelihood(const gsl_vector *params, Regression ®) { return log_lik; } -[[nodiscard]] static double -log_likelihood(const std::vector ¶ms, - Regression ®) { - const auto phi = logistic(params[reg.design.n_factors()]); - const auto one_minus_phi = 1.0 - phi; - - const auto n_groups = reg.n_groups(); - const auto &groups = reg.design.groups; - const auto &cumul = reg.cumul; - - // ADS: precompute log1p(phi * (k - 1.0)), reused for each group - get_cached_log1p_factors(reg, phi); - const auto &log1p_fact_v = reg.cache_log1p_factors; - - double log_lik = 0.0; - for (std::size_t g_idx = 0; g_idx < n_groups; ++g_idx) { - const auto p = get_p(groups[g_idx], params); - const auto one_minus_p = 1.0 - p; - - const auto term1 = one_minus_phi * p; - const auto &cumul_y = cumul[g_idx].m_counts; - std::size_t sz = std::size(cumul_y); - for (std::size_t k = 0; k < sz; ++k) - log_lik += cumul_y[k] * std::log(term1 + phi * k); - - const auto term2 = one_minus_phi * one_minus_p; - const auto &cumul_d = cumul[g_idx].d_counts; - sz = std::size(cumul_d); - for (std::size_t k = 0; k < sz; ++k) - log_lik += cumul_d[k] * std::log(term2 + phi * k); - - const auto &cumul_n = cumul[g_idx].r_counts; - sz = std::size(cumul_n); - for (std::size_t k = 0; k < sz; ++k) - log_lik -= cumul_n[k] * log1p_fact_v[k]; - } - return -log_lik; -} - static void gradient(const gsl_vector *params, Regression ®, gsl_vector *output) { @@ -220,111 +174,6 @@ gradient(const gsl_vector *params, Regression ®, gsl_vector_set(output, n_factors, disp_deriv * (phi * one_minus_phi)); } -static void -gradient(const std::vector ¶ms, Regression ®, - std::vector &output) { - const auto n_factors = reg.design.n_factors(); - const auto phi = logistic(params[n_factors]); - const auto one_minus_phi = 1.0 - phi; - - const auto n_groups = reg.n_groups(); - const auto &groups = reg.design.groups; - const auto &cumul = reg.cumul; - - auto &p_v = reg.p_v; - for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) - p_v[g_idx] = get_p(groups[g_idx], params); - - get_cached_dispersion_effect(reg, phi); // (k-1)/(1 + phi(k-1)) - const auto &dispersion_effect = reg.cache_dispersion_effect; - - // init output to zero for all factors - output.clear(); - output.resize(std::size(params), 0.0); - auto &data = output; - - double disp_deriv = 0.0; - for (std::size_t g_idx = 0; g_idx < n_groups; ++g_idx) { - const auto p = p_v[g_idx]; - const auto one_minus_p = 1.0 - p; - - // if (p == 0.0 || one_minus_p == 0.0 || phi == 0.0 || one_minus_phi == 0.0) - // { - // std::cout << p << '\t' << one_minus_p << std::endl; - // std::cout << phi << '\t' << one_minus_phi << std::endl; - // } - - double deriv = 0.0; - - const auto denom_term1 = one_minus_phi * p; - const auto &cumul_y = cumul[g_idx].m_counts; - const auto y_lim = std::size(cumul_y); - for (auto k = 0u; k < y_lim; ++k) { - const auto common_factor = cumul_y[k] / (denom_term1 + phi * k); - deriv += common_factor; - disp_deriv += (k - p) * common_factor; - // if (!std::isfinite(disp_deriv)) { - // // clang-format off - // std::cout << "common_factor: " << common_factor << '\n' - // << "cumul_y[k]: " << cumul_y[k] << '\n' - // << "denom_term1: " << denom_term1 << '\n' - // << "phi: " << phi << '\n' - // << "k: " << k << '\n' - // << "p: " << p << '\n' - // << "one_minus_phi: " << one_minus_phi << '\n' - // << "phi: " << phi << '\n' - // << "params[n_factors]: " << params[n_factors] << '\n' - // << "(denom_term1 + phi * k): " << (denom_term1 + phi * k) - // << '\n' - // << std::endl; - // throw std::runtime_error("bad disp_deriv 1"); - // // clang-format on - // } - } - - const auto denom_term2 = one_minus_phi * one_minus_p; - const auto &cumul_d = cumul[g_idx].d_counts; - const auto d_lim = std::size(cumul_d); - for (auto k = 0u; k < d_lim; ++k) { - const auto common_factor = cumul_d[k] / (denom_term2 + phi * k); - deriv -= common_factor; - disp_deriv += (k - one_minus_p) * common_factor; - // if (!std::isfinite(disp_deriv)) - // throw std::runtime_error("bad disp_deriv 2"); - } - - const auto &cumul_n = cumul[g_idx].r_counts; - const auto n_lim = std::size(cumul_n); - for (auto k = 0u; k < n_lim; ++k) { - disp_deriv -= cumul_n[k] * dispersion_effect[k]; - // if (!std::isfinite(disp_deriv)) - // throw std::runtime_error("bad disp_deriv 3"); - } - - const auto &g = groups[g_idx]; - const auto denom_term1_one_minus_p = denom_term1 * one_minus_p; - for (auto fact_idx = 0u; fact_idx < n_factors; ++fact_idx) { - const auto level = g[fact_idx]; - if (level == 0) - continue; - data[fact_idx] += deriv * (denom_term1_one_minus_p * level); - } - } - // if (!std::isfinite(disp_deriv)) - // throw std::runtime_error("bad disp_deriv"); - output[n_factors] = disp_deriv * (phi * one_minus_phi); - // for (auto i = 0u; i <= n_factors; ++i) { - // if (!std::isfinite(output[i])) - // throw std::runtime_error("bad2"); - // } - for (auto i = 0u; i <= n_factors; ++i) { - output[i] *= -1.0; - // if (!std::isfinite(output[i])) - // throw std::runtime_error("bad"); - } - // std::cout << "grad=" << disp_deriv * (phi * one_minus_phi) << std::endl; -} - [[nodiscard]] static double neg_loglik(const gsl_vector *params, void *object) { auto reg = static_cast *>(object); From f4b61216fd4e37de7850e4eb4be323e3df7d3921 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 19 Oct 2025 19:51:22 -0700 Subject: [PATCH 31/33] src/radmeth/radmeth_optimize_gamma.cpp: making bfgs2 the default and cleaning up --- src/radmeth/radmeth_optimize_gamma.cpp | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/radmeth/radmeth_optimize_gamma.cpp b/src/radmeth/radmeth_optimize_gamma.cpp index 592f5fb0..d4efc75b 100644 --- a/src/radmeth/radmeth_optimize_gamma.cpp +++ b/src/radmeth/radmeth_optimize_gamma.cpp @@ -13,8 +13,8 @@ * more details. */ -#include "radmeth_model.hpp" #include "radmeth_optimize_gamma.hpp" +#include "radmeth_model.hpp" #include "radmeth_optimize_params.hpp" #include @@ -268,15 +268,14 @@ template static void fit_regression_model(Regression &r, std::vector &p_estimates, double &dispersion_estimate) { - static constexpr auto init_dispersion_param = -1.0; + static constexpr auto init_dispersion_param = -2.5; const auto stepsize = radmeth_optimize_params::stepsize; const auto max_iter = radmeth_optimize_params::max_iter; const auto n_groups = r.n_groups(); r.cache.resize(n_groups); // make sure scratch space is allocated const std::size_t n_params = r.n_params(); - const auto tol = - std::sqrt(n_params) * r.n_samples() * radmeth_optimize_params::tolerance; + const auto tol = radmeth_optimize_params::tolerance; // set the parameters: zero for "p" parameters and the final one for // dispersion using the constant @@ -300,14 +299,13 @@ fit_regression_model(Regression &r, std::vector &p_estimates, int status = 0; std::size_t iter = 0; - double size{}; do { status = gsl_multimin_fminimizer_iterate(s); // one iter and get status if (status) break; - size = gsl_multimin_fminimizer_size(s); + const auto size = gsl_multimin_fminimizer_size(s); status = gsl_multimin_test_size(size, tol); } while (status == GSL_CONTINUE && ++iter < max_iter); // ADS: can't use (status != GSL_SUCCESS) @@ -332,7 +330,7 @@ template static void fit_regression_model_fdf(Regression &r, std::vector &p_estimates, double &dispersion_estimate) { - static constexpr auto init_dispersion_param = -1.0; + static constexpr auto init_dispersion_param = -2.5; const auto stepsize = radmeth_optimize_params::stepsize; const auto max_iter = radmeth_optimize_params::max_iter; const auto n_groups = r.n_groups(); @@ -362,7 +360,7 @@ fit_regression_model_fdf(Regression &r, std::vector &p_estimates, // - gsl_multimin_fdfminimizer_conjugate_fr // - gsl_multimin_fdfminimizer_vector_bfgs2 // - gsl_multimin_fdfminimizer_steepest_descent - const auto minimizer = gsl_multimin_fdfminimizer_conjugate_pr; + const auto minimizer = gsl_multimin_fdfminimizer_vector_bfgs2; auto s = gsl_multimin_fdfminimizer_alloc(minimizer, n_params); gsl_multimin_fdfminimizer_set(s, &loglik_bundle, params, stepsize, tol); From 0340295a38b8680410ba9343c1cd4f43c36d0c96 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Sun, 19 Oct 2025 19:52:14 -0700 Subject: [PATCH 32/33] src/radmeth/radmeth.cpp and radmeth_nano.cpp: making these two sources more similar and improving output messages --- src/radmeth/radmeth.cpp | 59 ++++++++++++++++++------------------ src/radmeth/radmeth_nano.cpp | 56 +++++++++++++++++----------------- 2 files changed, 57 insertions(+), 58 deletions(-) diff --git a/src/radmeth/radmeth.cpp b/src/radmeth/radmeth.cpp index 0da6bfc7..68eab100 100644 --- a/src/radmeth/radmeth.cpp +++ b/src/radmeth/radmeth.cpp @@ -76,12 +76,12 @@ enum class row_status : std::uint8_t { na_extreme_cnt, }; -template +template static void radmeth(const bool show_progress, const bool more_na_info, const std::uint32_t n_threads, const std::string &table_filename, - const std::string &outfile, const RegressionType &alt_model, - const RegressionType &null_model, const std::uint32_t test_factor_idx) { + const std::string &outfile, const Regression &alt_model, + const Regression &null_model, const std::uint32_t test_factor_idx) { static constexpr auto buf_size = 1024; static constexpr auto max_lines = 16384; @@ -121,8 +121,8 @@ that the design matrix and the proportion table are correctly formatted. std::vector n_bytes(max_lines, 0); std::vector lines(max_lines); - std::vector alt_models(n_threads, alt_model); - std::vector null_models(n_threads, null_model); + std::vector> alt_models(n_threads, alt_model); + std::vector> null_models(n_threads, null_model); const auto n_groups = alt_model.n_groups(); @@ -183,48 +183,54 @@ that the design matrix and the proportion table are correctly formatted. const auto status = std::get<1>(p_val_status); n_bytes[b] = [&] { + auto bufsize = std::size(bufs[b]); // clang-format off const int n_prefix_bytes = - std::sprintf(bufs[b].data(), "%s\t", t_alt_model.rowname.data()); + std::snprintf(bufs[b].data(), bufsize, + "%s\t", t_alt_model.rowname.data()); // clang-format on if (n_prefix_bytes < 0) return n_prefix_bytes; + bufsize -= n_prefix_bytes; auto cursor = bufs[b].data() + n_prefix_bytes; const int n_pval_bytes = [&] { if (status == row_status::ok) - return std::sprintf(cursor, "%.6g", p_val); + return std::snprintf(cursor, bufsize, "%.6g", p_val); if (!more_na_info || status == row_status::na) - return std::sprintf(cursor, "NA"); + return std::snprintf(cursor, bufsize, "NA"); if (status == row_status::na_extreme_cnt) - return std::sprintf(cursor, "NA_EXTREME_CNT"); + return std::snprintf(cursor, bufsize, "NA_EXTREME_CNT"); // if (status == row_status::na_low_cov) - return std::sprintf(cursor, "NA_LOW_COV"); + return std::snprintf(cursor, bufsize, "NA_LOW_COV"); }(); if (n_pval_bytes < 0) return n_pval_bytes; + bufsize -= n_pval_bytes; cursor += n_pval_bytes; const int n_param_bytes = [&] { - std::int32_t n_bytes = 0; + std::int32_t n_param_bytes = 0; if (p_estim_alt.empty()) p_estim_alt.resize(n_groups, 0.0); for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) { - const int n = std::sprintf(cursor, "\t%f", p_estim_alt[g_idx]); + const int n = + std::snprintf(cursor, bufsize, "\t%f", p_estim_alt[g_idx]); + bufsize -= n; cursor += n; if (n < 0) return n; - n_bytes += n; + n_param_bytes += n; } const auto od = overdispersion_factor(n_samples, phi_estim_alt); - const int n = std::sprintf(cursor, "\t%f\n", od); + const int n = std::snprintf(cursor, bufsize, "\t%f\n", od); if (n < 0) return n; - n_bytes += n; - return n_bytes; + n_param_bytes += n; + return n_param_bytes; }(); if (n_param_bytes < 0) @@ -314,9 +320,6 @@ main_radmeth(int argc, char *argv[]) { const std::string table_filename(leftover_args.back()); /****************** END COMMAND LINE OPTIONS *****************/ - if (verbose) - std::cerr << "design table filename: " << design_filename << "\n\n"; - Design design = Design::read_design(design_filename); // Check that provided test factor name exists and find its index. @@ -324,9 +327,6 @@ main_radmeth(int argc, char *argv[]) { ensure_sample_order(table_filename, design); - if (verbose) - std::cerr << "Alternate model:\n" << design << '\n'; - // verify that the design includes more than one level for the // test factor if (!design.has_two_values(test_factor_idx)) { @@ -336,22 +336,21 @@ main_radmeth(int argc, char *argv[]) { } const Design null_design = design.drop_factor(test_factor_idx); - if (verbose) - std::cerr << "Null model:\n" << null_design << '\n'; // clang-format off if (verbose) - std::cerr << "Output columns:\n" + std::cerr << "design table filename: " << design_filename << "\n\n" + << "Alternate model:\n" << design << '\n' + << "Null model:\n" << null_design << '\n' + << "Output columns:\n" << "(1) chrom\n" << "(2) position\n" << "(3) strand\n" << "(4) cytosine context\n" << "(5) p-value\n" - << "(6) estimated methylation 0\n" - << "(7) estimated methylation 1\n" - << "(8) overdispersion factor (variance above binomial)\n" - << "estimated methylation is for test factor value (0 or 1)\n" - << '\n'; + << "(6) overdispersion factor (variance above binomial)\n" + << "(7-) estimated methylation for each unique group\n" + << "groups are defined by combinations of factor levels\n\n"; // clang-format on Regression alt_model; diff --git a/src/radmeth/radmeth_nano.cpp b/src/radmeth/radmeth_nano.cpp index 4bb41a73..0932cb84 100644 --- a/src/radmeth/radmeth_nano.cpp +++ b/src/radmeth/radmeth_nano.cpp @@ -73,12 +73,12 @@ enum class row_status : std::uint8_t { na_extreme_cnt, }; +template static void radmeth(const bool show_progress, const bool more_na_info, const std::uint32_t n_threads, const std::string &table_filename, - const std::string &outfile, const Regression &alt_model, - const Regression &null_model, - const std::uint32_t test_factor_idx) { + const std::string &outfile, const Regression &alt_model, + const Regression &null_model, const std::uint32_t test_factor_idx) { static constexpr auto buf_size = 1024; static constexpr auto max_lines = 16384; @@ -118,8 +118,8 @@ that the design matrix and the proportion table are correctly formatted. std::vector n_bytes(max_lines, 0); std::vector lines(max_lines); - std::vector> alt_models(n_threads, alt_model); - std::vector> null_models(n_threads, null_model); + std::vector> alt_models(n_threads, alt_model); + std::vector> null_models(n_threads, null_model); const auto n_groups = alt_model.n_groups(); @@ -182,48 +182,54 @@ that the design matrix and the proportion table are correctly formatted. const auto status = std::get<1>(p_val_status); n_bytes[b] = [&] { + auto bufsize = std::size(bufs[b]); // clang-format off const int n_prefix_bytes = - std::sprintf(bufs[b].data(), "%s\t", t_alt_model.rowname.data()); + std::snprintf(bufs[b].data(), bufsize, + "%s\t", t_alt_model.rowname.data()); // clang-format on if (n_prefix_bytes < 0) return n_prefix_bytes; + bufsize -= n_prefix_bytes; auto cursor = bufs[b].data() + n_prefix_bytes; const int n_pval_bytes = [&] { if (status == row_status::ok) - return std::sprintf(cursor, "%.6g", p_val); + return std::snprintf(cursor, bufsize, "%.6g", p_val); if (!more_na_info || status == row_status::na) - return std::sprintf(cursor, "NA"); + return std::snprintf(cursor, bufsize, "NA"); if (status == row_status::na_extreme_cnt) - return std::sprintf(cursor, "NA_EXTREME_CNT"); + return std::snprintf(cursor, bufsize, "NA_EXTREME_CNT"); // if (status == row_status::na_low_cov) - return std::sprintf(cursor, "NA_LOW_COV"); + return std::snprintf(cursor, bufsize, "NA_LOW_COV"); }(); if (n_pval_bytes < 0) return n_pval_bytes; + bufsize -= n_pval_bytes; cursor += n_pval_bytes; const int n_param_bytes = [&] { - std::int32_t n_bytes = 0; + std::int32_t n_param_bytes = 0; if (p_estim_alt.empty()) p_estim_alt.resize(n_groups, 0.0); for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) { - const int n = std::sprintf(cursor, "\t%f", p_estim_alt[g_idx]); + const int n = + std::snprintf(cursor, bufsize, "\t%f", p_estim_alt[g_idx]); + bufsize -= n; cursor += n; if (n < 0) return n; - n_bytes += n; + n_param_bytes += n; } const auto od = overdispersion_factor(n_samples, phi_estim_alt); - const int n = std::sprintf(cursor, "\t%f\n", od); + const int n = std::snprintf(cursor, bufsize, "\t%f\n", od); if (n < 0) return n; - n_bytes += n; - return n_bytes; + n_param_bytes += n; + return n_param_bytes; }(); if (n_param_bytes < 0) @@ -313,9 +319,6 @@ main_radmeth_nano(int argc, char *argv[]) { const std::string table_filename(leftover_args.back()); /****************** END COMMAND LINE OPTIONS *****************/ - if (verbose) - std::cerr << "design table filename: " << design_filename << "\n\n"; - Design design = Design::read_design(design_filename); // Check that provided test factor name exists and find its index. @@ -323,9 +326,6 @@ main_radmeth_nano(int argc, char *argv[]) { ensure_sample_order(table_filename, design); - if (verbose) - std::cerr << "Alternate model:\n" << design << '\n'; - // verify that the design includes more than one level for the // test factor if (!design.has_two_values(test_factor_idx)) { @@ -338,18 +338,18 @@ main_radmeth_nano(int argc, char *argv[]) { // clang-format off if (verbose) - std::cerr << "Null model:\n" << null_design << '\n' + std::cerr << "design table filename: " << design_filename << "\n\n" + << "Alternate model:\n" << design << '\n' + << "Null model:\n" << null_design << '\n' << "Output columns:\n" << "(1) chrom\n" << "(2) position\n" << "(3) strand\n" << "(4) cytosine context\n" << "(5) p-value\n" - << "(6) estimated methylation 0\n" - << "(7) estimated methylation 1\n" - << "(8) overdispersion factor (variance above binomial)\n" - << "estimated methylation is for test factor value (0 or 1)\n" - << '\n'; + << "(6) overdispersion factor (variance above binomial)\n" + << "(7-) estimated methylation for each unique group\n" + << "groups are defined by combinations of factor levels\n\n"; // clang-format on Regression alt_model; From 205676e9120095566d3ef39fffbef6f673f5d9b9 Mon Sep 17 00:00:00 2001 From: Andrew D Smith Date: Wed, 22 Oct 2025 14:31:07 -0700 Subject: [PATCH 33/33] macOS clang fixes: from_chars for float to sscanf and sprintf to snprintf --- src/analysis/nanopore.cpp | 38 ++++++++++++++++++++--------------- src/radmeth/methdiff.cpp | 20 +++++++++--------- src/radmeth/radmeth_model.hpp | 8 +++++++- 3 files changed, 39 insertions(+), 27 deletions(-) diff --git a/src/analysis/nanopore.cpp b/src/analysis/nanopore.cpp index 6eb3ba53..cdf29b4c 100644 --- a/src/analysis/nanopore.cpp +++ b/src/analysis/nanopore.cpp @@ -845,7 +845,8 @@ struct read_processor { if (chrom_name == nullptr) throw std::runtime_error("failed to identify chrom for tid: " + std::to_string(tid)); - const auto chrom_name_offset = std::sprintf(buffer, "%s\t", chrom_name); + const auto chrom_name_offset = + std::snprintf(buffer, buf_size, "%s\t", chrom_name); if (chrom_name_offset < 0) throw std::runtime_error("failed to write to output buffer"); auto buffer_after_chrom = buffer + chrom_name_offset; @@ -869,14 +870,16 @@ struct read_processor { const double mods = counts[chrom_posn].get_mods(is_c) / denom; // clang-format off - const int r = std::sprintf(buffer_after_chrom, out_fmt, - chrom_posn, - is_c ? '+' : '-', - tag_values[the_tag], - mods, - n_reads, - hydroxy, - methyl); + const int r = std::snprintf(buffer_after_chrom, + buf_size - chrom_name_offset, + out_fmt, + chrom_posn, + is_c ? '+' : '-', + tag_values[the_tag], + mods, + n_reads, + hydroxy, + methyl); // clang-format on if (r < 0) @@ -900,7 +903,8 @@ struct read_processor { if (chrom_name == nullptr) throw std::runtime_error("failed to identify chrom for tid: " + std::to_string(tid)); - const auto chrom_name_offset = std::sprintf(buffer, "%s\t", chrom_name); + const auto chrom_name_offset = + std::snprintf(buffer, buf_size, "%s\t", chrom_name); if (chrom_name_offset < 0) throw std::runtime_error("failed to write to output buffer"); auto buffer_after_chrom = buffer + chrom_name_offset; @@ -935,12 +939,14 @@ struct read_processor { const auto mods = hydroxy + methyl; // clang-format off - const int r = std::sprintf(buffer_after_chrom, out_fmt, - chrom_posn - 1, // for previous position - mods, - n_reads, - hydroxy, - methyl); + const int r = std::snprintf(buffer_after_chrom, + buf_size - chrom_name_offset, + out_fmt, + chrom_posn - 1, // for previous position + mods, + n_reads, + hydroxy, + methyl); // clang-format on if (r < 0) diff --git a/src/radmeth/methdiff.cpp b/src/radmeth/methdiff.cpp index 2fe12816..33d4d474 100644 --- a/src/radmeth/methdiff.cpp +++ b/src/radmeth/methdiff.cpp @@ -105,16 +105,16 @@ write_methdiff_site(T &out, const MSite &a, const MSite &b, static char buffer[buf_size]; // clang-format off - const int r = std::sprintf(buffer, out_fmt, - a.chrom.data(), - a.pos, - a.strand, - a.context.data(), - diffscore, - a.n_meth(), - a.n_unmeth(), - b.n_meth(), - b.n_unmeth()); + const int r = std::snprintf(buffer, buf_size, out_fmt, + a.chrom.data(), + a.pos, + a.strand, + a.context.data(), + diffscore, + a.n_meth(), + a.n_unmeth(), + b.n_meth(), + b.n_unmeth()); // clang-format on if (r < 0) throw std::runtime_error("failed to write to output buffer"); diff --git a/src/radmeth/radmeth_model.hpp b/src/radmeth/radmeth_model.hpp index a202664e..9be0271a 100644 --- a/src/radmeth/radmeth_model.hpp +++ b/src/radmeth/radmeth_model.hpp @@ -135,9 +135,15 @@ Regression::parse(const std::string &line) { field_s = std::find_if(field_e + 1, line_end, not_sep); field_e = std::find_if(field_s + 1, line_end, is_sep); + { +#ifdef __APPLE__ + const int ret = std::sscanf(field_s, "%lf", &mc1.n_meth); + failed = failed || (ret < 1); +#else const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_meth); - failed = failed || (ptr != field_e); + failed = failed || (ec != std::errc{}); +#endif } mc.push_back(mc1);