add the outfreq, add sigmaSq, bug fix

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"),

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)
 }
 

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)
     }
 

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<double>();
+}
+
+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);

src/AnnoProb.hpp

@@ -70,7 +70,7 @@ class AnnoProb {
     void writeHeader(const vector<string> &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<string> &annoStrs);
     void close();
     void output();
+    VectorXd getSigmaSq();
 };
 
 #endif //SBRC_ANNO_PROB_H

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<Eigen::VectorXd> bhat, int numAnno, Rcpp::StringVector annoStrs, std::string mldmDir, double vary, Eigen::Map<Eigen::VectorXd> blkN, Eigen::Map<Eigen::VectorXd> cgamma, Eigen::Map<Eigen::VectorXd> 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<Eigen::VectorXd> bhat, int numAnno, Rcpp::StringVector annoStrs, std::string mldmDir, double vary, Eigen::Map<Eigen::VectorXd> blkN, Eigen::Map<Eigen::VectorXd> cgamma, Eigen::Map<Eigen::VectorXd> 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,15 +64,17 @@ 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
 }
 
 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}
 };
 

src/SBayesRC.cpp

@@ -22,7 +22,7 @@
 #include "Timer.hpp"
 
 
-SBayesRC::SBayesRC(int niter, int burn, VectorXf fbhat, int numAnno, vector<string> &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<string> &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, vector<stri
     this->outPrefix = 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, vector<stri
 
     if(bAnnot){
         snpPi.resize(nMarker, ndist);
-        anno = new AnnoProb(fileAnnot, numAnno, pi, snpPi, bOutDetail);
+        anno = new AnnoProb(fileAnnot, numAnno, pi, snpPi, bOutDetail, initAnnoSS);
         if(!outPrefix.empty()) anno->open(outPrefix, annoStrs);
     }
     //annoMat.resize(0, 0);
@@ -560,6 +561,7 @@ void SBayesRC::mcmc(){
     //NumericVector mean_par = as<NumericVector>(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;
+}

src/SBayesRC.hpp

@@ -24,7 +24,7 @@ using Eigen::MatrixXf;
 class SBayesRC{
 
 public:
-    SBayesRC(int niter, int burn, VectorXf fbhat, int numAnno, vector<string> &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<string> &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
 
 

src/dist.cpp

@@ -30,7 +30,7 @@ namespace InvChiSq{
     typedef boost::mt19937 random_engine;
     typedef boost::gamma_distribution<> gamma_distribution;
     typedef boost::variate_generator<random_engine&, gamma_distribution> 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<float> unif(0, 1);
+    static thread_local std::mt19937_64 rng;
+    static thread_local std::uniform_real_distribution<float> 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<random_engine&, gamma_distribution> 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<random_engine&, uniform_01> uniform01_generator;
     typedef boost::variate_generator<random_engine&, normal_distribution> 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){

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<Eigen::VectorXd> bhat, int numAnno, Rcpp::StringVector annoStrs, std::string mldmDir, double vary, Eigen::Map<Eigen::VectorXd> blkN, Eigen::Map<Eigen::VectorXd> cgamma, Eigen::Map<Eigen::VectorXd> 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<Eigen::VectorXd> bhat, int numAnno, Rcpp::StringVector annoStrs, std::string mldmDir, double vary, Eigen::Map<Eigen::VectorXd> blkN, Eigen::Map<Eigen::VectorXd> cgamma, Eigen::Map<Eigen::VectorXd> 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<float>();
@@ -43,7 +43,7 @@ List sbayesr_eigen_joint_annot(int niter, int burn, Eigen::Map<Eigen::VectorXd>
         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<Eigen::VectorXd>
     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<Eigen::VectorXd>
                 _["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
                 )); 
 }