From 017cd2e25d2120beccb1450b08870df5cf975f35 Mon Sep 17 00:00:00 2001 From: ZhiliZheng Date: Thu, 23 Feb 2023 03:16:56 +1000 Subject: [PATCH] add the outfreq, add sigmaSq, bug fix --- DESCRIPTION | 4 ++-- R/RcppExports.R | 4 ++-- R/sbr_eigen.r | 6 +++--- src/AnnoProb.cpp | 21 +++++++++++++++++---- src/AnnoProb.hpp | 3 ++- src/RcppExports.cpp | 10 ++++++---- src/SBayesRC.cpp | 9 +++++++-- src/SBayesRC.hpp | 5 ++++- src/dist.cpp | 14 +++++++------- src/sbr_eigen_c.cpp | 9 ++++++--- 10 files changed, 56 insertions(+), 29 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index a9bd4c3..3159ea7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: SBayesRC Type: Package Title: Run SBayesRC -Version: 0.1.4 -Date: 2023-01-02 +Version: 0.1.5 +Date: 2023-02-22 Description: SBayesRC Authors@R: c( person("Zhili","Zheng", role=c("aut","cre"), email="zhilizheng@outlook.com"), diff --git a/R/RcppExports.R b/R/RcppExports.R index a75864f..7113fde 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -14,7 +14,7 @@ impGa <- function(ldm, z, typedIndex, m, cutThresh = 1, diag_mod = 0.1) { .Call(`_SBayesRC_impGa`, ldm, z, typedIndex, m, cutThresh, diag_mod) } -sbayesr_eigen_joint_annot <- function(niter, burn, bhat, numAnno, annoStrs, mldmDir, vary, blkN, cgamma, startPi, starth2 = 0.01, cutThresh = 1, bOrigin = FALSE, outPrefix = "", samVe = "fixVe", resam_thresh = 1.1, bOutDetail = FALSE) { - .Call(`_SBayesRC_sbayesr_eigen_joint_annot`, niter, burn, bhat, numAnno, annoStrs, mldmDir, vary, blkN, cgamma, startPi, starth2, cutThresh, bOrigin, outPrefix, samVe, resam_thresh, bOutDetail) +sbayesr_eigen_joint_annot <- function(niter, burn, bhat, numAnno, annoStrs, mldmDir, vary, blkN, cgamma, startPi, starth2 = 0.01, cutThresh = 1, bOrigin = FALSE, outPrefix = "", samVe = "fixVe", resam_thresh = 1.1, bOutDetail = FALSE, outFreq = 10L, initAnnoSS = 1.0) { + .Call(`_SBayesRC_sbayesr_eigen_joint_annot`, niter, burn, bhat, numAnno, annoStrs, mldmDir, vary, blkN, cgamma, startPi, starth2, cutThresh, bOrigin, outPrefix, samVe, resam_thresh, bOutDetail, outFreq, initAnnoSS) } diff --git a/R/sbr_eigen.r b/R/sbr_eigen.r index 61be3d0..4bc8106 100644 --- a/R/sbr_eigen.r +++ b/R/sbr_eigen.r @@ -21,7 +21,7 @@ sbayesrc = function(file_summary, ld_folder, file_out, thresh=0.995, niter=3000, burn=1000, fileAnnot="", method="sbr_ori", sSamVe="allMixVe", twopq="nbsq", bOutDetail=FALSE, resam_thresh=1.1, - starth2=0.5, startPi=c(0.990, 0.005, 0.003, 0.001, 0.001), gamma=c(0, 0.001, 0.01, 0.1, 1), seed=22){ + starth2=0.5, startPi=c(0.990, 0.005, 0.003, 0.001, 0.001), gamma=c(0, 0.001, 0.01, 0.1, 1), seed=22, outFreq=10, annoSigmaScale=1.0){ cSamVe = "fixVe" if(sSamVe == "noReSamVe"){ cSamVe = "fixVe" @@ -78,7 +78,7 @@ sbayesrc = function(file_summary, ld_folder, file_out, thresh=0.995, niter=3000, if(file.exists(paste0(outfile, ".txt"))){ message("Don't need to run: ", outfile, " exists") - quit() + return; } if(submethod == "sbc"){ @@ -237,7 +237,7 @@ sbayesrc = function(file_summary, ld_folder, file_out, thresh=0.995, niter=3000, message("The parameter file exists, loading the parameter instead of a re-run: ", outRes) res = readRDS(outRes) }else{ - res = sbayesr_eigen_joint_annot(niter, burn, bhat, numAnno, annoStrings, ld_folder, vary, n, gamma, startPi, starth2, thresh, bOri, outfile, cSamVe, resam_thresh, bOutDetail) + res = sbayesr_eigen_joint_annot(niter, burn, bhat, numAnno, annoStrings, ld_folder, vary, n, gamma, startPi, starth2, thresh, bOri, outfile, cSamVe, resam_thresh, bOutDetail, outFreq, annoSigmaScale) saveRDS(res, file=outRes) } diff --git a/src/AnnoProb.cpp b/src/AnnoProb.cpp index 164e553..201225c 100644 --- a/src/AnnoProb.cpp +++ b/src/AnnoProb.cpp @@ -149,7 +149,7 @@ void AnnoProb::close(){ void AnnoProb::annoSigmaSS_sample() { // init sigmaSS: SetOnes(ssq.size()); const int df = 4; - const int scale = 1; + const int scale = 1; // changed from 1 float dfTilde = df + (numAnno - 1); // exclude the intercept VectorXf scaleTilde = ssq.array() + df * scale; @@ -218,6 +218,7 @@ void AnnoProb::annoEffect_sample_Gibbs(const MatrixXf &z) { // sample latent variables // # pragma omp parallel for schedule(dynamic) + #pragma omp parallel for for (int j = 0; j < numDP; ++j) { if (zi[j]) y[j] = TruncatedNormal::sample_lower_truncated(mean[j], 1.0, 0.0); @@ -301,7 +302,11 @@ void AnnoProb::computePiFromP(const MatrixXf &snpP, MatrixXf &snpPi) { } } -AnnoProb::AnnoProb(string fileAnnot, int numAnno, const VectorXf &Pi, MatrixXf &snpPi, bool bOutDetail) { +VectorXd AnnoProb::getSigmaSq(){ + return sigmaSq.cast(); +} + +AnnoProb::AnnoProb(string fileAnnot, int numAnno, const VectorXf &Pi, MatrixXf &snpPi, bool bOutDetail, double initSS) { this->bOutDetail = bOutDetail; numDist = Pi.size(); @@ -352,7 +357,7 @@ AnnoProb::AnnoProb(string fileAnnot, int numAnno, const VectorXf &Pi, MatrixXf & XPXiSD.resize(numAnno); XPXqiSD.resize(numAnno); - #pragma omp parellel for schedule(dynamic) + #pragma omp parallel for schedule(dynamic) for(int i = 0; i < numAnno; i++){ float XPX = annoMat.col(i).dot(annoMat.col(i)); if(XPX < 1e-6){ @@ -394,7 +399,15 @@ AnnoProb::AnnoProb(string fileAnnot, int numAnno, const VectorXf &Pi, MatrixXf & initP_Pi_annoAlpha(Pi, snpPi); // init sigmaSS - sigmaSq = VectorXf::Ones(numComp); + //sigmaSq = VectorXf::Ones(numComp).array() * (float) initSS; + //sigmaSq = VectorXf::LinSpaced(numComp, 0.1, 1.0) * (float) initSS; + sigmaSq = VectorXf::Ones(numComp) * (float)initSS; + /* + for(int i = 1; i < numComp; i++){ + sigmaSq[i] = sigmaSq[i-1]*2; + } + if(sigmaSq[numComp-1] < 1.0) sigmaSq[numComp-1] = 1.0; + */ //init ssq ssq.resize(numComp); diff --git a/src/AnnoProb.hpp b/src/AnnoProb.hpp index dc8321e..1486d59 100644 --- a/src/AnnoProb.hpp +++ b/src/AnnoProb.hpp @@ -70,7 +70,7 @@ class AnnoProb { void writeHeader(const vector &annoStrs); public: - AnnoProb(string fileAnnot, int numAnno, const VectorXf &Pi, MatrixXf &snpPi, bool bOutDetail); + AnnoProb(string fileAnnot, int numAnno, const VectorXf &Pi, MatrixXf &snpPi, bool bOutDetail, double initSS=1); void samplePi(const MatrixXf &z, MatrixXf &snpPi); void computeEnrichBin(const MatrixXf &vg_snp_comp, float vp); void computeEnrichQt(const VectorXf &vg_snp, float vg); @@ -82,6 +82,7 @@ class AnnoProb { void open(string prefix, const vector &annoStrs); void close(); void output(); + VectorXd getSigmaSq(); }; #endif //SBRC_ANNO_PROB_H diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 8de6c34..5d69eaf 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -42,8 +42,8 @@ BEGIN_RCPP END_RCPP } // sbayesr_eigen_joint_annot -List sbayesr_eigen_joint_annot(int niter, int burn, Eigen::Map bhat, int numAnno, Rcpp::StringVector annoStrs, std::string mldmDir, double vary, Eigen::Map blkN, Eigen::Map cgamma, Eigen::Map startPi, double starth2, double cutThresh, bool bOrigin, std::string outPrefix, std::string samVe, double resam_thresh, bool bOutDetail); -RcppExport SEXP _SBayesRC_sbayesr_eigen_joint_annot(SEXP niterSEXP, SEXP burnSEXP, SEXP bhatSEXP, SEXP numAnnoSEXP, SEXP annoStrsSEXP, SEXP mldmDirSEXP, SEXP varySEXP, SEXP blkNSEXP, SEXP cgammaSEXP, SEXP startPiSEXP, SEXP starth2SEXP, SEXP cutThreshSEXP, SEXP bOriginSEXP, SEXP outPrefixSEXP, SEXP samVeSEXP, SEXP resam_threshSEXP, SEXP bOutDetailSEXP) { +List sbayesr_eigen_joint_annot(int niter, int burn, Eigen::Map bhat, int numAnno, Rcpp::StringVector annoStrs, std::string mldmDir, double vary, Eigen::Map blkN, Eigen::Map cgamma, Eigen::Map startPi, double starth2, double cutThresh, bool bOrigin, std::string outPrefix, std::string samVe, double resam_thresh, bool bOutDetail, int outFreq, double initAnnoSS); +RcppExport SEXP _SBayesRC_sbayesr_eigen_joint_annot(SEXP niterSEXP, SEXP burnSEXP, SEXP bhatSEXP, SEXP numAnnoSEXP, SEXP annoStrsSEXP, SEXP mldmDirSEXP, SEXP varySEXP, SEXP blkNSEXP, SEXP cgammaSEXP, SEXP startPiSEXP, SEXP starth2SEXP, SEXP cutThreshSEXP, SEXP bOriginSEXP, SEXP outPrefixSEXP, SEXP samVeSEXP, SEXP resam_threshSEXP, SEXP bOutDetailSEXP, SEXP outFreqSEXP, SEXP initAnnoSSSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -64,7 +64,9 @@ BEGIN_RCPP Rcpp::traits::input_parameter< std::string >::type samVe(samVeSEXP); Rcpp::traits::input_parameter< double >::type resam_thresh(resam_threshSEXP); Rcpp::traits::input_parameter< bool >::type bOutDetail(bOutDetailSEXP); - rcpp_result_gen = Rcpp::wrap(sbayesr_eigen_joint_annot(niter, burn, bhat, numAnno, annoStrs, mldmDir, vary, blkN, cgamma, startPi, starth2, cutThresh, bOrigin, outPrefix, samVe, resam_thresh, bOutDetail)); + Rcpp::traits::input_parameter< int >::type outFreq(outFreqSEXP); + Rcpp::traits::input_parameter< double >::type initAnnoSS(initAnnoSSSEXP); + rcpp_result_gen = Rcpp::wrap(sbayesr_eigen_joint_annot(niter, burn, bhat, numAnno, annoStrs, mldmDir, vary, blkN, cgamma, startPi, starth2, cutThresh, bOrigin, outPrefix, samVe, resam_thresh, bOutDetail, outFreq, initAnnoSS)); return rcpp_result_gen; END_RCPP } @@ -72,7 +74,7 @@ END_RCPP static const R_CallMethodDef CallEntries[] = { {"_SBayesRC_cutLD", (DL_FUNC) &_SBayesRC_cutLD, 3}, {"_SBayesRC_impGa", (DL_FUNC) &_SBayesRC_impGa, 6}, - {"_SBayesRC_sbayesr_eigen_joint_annot", (DL_FUNC) &_SBayesRC_sbayesr_eigen_joint_annot, 17}, + {"_SBayesRC_sbayesr_eigen_joint_annot", (DL_FUNC) &_SBayesRC_sbayesr_eigen_joint_annot, 19}, {NULL, NULL, 0} }; diff --git a/src/SBayesRC.cpp b/src/SBayesRC.cpp index 2f92387..b6ed051 100644 --- a/src/SBayesRC.cpp +++ b/src/SBayesRC.cpp @@ -22,7 +22,7 @@ #include "Timer.hpp" -SBayesRC::SBayesRC(int niter, int burn, VectorXf fbhat, int numAnno, vector &annoStrs, std::string mldmDir, double vary, VectorXf n, VectorXf fgamma, VectorXf pi, double starth2, double cutThresh, bool bOrigin, std::string outPrefix, std::string samVe, double resam_thresh, bool bOutDetail){ +SBayesRC::SBayesRC(int niter, int burn, VectorXf fbhat, int numAnno, vector &annoStrs, std::string mldmDir, double vary, VectorXf n, VectorXf fgamma, VectorXf pi, double starth2, double cutThresh, bool bOrigin, std::string outPrefix, std::string samVe, double resam_thresh, bool bOutDetail, int outFreq, double initAnnoSS){ this->niter = niter; this->burn = burn; @@ -35,6 +35,7 @@ SBayesRC::SBayesRC(int niter, int burn, VectorXf fbhat, int numAnno, vectoroutPrefix = outPrefix; this->curSamVe = samVe; this->fgamma = fgamma; + this->outFreq = outFreq; int nMarker = fbhat.size(); ndist = fgamma.size(); // number of mixture distribution components @@ -52,7 +53,7 @@ SBayesRC::SBayesRC(int niter, int burn, VectorXf fbhat, int numAnno, vectoropen(outPrefix, annoStrs); } //annoMat.resize(0, 0); @@ -560,6 +561,7 @@ void SBayesRC::mcmc(){ //NumericVector mean_par = as(wrap((sum(iter_infos, 0) / iter_infos.n_rows).as_col())); // clean if(anno){ + annoSS = anno->getSigmaSq(); if(!outPrefix.empty())anno->close(); delete anno; } @@ -615,3 +617,6 @@ MatrixXd SBayesRC::get_hsq_infos(){ return hsq_infos; } +VectorXd SBayesRC::get_anno_ss(){ + return annoSS; +} diff --git a/src/SBayesRC.hpp b/src/SBayesRC.hpp index 95fc5db..996d0aa 100644 --- a/src/SBayesRC.hpp +++ b/src/SBayesRC.hpp @@ -24,7 +24,7 @@ using Eigen::MatrixXf; class SBayesRC{ public: - SBayesRC(int niter, int burn, VectorXf fbhat, int numAnno, vector &annoStrs, std::string mldmDir, double vary, VectorXf n, VectorXf fgamma, VectorXf pi, double starth2=0.01, double cutThresh=1, bool bOrigin = false, std::string outPrefix="", std::string samVe = "fixVe", double resam_thresh=1.1, bool bOutDetail=false); + SBayesRC(int niter, int burn, VectorXf fbhat, int numAnno, vector &annoStrs, std::string mldmDir, double vary, VectorXf n, VectorXf fgamma, VectorXf pi, double starth2=0.01, double cutThresh=1, bool bOrigin = false, std::string outPrefix="", std::string samVe = "fixVe", double resam_thresh=1.1, bool bOutDetail=false, int outFreq=10, double initAnnoSS=1.0); void mcmc(); VectorXd get_mean_par_vec(); VectorXf get_betaMean_vec(); @@ -38,6 +38,8 @@ class SBayesRC{ MatrixXd get_vare_infos(); MatrixXd get_hsq_infos(); MatrixXd get_ssq_infos(); + VectorXd get_anno_ss(); + void setOutFreq(); private: int ndist; // number of mixture distribution components @@ -57,6 +59,7 @@ class SBayesRC{ double resam_thresh; VectorXf n; AnnoProb * anno = NULL; // annoation + VectorXd annoSS; BlockLDeig blockLDeig; // block LD eigen diff --git a/src/dist.cpp b/src/dist.cpp index 3f615bb..2d7fadd 100644 --- a/src/dist.cpp +++ b/src/dist.cpp @@ -30,7 +30,7 @@ namespace InvChiSq{ typedef boost::mt19937 random_engine; typedef boost::gamma_distribution<> gamma_distribution; typedef boost::variate_generator gamma_generator; - static random_engine engine; + static thread_local random_engine engine; float sample(const float df, const float scale){ gamma_generator sgamma(engine, gamma_distribution(0.5f*df, 1)); return scale/(2.0f*sgamma()); @@ -44,8 +44,8 @@ namespace Bernoulli{ //static random_engine engine; //static uniform01_generator ranf(engine, uniform_01()); - static std::mt19937_64 rng; - static std::uniform_real_distribution unif(0, 1); + static thread_local std::mt19937_64 rng; + static thread_local std::uniform_real_distribution unif(0, 1); int sample(const VectorXf &p){ float cum = 0; @@ -68,7 +68,7 @@ namespace Gamma{ typedef boost::mt19937 random_engine; typedef boost::gamma_distribution<> gamma_distribution; typedef boost::variate_generator gamma_generator; - static random_engine engine; + static thread_local random_engine engine; float sample(const float shape, const float scale){ gamma_generator sgamma(engine, gamma_distribution(shape, scale)); return sgamma(); @@ -90,13 +90,13 @@ namespace Dirichlet{ namespace Normal{ typedef boost::mt19937 random_engine; - static random_engine engine; + static thread_local random_engine engine; typedef boost::uniform_01<> uniform_01; typedef boost::normal_distribution<> normal_distribution; typedef boost::variate_generator uniform01_generator; typedef boost::variate_generator normal_generator; - static uniform01_generator ranf(engine, uniform_01()); - static normal_generator snorm(engine, normal_distribution(0,1)); // standard normal + static thread_local uniform01_generator ranf(engine, uniform_01()); + static thread_local normal_generator snorm(engine, normal_distribution(0,1)); // standard normal boost::math::normal_distribution <> d = boost::math::normal_distribution <> (0 ,1); double quantile_01(double value){ diff --git a/src/sbr_eigen_c.cpp b/src/sbr_eigen_c.cpp index 699c179..2ef3be3 100644 --- a/src/sbr_eigen_c.cpp +++ b/src/sbr_eigen_c.cpp @@ -30,7 +30,7 @@ using Eigen::MatrixXd; // [[Rcpp::export]] -List sbayesr_eigen_joint_annot(int niter, int burn, Eigen::Map bhat, int numAnno, Rcpp::StringVector annoStrs, std::string mldmDir, double vary, Eigen::Map blkN, Eigen::Map cgamma, Eigen::Map startPi, double starth2=0.01, double cutThresh=1, bool bOrigin = false, std::string outPrefix="", std::string samVe = "fixVe", double resam_thresh=1.1, bool bOutDetail=false){ +List sbayesr_eigen_joint_annot(int niter, int burn, Eigen::Map bhat, int numAnno, Rcpp::StringVector annoStrs, std::string mldmDir, double vary, Eigen::Map blkN, Eigen::Map cgamma, Eigen::Map startPi, double starth2=0.01, double cutThresh=1, bool bOrigin = false, std::string outPrefix="", std::string samVe = "fixVe", double resam_thresh=1.1, bool bOutDetail=false, int outFreq=10, double initAnnoSS=1.0){ string curSamVe = samVe; VectorXf fbhat = bhat.cast(); @@ -43,7 +43,7 @@ List sbayesr_eigen_joint_annot(int niter, int burn, Eigen::Map annoStrings[i] = annoStrs[i]; } - SBayesRC sbr(niter, burn, fbhat, numAnno, annoStrings, mldmDir, vary, n, fgamma, pi, starth2, cutThresh, bOrigin, outPrefix, samVe, resam_thresh, bOutDetail); + SBayesRC sbr(niter, burn, fbhat, numAnno, annoStrings, mldmDir, vary, n, fgamma, pi, starth2, cutThresh, bOrigin, outPrefix, samVe, resam_thresh, bOutDetail, outFreq, initAnnoSS); sbr.mcmc(); // return @@ -74,6 +74,8 @@ List sbayesr_eigen_joint_annot(int niter, int burn, Eigen::Map MatrixXf pip = sbr.get_pip(); NumericMatrix pip_r(pip.rows(), pip.cols(), pip.data()); //NumericVector weightQ_r(weightQ.data(), weightQ.data() + weightQ.size()); + Eigen::VectorXd sigmaAnno = sbr.get_anno_ss(); + NumericVector sigmaAnno_r(sigmaAnno.data(), sigmaAnno.data() + sigmaAnno.size()); return(List::create(_["par"]=mean_par, _["betaMean"] = betaMean, @@ -86,6 +88,7 @@ List sbayesr_eigen_joint_annot(int niter, int burn, Eigen::Map _["hsq_block_hist"] = sbr.get_hsq_infos(), _["pip"] = pip_r, _["vg_comp_hist"] = vg_comp_mcmc_r, - _["n_comp_hist"] = n_comp_mcmc_r + _["n_comp_hist"] = n_comp_mcmc_r, + _["sigma_anno"] = sigmaAnno_r )); }