From c9a69a678e259ba8f5d5b2f5b6bcbab39627aeaa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tommi=20M=C3=A4klin?= Date: Fri, 27 Oct 2023 18:49:50 +0300 Subject: [PATCH 01/11] Remove deprecated MPI implementation. --- include/CppRateRes_mpi.hpp | 220 ------------------------------------- src/cpprate_mpi.cpp | 172 ----------------------------- 2 files changed, 392 deletions(-) delete mode 100644 include/CppRateRes_mpi.hpp delete mode 100644 src/cpprate_mpi.cpp diff --git a/include/CppRateRes_mpi.hpp b/include/CppRateRes_mpi.hpp deleted file mode 100644 index 1f0b6c9..0000000 --- a/include/CppRateRes_mpi.hpp +++ /dev/null @@ -1,220 +0,0 @@ -// cpprate: Variable Selection in Black Box Methods with RelATive cEntrality (RATE) Measures -// https://github.com/tmaklin/cpprate -// Copyright (c) 2023 Tommi Mäklin (tommi@maklin.fi) -// -// BSD-3-Clause license -// -// Redistribution and use in source and binary forms, with or without -// modification, are permitted provided that the following conditions are -// met: -// -// (1) Redistributions of source code must retain the above copyright -// notice, this list of conditions and the following disclaimer. -// -// (2) Redistributions in binary form must reproduce the above copyright -// notice, this list of conditions and the following disclaimer in -// the documentation and/or other materials provided with the -// distribution. -// -// (3)The name of the author may not be used to -// endorse or promote products derived from this software without -// specific prior written permission. -// -// THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR -// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED -// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -// DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, -// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING -// IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -// POSSIBILITY OF SUCH DAMAGE. -// -#ifndef CPPRATE_CPPRATERES_MPI_HPP -#define CPPRATE_CPPRATERES_MPI_HPP - -#include "CppRateRes.hpp" - -#include -#include - -#include "cpprate_mpi_config.hpp" - -inline RATEd RATE_lowrank_mpi(Eigen::MatrixXd &f_draws, Eigen::SparseMatrix &design_matrix, const size_t n_snps, const size_t svd_rank, const double prop_var) { - // ## WARNING: Do not compile with -ffast-math - - // Setup MPI - int rank; - int n_tasks; - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - MPI_Comm_size(MPI_COMM_WORLD, &n_tasks); - - Eigen::VectorXd col_means_beta(0); - Eigen::MatrixXd v_Sigma_star(0, 0); - Eigen::MatrixXd Lambda_f(0, 0); - Eigen::MatrixXd svd_design_matrix_v(0, 0); - std::vector flat_Lambda(0); - if (rank == 0) { - Eigen::MatrixXd u; - decompose_design_matrix(design_matrix, svd_rank, prop_var, &u, &svd_design_matrix_v); - design_matrix.resize(0, 0); - col_means_beta = approximate_beta_means(f_draws, u, svd_design_matrix_v); - Eigen::MatrixXd Sigma_star = project_f_draws(f_draws, u); - u.resize(0, 0); - - v_Sigma_star = svd_design_matrix_v * Sigma_star.triangularView(); - Eigen::MatrixXd Lambda_chol = decompose_covariance_approximation(Sigma_star, svd_design_matrix_v, svd_rank); - f_draws.resize(0, 0); - Sigma_star.resize(0, 0); - - Eigen::MatrixXd Lambda = Eigen::MatrixXd::Zero(n_snps, n_snps); - Lambda.template selfadjointView().rankUpdate(Lambda_chol); - Lambda_f = Lambda.triangularView() * v_Sigma_star; - - flat_Lambda = flatten_triangular(Lambda); - -#pragma omp parallel for schedule(static) - for (size_t i = 0; i < v_Sigma_star.cols(); ++i) { - for (size_t j = 0; j < v_Sigma_star.rows(); ++j) { - v_Sigma_star(j, i) = std::log(std::abs(v_Sigma_star(j, i)) + 1e-16) + std::log(std::abs(Lambda_chol(j, i)) + 1e-16); - } - } - - svd_design_matrix_v.transposeInPlace(); - -#pragma omp parallel for schedule(static) - for (size_t i = 0; i < Lambda_f.cols(); ++i) { - for (size_t j = 0; j < Lambda_f.rows(); ++j) { - Lambda_f(j, i) = std::log(std::abs(Lambda_f(j, i)) + 1e-16); - } - } - } - f_draws.resize(0, 0); - design_matrix.resize(0, 0); - - // Already known (argument) - size_t flat_Lambda_size = n_snps * (n_snps + 1)/2; - - size_t Sigma_star_rows = v_Sigma_star.rows(); - size_t Sigma_star_cols = v_Sigma_star.cols(); - size_t Lambda_f_rows = Lambda_f.rows(); - size_t Lambda_f_cols = Lambda_f.cols(); - size_t svd_design_matrix_v_rows = svd_design_matrix_v.rows(); - size_t svd_design_matrix_v_cols = svd_design_matrix_v.cols(); - - // Broadcast sizes - MPI_Bcast(&Sigma_star_rows, 1, MPI_UNSIGNED_LONG_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast(&Sigma_star_cols, 1, MPI_UNSIGNED_LONG_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast(&Lambda_f_rows, 1, MPI_UNSIGNED_LONG_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast(&Lambda_f_cols, 1, MPI_UNSIGNED_LONG_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast(&svd_design_matrix_v_rows, 1, MPI_UNSIGNED_LONG_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast(&svd_design_matrix_v_cols, 1, MPI_UNSIGNED_LONG_LONG, 0, MPI_COMM_WORLD); - - // Initialize ranges for MPI - size_t n_snps_per_task = std::floor(n_snps/(double)n_tasks); - size_t start_id = rank * n_snps_per_task; - size_t end_id = std::min(n_snps, (rank + 1) * n_snps_per_task); - if (rank == (n_tasks - 1)) { - n_snps_per_task += n_snps - n_snps_per_task*n_tasks; - end_id = std::min(n_snps, (rank + 1) * n_snps_per_task); - } - - if (rank != 0) { - v_Sigma_star.resize(Sigma_star_rows, Sigma_star_cols); - Lambda_f.resize(Lambda_f_rows, Lambda_f_cols); - svd_design_matrix_v.resize(svd_design_matrix_v_rows, svd_design_matrix_v_cols); - col_means_beta.resize(n_snps_per_task); - flat_Lambda.resize(flat_Lambda_size); - } - - { - // Initializes the displacements and bufcounts. - int displacements[1024]; - int bufcounts[1024] = { 0 }; - - uint32_t sent_so_far = 0; - uint32_t n_obs_per_task = std::floor(n_snps/n_tasks); - for (uint16_t i = 0; i < n_tasks - 1; ++i) { - displacements[i] = sent_so_far; - bufcounts[i] = n_obs_per_task; - sent_so_far += bufcounts[i]; - } - displacements[n_tasks - 1] = sent_so_far; - bufcounts[n_tasks - 1] = n_snps - sent_so_far; - bufcounts[0] = 0; - - MPI_Scatterv(col_means_beta.data(), bufcounts, displacements, MPI_DOUBLE, col_means_beta.data(), n_snps_per_task, MPI_DOUBLE, 0, MPI_COMM_WORLD); - } - - { - // Initializes the displacements and bufcounts. - int displacements[1024]; - int bufcounts[1024] = { 0 }; - - uint32_t sent_so_far = 0; - for (uint16_t i = 0; i < n_tasks - 1; ++i) { - displacements[i] = sent_so_far; - bufcounts[i] = n_snps_per_task * svd_design_matrix_v_rows; - sent_so_far += bufcounts[i]; - } - displacements[n_tasks - 1] = sent_so_far; - bufcounts[n_tasks - 1] = (n_snps * svd_design_matrix_v_rows) - sent_so_far; - bufcounts[0] = 0; - - MPI_Scatterv(svd_design_matrix_v.data(), bufcounts, displacements, MPI_DOUBLE, svd_design_matrix_v.data(), n_snps_per_task * svd_design_matrix_v_rows, MPI_DOUBLE, 0, MPI_COMM_WORLD); - } - - // Broadcast variables needed by all processes - MPI_Bcast(v_Sigma_star.data(), v_Sigma_star.rows()*v_Sigma_star.cols(), MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast(Lambda_f.data(), Lambda_f.rows()*Lambda_f.cols(), MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast(flat_Lambda.data(), flat_Lambda_size, MPI_DOUBLE, 0, MPI_COMM_WORLD); - - std::vector log_KLD_partial(n_snps_per_task); - for (size_t i = 0; i < n_snps_per_task; ++i) { - log_KLD_partial[i] = dropped_predictor_kld_lowrank(flat_Lambda, Lambda_f, v_Sigma_star, svd_design_matrix_v.col(i), col_means_beta[i], start_id + i); - } - - std::vector KLD_partial; - std::transform(log_KLD_partial.begin(), log_KLD_partial.end(), std::back_inserter(KLD_partial), static_cast(std::exp)); - double KLD_sum_local = std::accumulate(KLD_partial.begin(), KLD_partial.end(), 0.0); - double KLD_sum_global = 0.0; - MPI_Allreduce(&KLD_sum_local, &KLD_sum_global, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - - std::vector RATE_partial = rate_from_kld(log_KLD_partial, KLD_sum_global); - - std::vector KLD(0); - std::vector RATE(0); - if (rank == 0) { - KLD.resize(n_snps); - RATE.resize(n_snps); - } - - { - // Initializes the displacements and bufcounts. - int displacements[1024]; - int bufcounts[1024] = { 0 }; - - uint32_t sent_so_far = 0; - uint32_t n_obs_per_task = std::floor(n_snps/n_tasks); - for (uint16_t i = 0; i < n_tasks - 1; ++i) { - displacements[i] = sent_so_far; - bufcounts[i] = n_obs_per_task; - sent_so_far += bufcounts[i]; - } - displacements[n_tasks - 1] = sent_so_far; - bufcounts[n_tasks - 1] = n_snps - sent_so_far; - - MPI_Gatherv(&KLD_partial.front(), n_snps_per_task, MPI_DOUBLE, &KLD.front(), bufcounts, displacements, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Gatherv(&RATE_partial.front(), n_snps_per_task, MPI_DOUBLE, &RATE.front(), bufcounts, displacements, MPI_DOUBLE, 0, MPI_COMM_WORLD); - } - - if (rank == 0) { - return RATEd(KLD, RATE); - } else { - return RATEd(); - } -} - -#endif diff --git a/src/cpprate_mpi.cpp b/src/cpprate_mpi.cpp deleted file mode 100644 index 4a781eb..0000000 --- a/src/cpprate_mpi.cpp +++ /dev/null @@ -1,172 +0,0 @@ -// cpprate: Variable Selection in Black Box Methods with RelATive cEntrality (RATE) Measures -// https://github.com/tmaklin/cpprate -// Copyright (c) 2023 Tommi Mäklin (tommi@maklin.fi) -// -// BSD-3-Clause license -// -// Redistribution and use in source and binary forms, with or without -// modification, are permitted provided that the following conditions are -// met: -// -// (1) Redistributions of source code must retain the above copyright -// notice, this list of conditions and the following disclaimer. -// -// (2) Redistributions in binary form must reproduce the above copyright -// notice, this list of conditions and the following disclaimer in -// the documentation and/or other materials provided with the -// distribution. -// -// (3)The name of the author may not be used to -// endorse or promote products derived from this software without -// specific prior written permission. -// -// THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR -// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED -// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -// DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, -// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING -// IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -// POSSIBILITY OF SUCH DAMAGE. -// -#include "cxxargs.hpp" - -#include -#include -#include - -#include "bxzstr.hpp" - -#include "CppRateRes.hpp" -#include "CppRateRes_mpi.hpp" - -#include "cpprate_mpi_config.hpp" - -bool CmdOptionPresent(char **begin, char **end, const std::string &option) { - return (std::find(begin, end, option) != end); -} - -void parse_args(int argc, char* argv[], cxxargs::Arguments &args) { - args.add_short_argument('f', "f-draws file (comma separated)"); - args.add_short_argument('x', "design matrix (comma separated)"); - args.add_short_argument('n', "Number of observations (rows in design matrix; columns in f-draws"); - args.add_short_argument('d', "Number of SNPs tested (cols in design matrix"); - args.add_short_argument('m', "Number of posterior samples (rows in f-draws)"); - args.add_short_argument('t', "Number of threads to use (default: 1)", 1); - args.add_long_argument("prop-var", "Proportion of variance to explain in lowrank factorization (default: 100%)", 1.1); - args.add_long_argument("low-rank", "Rank of the low-rank factorization (default: min(-n, -d))", 0); - args.add_long_argument("fullrank", "Run fullrank algorithm (default: false)", false); - args.add_long_argument("lowrank-dim", "Dimension of the lowrank approximation (default: min(n, d))", 0); - args.add_long_argument("help", "Print the help message.", false); - if (CmdOptionPresent(argv, argv+argc, "--help")) { - std::cout << "\n" + args.help() << '\n' << std::endl; - } - args.parse(argc, argv); -} - -int main(int argc, char* argv[]) { - cxxargs::Arguments args("cpprate", "Usage: cpprate -f -x -n -d -m "); - - int rc = MPI_Init(&argc, &argv); - int rank; - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - - size_t num_threads = 0; - if (rank == 0) { - try { - parse_args(argc, argv, args); - } catch (std::exception &e) { - std::cerr << "Parsing arguments failed:\n" - << std::string("\t") + std::string(e.what()) + "\n" - << "\trun cpprate with the --help option for usage instructions.\n"; - std::cerr << std::endl; - return 1; - } - num_threads = args.value('t'); - } - MPI_Bcast(&num_threads, 1, MPI_UNSIGNED_LONG_LONG, 0, MPI_COMM_WORLD); - omp_set_num_threads(num_threads); - - size_t n_snps; - size_t n_obs; - - Eigen::MatrixXd f_draws_mat; - Eigen::SparseMatrix design_matrix; - if (rank == 0) { - n_obs = args.value('n'); - n_snps = args.value('d'); - size_t n_f_draws = args.value('m'); - - std::vector f_draws; - bxz::ifstream in(args.value('f')); - std::string line; - while (std::getline(in, line)) { - std::stringstream parts(line); - std::string part; - while(std::getline(parts, part, ',')) { - f_draws.emplace_back(std::stold(part)); - } - } - in.close(); - - f_draws_mat = std::move(vec_to_dense_matrix(f_draws, n_f_draws, n_obs)); - - std::vector X; - bxz::ifstream in2(args.value('x')); - while (std::getline(in2, line)) { - std::stringstream parts(line); - std::string part; - while(std::getline(parts, part, ',')) { - X.emplace_back((bool)std::stol(part)); - } - } - in2.close(); - - design_matrix = std::move(vec_to_sparse_matrix(X, n_obs, n_snps)); - } - - MPI_Bcast(&n_snps, 1, MPI_UNSIGNED_LONG_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast(&n_obs, 1, MPI_UNSIGNED_LONG_LONG, 0, MPI_COMM_WORLD); - - RATEd res; - int fullrank = false; - if (rank == 0) { - fullrank = args.value("fullrank"); - } - MPI_Bcast(&fullrank, 1, MPI_INT, 0, MPI_COMM_WORLD); - - MPI_Barrier(MPI_COMM_WORLD); // Wait for all ranks - if (fullrank) { - res = RATE_fullrank(f_draws_mat, design_matrix, n_snps); - } else { - size_t svd_rank; - double prop_var; - if (rank == 0) { - svd_rank = args.value("low-rank"); - svd_rank = svd_rank == 0 ? std::min(n_obs, n_snps) : svd_rank; - - prop_var = args.value("prop-var"); - } - MPI_Bcast(&svd_rank, 1, MPI_UNSIGNED_LONG_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast(&prop_var, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); - - res = RATE_lowrank_mpi(f_draws_mat, design_matrix, n_snps, svd_rank, prop_var); - } - - if (rank == 0) { - std::cout << "#ESS: " << res.ESS << '\n'; - std::cout << "#Delta: " << res.Delta << '\n'; - std::cout << "#snp_id\tRATE\tKLD\n"; - for (size_t i = 0; i < n_snps; ++i) { - std::cout << i << '\t' << res.RATE[i] << '\t' << res.KLD[i] << '\n'; - } - std::cout << std::endl; - } - - rc = MPI_Finalize(); - - return 0; -} From 15ccef45c05653b4ec60e6da74fd36a9b3c55453 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tommi=20M=C3=A4klin?= Date: Fri, 27 Oct 2023 18:55:40 +0300 Subject: [PATCH 02/11] Reve MPI build config. --- CMakeLists.txt | 41 ---------------------------- config/cpprate_mpi_config.hpp.in | 46 -------------------------------- 2 files changed, 87 deletions(-) delete mode 100644 config/cpprate_mpi_config.hpp.in diff --git a/CMakeLists.txt b/CMakeLists.txt index c83b055..ec1969d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -52,23 +52,6 @@ else() endif() configure_file(${CMAKE_CURRENT_SOURCE_DIR}/config/cpprate_blas_config.hpp.in ${CMAKE_CURRENT_BINARY_DIR}/include/cpprate_blas_config.hpp @ONLY) -### MPI -if (CMAKE_ENABLE_MPI_SUPPORT) - find_package(MPI REQUIRED) - set(CPPRATE_MPI_SUPPORT 1) - message(STATUS "Found MPI includes in: ${MPI_C_INCLUDE_DIRS}") - include_directories(MPI_C_INCLUDE_DIRS) - if (CMAKE_MPI_MAX_PROCESSES) - set(CPPRATE_MPI_MAX_PROCESSES ${CMAKE_MPI_MAX_PROCESSES}) - else() - set(CPPRATE_MPI_MAX_PROCESSES 1024) - endif() -else() - set(CPPRATE_MPI_SUPPORT 0) -endif() -## Configure MPI if it's supported on the system. -configure_file(${CMAKE_CURRENT_SOURCE_DIR}/config/cpprate_mpi_config.hpp.in ${CMAKE_CURRENT_BINARY_DIR}/include/cpprate_mpi_config.hpp @ONLY) - ## eigen if (DEFINED CMAKE_EIGEN_HEADERS) message(STATUS "eigen headers provided in: ${CMAKE_EIGEN_HEADERS}") @@ -132,9 +115,6 @@ if(CMAKE_BUILD_TESTS) if (OPENMP_FOUND) target_link_libraries(runTests OpenMP::OpenMP_CXX) endif() - if (CMAKE_ENABLE_MPI_SUPPORT) - target_link_libraries(runTests MPI::MPI_CXX) - endif() if (BLAS_FOUND) target_link_libraries(runTests ${BLAS_LIBRARIES}) endif() @@ -164,20 +144,8 @@ include_directories(${CMAKE_BXZSTR_HEADERS}) if(CMAKE_BUILD_EXECUTABLE) set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/bin) add_executable(cpprate ${CMAKE_CURRENT_SOURCE_DIR}/src/cpprate.cpp) - if (CMAKE_ENABLE_MPI_SUPPORT) - add_executable(cpprate-mpi ${CMAKE_CURRENT_SOURCE_DIR}/src/cpprate_mpi.cpp) - endif() if (OPENMP_FOUND) target_link_libraries(cpprate OpenMP::OpenMP_CXX) - if (CMAKE_ENABLE_MPI_SUPPORT) - target_link_libraries(cpprate-mpi OpenMP::OpenMP_CXX) - endif() - endif() - if (CMAKE_ENABLE_MPI_SUPPORT) - target_link_libraries(cpprate-mpi MPI::MPI_CXX) - if (BLAS_FOUND) - target_link_libraries(cpprate-mpi ${BLAS_LIBRARIES}) - endif() endif() if (BLAS_FOUND) target_link_libraries(cpprate ${BLAS_LIBRARIES}) @@ -187,25 +155,16 @@ if(CMAKE_BUILD_EXECUTABLE) if (BZIP2_FOUND) include_directories(${BZIP2_INCLUDE_DIRS}) target_link_libraries(cpprate ${BZIP2_LIBRARIES}) - if (CMAKE_ENABLE_MPI_SUPPORT) - target_link_libraries(cpprate-mpi ${BZIP2_LIBRARIES}) - endif() endif() find_package(LibLZMA) if (LIBLZMA_FOUND) include_directories(${LIBLZMA_INCLUDE_DIRS}) target_link_libraries(cpprate ${LIBLZMA_LIBRARIES}) - if (CMAKE_ENABLE_MPI_SUPPORT) - target_link_libraries(cpprate-mpi ${LIBLZMA_LIBRARIES}) - endif() endif() find_package(ZLIB) if (ZLIB_FOUND) include_directories(${ZLIB_INCLUDE_DIRS}) target_link_libraries(cpprate ${ZLIB_LIBRARIES}) - if (CMAKE_ENABLE_MPI_SUPPORT) - target_link_libraries(cpprate-mpi ${ZLIB_LIBRARIES}) - endif() endif() endif() diff --git a/config/cpprate_mpi_config.hpp.in b/config/cpprate_mpi_config.hpp.in deleted file mode 100644 index d581831..0000000 --- a/config/cpprate_mpi_config.hpp.in +++ /dev/null @@ -1,46 +0,0 @@ -// cpprate: Variable Selection in Black Box Methods with RelATive cEntrality (RATE) Measures -// https://github.com/tmaklin/cpprate -// Copyright (c) 2023 Tommi Mäklin (tommi@maklin.fi) -// -// BSD-3-Clause license -// -// Redistribution and use in source and binary forms, with or without -// modification, are permitted provided that the following conditions are -// met: -// -// (1) Redistributions of source code must retain the above copyright -// notice, this list of conditions and the following disclaimer. -// -// (2) Redistributions in binary form must reproduce the above copyright -// notice, this list of conditions and the following disclaimer in -// the documentation and/or other materials provided with the -// distribution. -// -// (3)The name of the author may not be used to -// endorse or promote products derived from this software without -// specific prior written permission. -// -// THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR -// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED -// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -// DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, -// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) -// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING -// IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -// POSSIBILITY OF SUCH DAMAGE. -// -#ifndef CPPRATE_MPI_CONFIG_HPP -#define CPPRATE_MPI_CONFIG_HPP - -#define CPPRATE_MPI_SUPPORT @CPPRATE_MPI_SUPPORT@ - -#if defined(CPPRATE_MPI_SUPPORT) && (CPPRATE_MPI_SUPPORT) == 1 -#define OMPI_SKIP_MPICXX 1 // See https://github.com/open-mpi/ompi/issues/5157 -#include -#endif - - -#endif From 57f48e27b9a88cfe617c0a90e383c0b3ee9b97e4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tommi=20M=C3=A4klin?= Date: Fri, 23 Feb 2024 16:11:32 +0200 Subject: [PATCH 03/11] Update readme with links to prebuilt binaries. --- README.md | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 6f14d7f..5838150 100644 --- a/README.md +++ b/README.md @@ -2,9 +2,9 @@ cpprate is a reimplementation of https://github.com/lorinanthony/rate in C++. # Installation -This is a header-only library so simply include "CppRateRes.hpp" or -"CppRateRes_mpi.hpp" in your project. An optional command-line -executable is also provided, see instructions below for compiling. +## Prebuilt binaries +Prebuilt binaries are available from the [releases page](https://github.com/tmaklin/cpprate/releases) for linux\_x86-64, macOS\_x86-64, and macOS\_arm64. + ## Compiling the cpprate executable from source ### Dependencies - cmake >= v3.1 @@ -19,6 +19,10 @@ cmake -DCMAKE_BUILD_EXECUTABLE=1 .. make -j ``` +## API +This is a header-only library so simply include "CppRateRes.hpp" or +"CppRateRes_mpi.hpp" header in your project. + # Usage Input files can optionally be compressed with gzip/bzip2/xz. The format is detected automatically. ## Nonlinear coefficients From 0942d80191af473048ea4ec4a4537b239550336e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tommi=20M=C3=A4klin?= Date: Fri, 8 Mar 2024 07:44:50 +0000 Subject: [PATCH 04/11] Create make_release.yml --- .github/workflows/make_release.yml | 127 +++++++++++++++++++++++++++++ 1 file changed, 127 insertions(+) create mode 100644 .github/workflows/make_release.yml diff --git a/.github/workflows/make_release.yml b/.github/workflows/make_release.yml new file mode 100644 index 0000000..9843f09 --- /dev/null +++ b/.github/workflows/make_release.yml @@ -0,0 +1,127 @@ +name: Build cpprate binaries +on: + push: + tags: + - "v*.*.*" + branches: + - cpprate-release-testing + +jobs: + build_linux-x86_64: + runs-on: ubuntu-latest + container: phusion/holy-build-box-64:3.0.2 + steps: + - name: Install wget + id: install-wget + run: yum install -y wget + + - name: Create io directory + id: mkdir-io + run: mkdir /io && cd /io + + - name: Download build script + id: dl-build-script + run: wget https://raw.githubusercontent.com/tmaklin/biobins/master/linux/cpprate/build.sh + + - name: Compile binary in Holy Build Box container + id: compile-in-container + run: chmod +x build.sh && ./build.sh ${{ github.ref_name }} + + - name: Upload linux-x86_64 binary + if: success() + uses: actions/upload-artifact@v3 + with: + name: cpprate-${{ github.ref_name }}-x86_64-redhat-linux + path: /io/cpprate-${{ github.ref_name }}-x86_64-redhat-linux.tar.gz + + build_macOS-x86_64: + runs-on: ubuntu-latest + container: ghcr.io/shepherdjerred/macos-cross-compiler:latest + steps: + - name: Install wget + id: install-wget + run: apt install -y wget + + - name: Create io directory + id: mkdir-io + run: mkdir /io && cd /io + + - name: Download toolchain file + id: dl-toolchain-file + run: wget https://raw.githubusercontent.com/tmaklin/biobins/master/macOS/x86-64-toolchain.cmake && cp x86-64-toolchain.cmake /io/x86-64-toolchain.cmake && cp x86-64-toolchain.cmake /x86-64-toolchain.cmake + + - name: Download build script + id: dl-build-script + run: wget https://raw.githubusercontent.com/tmaklin/biobins/master/linux/cpprate/build.sh + + - name: Compile binary in macOS Cross Compiler container + id: compile-in-container + run: chmod +x build.sh && ./build.sh ${{ github.ref_name }} x86-64 + + - name: Upload macOS-x86_64 binary + if: success() + uses: actions/upload-artifact@v3 + with: + name: cpprate-${{ github.ref_name }}-x86_64-apple-darwin22 + path: /io/cpprate-${{ github.ref_name }}-x86_64-apple-darwin22.tar.gz + + build_macOS-arm64: + runs-on: ubuntu-latest + container: ghcr.io/shepherdjerred/macos-cross-compiler:latest + steps: + - name: Install wget + id: install-wget + run: apt install -y wget + + - name: Create io directory + id: mkdir-io + run: mkdir /io && cd /io + + - name: Download toolchain file + id: dl-toolchain-file + run: wget https://raw.githubusercontent.com/tmaklin/biobins/master/macOS/arm64-toolchain.cmake && cp arm64-toolchain.cmake /io/arm64-toolchain.cmake && cp arm64-toolchain.cmake /arm64-toolchain.cmake + + - name: Download build script + id: dl-build-script + run: wget https://raw.githubusercontent.com/tmaklin/biobins/master/linux/cpprate/build.sh + + - name: Compile binary in macOS Cross Compiler container + id: compile-in-container + run: chmod +x build.sh && ./build.sh ${{ github.ref_name }} arm64 + + - name: Upload macOS-arm64 binary + if: success() + uses: actions/upload-artifact@v3 + with: + name: cpprate-${{ github.ref_name }}-arm64-apple-darwin22 + path: /io/cpprate-${{ github.ref_name }}-arm64-apple-darwin22.tar.gz + + create-release: + runs-on: ubuntu-latest + + needs: [ build_linux-x86_64, build_macOS-x86_64, build_macOS-arm64 ] + + steps: + - uses: actions/checkout@v2 + + - uses: actions/download-artifact@v2 + with: + path: build + + - name: Organise files + shell: bash + run: | + cp build/cpprate-${{ github.ref_name }}-arm64-apple-darwin22/cpprate-${{ github.ref_name }}-arm64-apple-darwin22.tar.gz . + cp build/cpprate-${{ github.ref_name }}-x86_64-apple-darwin22/cpprate-${{ github.ref_name }}-x86_64-apple-darwin22.tar.gz . + cp build/cpprate-${{ github.ref_name }}-x86_64-redhat-linux/cpprate-${{ github.ref_name }}-x86_64-redhat-linux.tar.gz . + - name: Create release + id: create_release + uses: softprops/action-gh-release@v1 + with: + name: Release ${{ github.ref_name }} + draft: true + prerelease: false + fail_on_unmatched_files: true + generate_release_notes: true + files: | + cpprate-*.tar.gz From 1130baa05a29b7b11e27485173db10c4331009ba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tommi=20M=C3=A4klin?= Date: Fri, 8 Mar 2024 08:11:32 +0000 Subject: [PATCH 05/11] Update make_release.yml Fix linux->macOS in downloading build scripts for mac. --- .github/workflows/make_release.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/make_release.yml b/.github/workflows/make_release.yml index 9843f09..f21511d 100644 --- a/.github/workflows/make_release.yml +++ b/.github/workflows/make_release.yml @@ -52,7 +52,7 @@ jobs: - name: Download build script id: dl-build-script - run: wget https://raw.githubusercontent.com/tmaklin/biobins/master/linux/cpprate/build.sh + run: wget https://raw.githubusercontent.com/tmaklin/biobins/master/macOS/cpprate/build.sh - name: Compile binary in macOS Cross Compiler container id: compile-in-container @@ -83,7 +83,7 @@ jobs: - name: Download build script id: dl-build-script - run: wget https://raw.githubusercontent.com/tmaklin/biobins/master/linux/cpprate/build.sh + run: wget https://raw.githubusercontent.com/tmaklin/biobins/master/macOS/cpprate/build.sh - name: Compile binary in macOS Cross Compiler container id: compile-in-container From 4de8a125a86d5e1af676e9c3062754f450712a07 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tommi=20M=C3=A4klin?= Date: Fri, 8 Mar 2024 11:26:20 +0200 Subject: [PATCH 06/11] Exit with an error message if input files don't exist (resolves #10) --- src/cpprate.cpp | 43 +++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 41 insertions(+), 2 deletions(-) diff --git a/src/cpprate.cpp b/src/cpprate.cpp index 811598a..6428f2f 100644 --- a/src/cpprate.cpp +++ b/src/cpprate.cpp @@ -39,6 +39,7 @@ #include #include #include +#include #include "bxzstr.hpp" #include "BS_thread_pool.hpp" @@ -51,6 +52,14 @@ bool CmdOptionPresent(char **begin, char **end, const std::string &option) { return (std::find(begin, end, option) != end); } +bool file_exists(const std::string &file_path) { + std::filesystem::path check_file{ file_path }; + if (!std::filesystem::exists(check_file)) { + throw std::runtime_error(file_path + " does not exist."); + } + return true; +} + bool parse_args(int argc, char* argv[], cxxargs::Arguments &args) { args.add_short_argument('f', "f-draws file (comma separated)"); args.add_short_argument('x', "design matrix (comma separated)"); @@ -203,7 +212,16 @@ int main(int argc, char* argv[]) { Eigen::MatrixXd posterior_draws; size_t n_draws = 0; size_t n_obs_draws = 0; - bxz::ifstream posterior_draws_in((from_beta_draws ? args.value("beta-draws") : args.value('f'))); + + std::string posterior_draws_path = (from_beta_draws ? args.value("beta-draws") : args.value('f')); + try { + file_exists(posterior_draws_path); + } catch (std::exception &e) { + std::cerr << std::string("cpprate: ") + e.what() << std::endl; + return 1; + } + + bxz::ifstream posterior_draws_in(posterior_draws_path); read_posterior_draws(&posterior_draws_in, &n_draws, &n_obs_draws, &posterior_draws); posterior_draws_in.close(); @@ -212,6 +230,12 @@ int main(int argc, char* argv[]) { // Read in the design matrix size_t n_snps = 0; size_t n_obs_X = 0; + try { + file_exists(args.value('x')); + } catch (std::exception &e) { + std::cerr << std::string("cpprate: ") + e.what() << std::endl; + return 1; + } bxz::ifstream design_matrix_in(args.value('x')); // Next call will overwrite the f_draws currently stored in posterior_draws // with beta draws = f_draws * pseudoinverse(design_matrix) @@ -235,7 +259,16 @@ int main(int argc, char* argv[]) { Eigen::MatrixXd posterior_draws; size_t n_draws = 0; size_t n_obs_draws = 0; - bxz::ifstream posterior_draws_in((from_beta_draws ? args.value("beta-draws") : args.value('f'))); + + std::string posterior_draws_path = (from_beta_draws ? args.value("beta-draws") : args.value('f')); + try { + file_exists(posterior_draws_path); + } catch (std::exception &e) { + std::cerr << std::string("cpprate: ") + e.what() << std::endl; + return 1; + } + + bxz::ifstream posterior_draws_in(posterior_draws_path); read_posterior_draws(&posterior_draws_in, &n_draws, &n_obs_draws, &posterior_draws); posterior_draws_in.close(); @@ -244,6 +277,12 @@ int main(int argc, char* argv[]) { Eigen::MatrixXd svd_design_matrix_U; size_t n_obs; + try { + file_exists(args.value('x')); + } catch (std::exception &e) { + std::cerr << std::string("cpprate: ") + e.what() << std::endl; + return 1; + } if (from_beta_draws) { const Eigen::SparseMatrix &design_matrix = read_decomposition(args.value('x'), low_rank_rank, args.value("prop-var"), &n_snps, &n_obs, &svd_design_matrix_U, static_cast(cov_beta_ptr.get())->get_svd_V_p()).transpose(); posterior_draws *= design_matrix ; From 01745dadc2e12fda206ee0f276e368c1f1fbd5cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tommi=20M=C3=A4klin?= Date: Fri, 8 Mar 2024 11:45:51 +0200 Subject: [PATCH 07/11] Catch and print exceptions from the cpprate executable (resolves #11) --- src/cpprate.cpp | 124 +++++++++++++++++++++++++++++++----------------- 1 file changed, 80 insertions(+), 44 deletions(-) diff --git a/src/cpprate.cpp b/src/cpprate.cpp index 6428f2f..49a7188 100644 --- a/src/cpprate.cpp +++ b/src/cpprate.cpp @@ -221,9 +221,14 @@ int main(int argc, char* argv[]) { return 1; } - bxz::ifstream posterior_draws_in(posterior_draws_path); - read_posterior_draws(&posterior_draws_in, &n_draws, &n_obs_draws, &posterior_draws); - posterior_draws_in.close(); + try { + bxz::ifstream posterior_draws_in(posterior_draws_path); + read_posterior_draws(&posterior_draws_in, &n_draws, &n_obs_draws, &posterior_draws); + posterior_draws_in.close(); + } catch (std::exception &e) { + std::cerr << std::string("cpprate: Error in reading posterior draws: ") + e.what() << std::endl; + return 1; + } // If running fullrank algorithm on f draws we also need to read the design matrix if (!from_beta_draws) { @@ -236,20 +241,31 @@ int main(int argc, char* argv[]) { std::cerr << std::string("cpprate: ") + e.what() << std::endl; return 1; } - bxz::ifstream design_matrix_in(args.value('x')); - // Next call will overwrite the f_draws currently stored in posterior_draws - // with beta draws = f_draws * pseudoinverse(design_matrix) - read_nonlinear_coefficients(&design_matrix_in, &n_snps, &n_obs_X, &posterior_draws); + try { + bxz::ifstream design_matrix_in(args.value('x')); + // Next call will overwrite the f_draws currently stored in posterior_draws + // with beta draws = f_draws * pseudoinverse(design_matrix) + read_nonlinear_coefficients(&design_matrix_in, &n_snps, &n_obs_X, &posterior_draws); + } catch (std::exception &e) { + std::cerr << std::string("cpprate: Error in computing `f_draws * pseudoinverse(design_matrix)`: ") + e.what() << std::endl; + return 1; + } // Check that the input dimensions are correct if (n_obs_X != n_obs_draws) { - throw std::runtime_error("Number of rows in file " + args.value('x') + " (" + std::to_string(n_obs_X) + ") does not match number of columns in file " + args.value('f') + " (" + std::to_string(n_obs_draws) + ")."); + std::cerr << std::string("cpprate: Number of rows in file " + args.value('x') + " (" + std::to_string(n_obs_X) + ") does not match number of columns in file " + args.value('f') + " (" + std::to_string(n_obs_draws) + ").") << std::endl; + return 1; } } - cov_beta_ptr.reset(new FullrankCovMat()); - static_cast(cov_beta_ptr.get())->fill(posterior_draws); - col_means_beta = std::move(col_means2(posterior_draws)); + try { + cov_beta_ptr.reset(new FullrankCovMat()); + static_cast(cov_beta_ptr.get())->fill(posterior_draws); + col_means_beta = std::move(col_means2(posterior_draws)); + } catch (std::exception &e) { + std::cerr << std::string("cpprate: Error in computing fullrank covariance matrix: ") + e.what() << std::endl; + return 1; + } } // Otherwise running lowrank algorithm @@ -268,45 +284,60 @@ int main(int argc, char* argv[]) { return 1; } - bxz::ifstream posterior_draws_in(posterior_draws_path); - read_posterior_draws(&posterior_draws_in, &n_draws, &n_obs_draws, &posterior_draws); - posterior_draws_in.close(); - - cov_beta_ptr.reset(new LowrankCovMat()); - // Decompose the design matrix - Eigen::MatrixXd svd_design_matrix_U; - size_t n_obs; - try { - file_exists(args.value('x')); + bxz::ifstream posterior_draws_in(posterior_draws_path); + read_posterior_draws(&posterior_draws_in, &n_draws, &n_obs_draws, &posterior_draws); + posterior_draws_in.close(); } catch (std::exception &e) { - std::cerr << std::string("cpprate: ") + e.what() << std::endl; + std::cerr << std::string("cpprate: Error in reading posterior draws: ") + e.what() << std::endl; return 1; } - if (from_beta_draws) { - const Eigen::SparseMatrix &design_matrix = read_decomposition(args.value('x'), low_rank_rank, args.value("prop-var"), &n_snps, &n_obs, &svd_design_matrix_U, static_cast(cov_beta_ptr.get())->get_svd_V_p()).transpose(); - posterior_draws *= design_matrix ; - } else { - read_decomposition(args.value('x'), low_rank_rank, args.value("prop-var"), &n_snps, &n_obs, &svd_design_matrix_U, static_cast(cov_beta_ptr.get())->get_svd_V_p()) ; - } - // Check that the input dimensions are correct - if (n_snps != n_obs_draws && from_beta_draws) { - throw std::runtime_error("Number of columns in file " + args.value('x') + " (" + std::to_string(n_snps) + ") does not match number of columns in file " + args.value("beta-draws") + " (" + std::to_string(n_obs_draws) + ")."); - } else if (n_obs != n_obs_draws && !from_beta_draws) { - throw std::runtime_error("Number of rows in file " + args.value('x') + " (" + std::to_string(n_obs) + ") does not match number of columns in file " + args.value('f') + " (" + std::to_string(n_obs_draws) + ")."); - } + try { + cov_beta_ptr.reset(new LowrankCovMat()); + // Decompose the design matrix + Eigen::MatrixXd svd_design_matrix_U; + size_t n_obs; - // If lowrank decomposition rank was supplied set it to minimum of the dimension of design_matrix - low_rank_rank = low_rank_rank == 0 ? std::min(n_snps, n_obs) : low_rank_rank; + try { + file_exists(args.value('x')); + } catch (std::exception &e) { + std::cerr << std::string("cpprate: ") + e.what() << std::endl; + return 1; + } + if (from_beta_draws) { + const Eigen::SparseMatrix &design_matrix = read_decomposition(args.value('x'), low_rank_rank, args.value("prop-var"), &n_snps, &n_obs, &svd_design_matrix_U, static_cast(cov_beta_ptr.get())->get_svd_V_p()).transpose(); + posterior_draws *= design_matrix ; + } else { + read_decomposition(args.value('x'), low_rank_rank, args.value("prop-var"), &n_snps, &n_obs, &svd_design_matrix_U, static_cast(cov_beta_ptr.get())->get_svd_V_p()) ; + } - static_cast(cov_beta_ptr.get())->construct(project_f_draws(posterior_draws, svd_design_matrix_U), low_rank_rank); - col_means_beta = std::move(approximate_beta_means(posterior_draws, svd_design_matrix_U, static_cast(cov_beta_ptr.get())->get_svd_V())); + // Check that the input dimensions are correct + if (n_snps != n_obs_draws && from_beta_draws) { + throw std::runtime_error("Number of columns in file " + args.value('x') + " (" + std::to_string(n_snps) + ") does not match number of columns in file " + args.value("beta-draws") + " (" + std::to_string(n_obs_draws) + ")."); + } else if (n_obs != n_obs_draws && !from_beta_draws) { + throw std::runtime_error("Number of rows in file " + args.value('x') + " (" + std::to_string(n_obs) + ") does not match number of columns in file " + args.value('f') + " (" + std::to_string(n_obs_draws) + ")."); + } + + // If lowrank decomposition rank was supplied set it to minimum of the dimension of design_matrix + low_rank_rank = low_rank_rank == 0 ? std::min(n_snps, n_obs) : low_rank_rank; + + static_cast(cov_beta_ptr.get())->construct(project_f_draws(posterior_draws, svd_design_matrix_U), low_rank_rank); + col_means_beta = std::move(approximate_beta_means(posterior_draws, svd_design_matrix_U, static_cast(cov_beta_ptr.get())->get_svd_V())); + } catch (std::exception &e) { + std::cerr << std::string("cpprate: Error in computing lowrank covariance matrix: ") + e.what() << std::endl; + return 1; + } } if (!run_fullrank) { - static_cast(cov_beta_ptr.get())->logarithmize_lambda(); - static_cast(cov_beta_ptr.get())->logarithmize_svd_V(); + try { + static_cast(cov_beta_ptr.get())->logarithmize_lambda(); + static_cast(cov_beta_ptr.get())->logarithmize_svd_V(); + } catch (std::exception &e) { + std::cerr << std::string("cpprate: Error in logarithmizing lowrank covariance matrix `lambda` and `svd_V`: ") + e.what() << std::endl; + return 1; + } } BS::thread_pool pool(n_ranks); @@ -334,10 +365,15 @@ int main(int argc, char* argv[]) { size_t global_index = id_start; std::vector log_KLD_global(n_snps, -36.84136); for (size_t thread_id = 0; thread_id < n_ranks; ++thread_id) { - const std::vector &log_KLD_local = thread_futures[thread_id].get(); - for (size_t i = 0; i < log_KLD_local.size(); ++i) { - log_KLD_global[global_index] = log_KLD_local[i]; - ++global_index; + try { + const std::vector &log_KLD_local = thread_futures[thread_id].get(); + for (size_t i = 0; i < log_KLD_local.size(); ++i) { + log_KLD_global[global_index] = log_KLD_local[i]; + ++global_index; + } + } catch (std::exception &e) { + std::cerr << std::string("cpprate: thread ") + std::to_string(thread_id) + std::string(": Error in computing KLD for regression coefficent: ") + std::to_string(global_index) + std::string(": ") + e.what() << std::endl; + return 1; } } From fda6983d344daa47f3ee2a2970c280a3ea1e0faf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tommi=20M=C3=A4klin?= Date: Fri, 8 Mar 2024 11:59:40 +0200 Subject: [PATCH 08/11] Add more checks for input dimensions (resolves #7). --- src/cpprate.cpp | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/src/cpprate.cpp b/src/cpprate.cpp index 49a7188..a222a7d 100644 --- a/src/cpprate.cpp +++ b/src/cpprate.cpp @@ -230,6 +230,17 @@ int main(int argc, char* argv[]) { return 1; } + if (n_obs_draws <= 1) { + std::cerr << + "cpprate: " + + (from_beta_draws ? + "Refusing to compute RATE for 1 regression coefficient (result is not meaningful)" + : + "Cannot compute RATE with only 1 observation in " + posterior_draws_path) + << std::endl; + return 1; + } + // If running fullrank algorithm on f draws we also need to read the design matrix if (!from_beta_draws) { // Read in the design matrix @@ -293,6 +304,17 @@ int main(int argc, char* argv[]) { return 1; } + if (n_obs_draws <= 1) { + std::cerr << + "cpprate: " + + (from_beta_draws ? + "Refusing to compute RATE for 1 regression coefficient (result is not meaningful)" + : + "Cannot compute RATE with only 1 observation in " + posterior_draws_path) + << std::endl; + return 1; + } + try { cov_beta_ptr.reset(new LowrankCovMat()); // Decompose the design matrix @@ -340,6 +362,11 @@ int main(int argc, char* argv[]) { } } + if (n_snps <= 1) { + std::cerr << "cpprate: Refusing to compute RATE for 1 regression coefficient (result is not meaningful)" << std::endl; + return 1; + } + BS::thread_pool pool(n_ranks); std::vector>> thread_futures; From ccbd828a907f1af63bb1b73fe721a46b46840611 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tommi=20M=C3=A4klin?= Date: Fri, 8 Mar 2024 12:31:15 +0200 Subject: [PATCH 09/11] Add progress indicator bar to run_RATE. --- include/CppRateRes.hpp | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/include/CppRateRes.hpp b/include/CppRateRes.hpp index 6f25f55..c30073b 100644 --- a/include/CppRateRes.hpp +++ b/include/CppRateRes.hpp @@ -56,6 +56,18 @@ #include "CovarianceMatrix.hpp" #include "RATE_res.hpp" +// Progress indicator +#define cpprate_PBSTR "||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||" +#define cpprate_PBWIDTH 60 + +inline void print_progress(double percentage) { + int val = (int) (percentage * 100); + int lpad = (int) (percentage * cpprate_PBWIDTH); + int rpad = cpprate_PBWIDTH - lpad; + printf("\r%3d%% [%.*s%*s]", val, lpad, cpprate_PBSTR, rpad, ""); + fflush(stderr); +} + inline void decompose_design_matrix(const Eigen::SparseMatrix &design_matrix, const size_t svd_rank, const double prop_var, Eigen::MatrixXd *u, Eigen::MatrixXd *v) { // Calculate the singular value decomposition of `design_matrix` @@ -240,7 +252,7 @@ Eigen::SparseMatrix vec_to_sparse_matrix(const std::vector &vec, const siz inline std::vector run_RATE(const std::vector &col_means_beta, const std::shared_ptr &cov_beta_ptr, const std::vector &ids_to_test, const size_t id_start, const size_t id_end, const size_t n_snps, - const size_t n_threads = 0) { + const size_t n_threads = 0, bool progress = false) { std::vector log_KLD(n_snps, -36.84136); // log(1e-16) = -36.84136 #if defined(CPPRATE_OPENMP_SUPPORT) && (CPPRATE_OPENMP_SUPPORT) == 1 if (n_threads > 0) { @@ -258,6 +270,13 @@ inline std::vector run_RATE(const std::vector &col_means_beta, c const double log_m = std::log(std::abs(col_means_beta[snp_id]) + 1e-16); log_KLD[log_KLD_index] = log_m + log_m + log_alpha + std::log(0.5); ++log_KLD_index; + + if (progress) { + print_progress((double)i/end); + } + } + if (progress) { + std::cerr << std::endl; } return log_KLD; } From a5526d9872725681a18898cfac2b3ac4cbeb0667 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tommi=20M=C3=A4klin?= Date: Fri, 8 Mar 2024 12:33:00 +0200 Subject: [PATCH 10/11] Print progress from first rank with --print-progress (resolves #12) --- src/cpprate.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/cpprate.cpp b/src/cpprate.cpp index a222a7d..e161c0e 100644 --- a/src/cpprate.cpp +++ b/src/cpprate.cpp @@ -72,6 +72,7 @@ bool parse_args(int argc, char* argv[], cxxargs::Arguments &args) { args.add_long_argument("prop-var", "Proportion of variance to explain in lowrank factorization (default: 100%)", 1.1); args.add_long_argument("low-rank", "Rank of the low-rank factorization (default: min(design_matrix.rows(), design_matrix.cols()))", 0); args.add_long_argument("fullrank", "Run fullrank algorithm (default: false)", false); + args.add_long_argument("print-progress", "Print progress bar from the first rank (default: false)", false); args.add_long_argument("help", "Print the help message.", false); if (CmdOptionPresent(argv, argv+argc, "--help")) { std::cout << "\n" + args.help() << '\n' << std::endl; @@ -386,7 +387,7 @@ int main(int argc, char* argv[]) { n_processed += snps_per_rank[thread_id]; size_t thread_n_snps = thread_stop_at - thread_start_at; thread_futures.emplace_back(pool.submit(run_RATE, col_means_beta, cov_beta_ptr, - args.value>("ids-to-test"), thread_start_at, thread_stop_at, thread_n_snps, n_threads)); + args.value>("ids-to-test"), thread_start_at, thread_stop_at, thread_n_snps, n_threads, thread_id == 0 && args.value("print-progress"))); } size_t global_index = id_start; From 722d7b0ced095b35c5b78cc51ded2debef88bc95 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tommi=20M=C3=A4klin?= Date: Fri, 8 Mar 2024 10:53:02 +0000 Subject: [PATCH 11/11] Update make_release.yml add permissions: contents: write --- .github/workflows/make_release.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/make_release.yml b/.github/workflows/make_release.yml index f21511d..6f5308f 100644 --- a/.github/workflows/make_release.yml +++ b/.github/workflows/make_release.yml @@ -98,6 +98,8 @@ jobs: create-release: runs-on: ubuntu-latest + permissions: + contents: write needs: [ build_linux-x86_64, build_macOS-x86_64, build_macOS-arm64 ]