diff --git a/DESCRIPTION b/DESCRIPTION index c7bb580b..b3d55171 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Package: BASiCS Type: Package Title: Bayesian Analysis of Single-Cell Sequencing data Version: 2.1.8 -Date: 2020-09-28 +Date: 2020-09-29 Authors@R: c(person("Catalina", "Vallejos", role=c("aut", "cre"), email="catalina.vallejos@igmm.ed.ac.uk"), person("Nils", "Eling", role=c("aut")), diff --git a/NEWS b/NEWS index 7ec00189..0d5e0cff 100644 --- a/NEWS +++ b/NEWS @@ -1519,6 +1519,14 @@ version 2.1.7 (2020-09-23) Also use a slightly more principled calculation internally. - See issues #182, #39, #91, #169, #181, #178, #202, #190, #173 + version 2.1.8 (2020-09-28) - Minor edit in the vignette to recommend the use of a seed for reproducible -results when using `BASiCS_MCMC`. \ No newline at end of file +results when using `BASiCS_MCMC`. + + +version 2.1.9 (2020-09-28) +- Add `Threads` argument to `BASiCS_MCMC`, allowing parallelisation of MCMC + updates across cells or genes using openMP. +- Reduce BatchInfo requirement for no spikes sampler to a warning, from an + error. diff --git a/R/BASiCS_MCMC.R b/R/BASiCS_MCMC.R index 74cdfa65..1312ff4a 100644 --- a/R/BASiCS_MCMC.R +++ b/R/BASiCS_MCMC.R @@ -196,10 +196,11 @@ BASiCS_MCMC <- function( Burn, Regression, WithSpikes = TRUE, + Threads = getOption("mc.cores", default = 1L), ...) { # Checks to ensure input arguments are valid - .BASiCS_MCMC_InputCheck(Data, N, Thin, Burn, Regression, WithSpikes) + .BASiCS_MCMC_InputCheck(Data, N, Thin, Burn, Regression, WithSpikes, Threads) # Some global values used throughout the MCMC algorithm and checks # Numbers of cells/genes/batches + count table + design matrix for batch @@ -265,7 +266,8 @@ BASiCS_MCMC <- function( mintol_nu = ArgsDef$mintol_nu, mintol_theta = ArgsDef$mintol_theta, geneExponent = PriorParam$GeneExponent, - cellExponent = PriorParam$CellExponent + cellExponent = PriorParam$CellExponent, + threads = Threads ) diff --git a/R/BASiCS_Package.R b/R/BASiCS_Package.R index 5877f6f6..22482f40 100644 --- a/R/BASiCS_Package.R +++ b/R/BASiCS_Package.R @@ -15,7 +15,7 @@ #' @importFrom matrixStats colMeans2 colMedians colVars #' @importFrom matrixStats rowMeans2 rowMedians rowVars #' @importFrom methods .hasSlot is new show slotNames Summary -#' @importFrom Rcpp evalCpp +#' @importFrom Rcpp evalCpp sourceCpp #' @importFrom S4Vectors DataFrame #' @importFrom scran calculateSumFactors #' @importFrom stats acf median model.matrix rgamma rpois runif var diff --git a/R/RcppExports.R b/R/RcppExports.R index 63b6964d..23729456 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,32 +1,24 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 -.BASiCS_MCMCcpp <- function(N, Thin, Burn, Counts, BatchDesign, muSpikes, mu0, delta0, phi0, s0, nu0, theta0, mu_mu, s2mu, adelta, bdelta, s2delta, prior_delta, aphi, as, bs, atheta, btheta, ar, LSmu0, LSdelta0, LSphi0, LSnu0, LStheta0, sumByCellBio, sumByGeneAll, sumByGeneBio, StoreAdapt, EndAdapt, PrintProgress, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent) { - .Call('_BASiCS_BASiCS_MCMCcpp', PACKAGE = 'BASiCS', N, Thin, Burn, Counts, BatchDesign, muSpikes, mu0, delta0, phi0, s0, nu0, theta0, mu_mu, s2mu, adelta, bdelta, s2delta, prior_delta, aphi, as, bs, atheta, btheta, ar, LSmu0, LSdelta0, LSphi0, LSnu0, LStheta0, sumByCellBio, sumByGeneAll, sumByGeneBio, StoreAdapt, EndAdapt, PrintProgress, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent) -} - -.BASiCS_MCMCcppReg <- function(N, Thin, Burn, Counts, BatchDesign, muSpikes, mu0, delta0, phi0, s0, nu0, theta0, mu_mu, s2mu, aphi, as, bs, atheta, btheta, k, m0, V0, sigma2_a0, sigma2_b0, beta0, sigma20, eta0, lambda0, variance, ar, LSmu0, LSdelta0, LSphi0, LSnu0, LStheta0, sumByCellBio, sumByGeneAll, sumByGeneBio, StoreAdapt, EndAdapt, PrintProgress, FixLocations, RBFMinMax, RBFLocations, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent) { - .Call('_BASiCS_BASiCS_MCMCcppReg', PACKAGE = 'BASiCS', N, Thin, Burn, Counts, BatchDesign, muSpikes, mu0, delta0, phi0, s0, nu0, theta0, mu_mu, s2mu, aphi, as, bs, atheta, btheta, k, m0, V0, sigma2_a0, sigma2_b0, beta0, sigma20, eta0, lambda0, variance, ar, LSmu0, LSdelta0, LSphi0, LSnu0, LStheta0, sumByCellBio, sumByGeneAll, sumByGeneBio, StoreAdapt, EndAdapt, PrintProgress, FixLocations, RBFMinMax, RBFLocations, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent) -} - -.BASiCS_MCMCcppNoSpikes <- function(N, Thin, Burn, Counts, BatchDesign, mu0, delta0, s0, nu0, theta0, mu_mu, s2mu, adelta, bdelta, s2delta, prior_delta, as, bs, atheta, btheta, Constrain, Index, RefGene, RefGenes, ConstrainGene, NotConstrainGene, StochasticRef, ar, LSmu0, LSdelta0, LSnu0, LStheta0, sumByCellAll, sumByGeneAll, StoreAdapt, EndAdapt, PrintProgress, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent) { - .Call('_BASiCS_BASiCS_MCMCcppNoSpikes', PACKAGE = 'BASiCS', N, Thin, Burn, Counts, BatchDesign, mu0, delta0, s0, nu0, theta0, mu_mu, s2mu, adelta, bdelta, s2delta, prior_delta, as, bs, atheta, btheta, Constrain, Index, RefGene, RefGenes, ConstrainGene, NotConstrainGene, StochasticRef, ar, LSmu0, LSdelta0, LSnu0, LStheta0, sumByCellAll, sumByGeneAll, StoreAdapt, EndAdapt, PrintProgress, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent) +.BASiCS_DenoisedRates <- function(CountsBio, Mu, TransInvDelta, PhiNu, N, q0, n) { + .Call('_BASiCS_BASiCS_DenoisedRates', PACKAGE = 'BASiCS', CountsBio, Mu, TransInvDelta, PhiNu, N, q0, n) } -.BASiCS_MCMCcppRegNoSpikes <- function(N, Thin, Burn, Counts, BatchDesign, mu0, delta0, s0, nu0, theta0, mu_mu, s2mu, as, bs, atheta, btheta, k, m0, V0, sigma2_a0, sigma2_b0, beta0, sigma20, eta0, lambda0, variance, Constrain, Index, RefGene, RefGenes, ConstrainGene, NotConstrainGene, StochasticRef, ar, LSmu0, LSdelta0, LSnu0, LStheta0, sumByCellAll, sumByGeneAll, StoreAdapt, EndAdapt, PrintProgress, RBFMinMax, FixLocations, RBFLocations, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent) { - .Call('_BASiCS_BASiCS_MCMCcppRegNoSpikes', PACKAGE = 'BASiCS', N, Thin, Burn, Counts, BatchDesign, mu0, delta0, s0, nu0, theta0, mu_mu, s2mu, as, bs, atheta, btheta, k, m0, V0, sigma2_a0, sigma2_b0, beta0, sigma20, eta0, lambda0, variance, Constrain, Index, RefGene, RefGenes, ConstrainGene, NotConstrainGene, StochasticRef, ar, LSmu0, LSdelta0, LSnu0, LStheta0, sumByCellAll, sumByGeneAll, StoreAdapt, EndAdapt, PrintProgress, RBFMinMax, FixLocations, RBFLocations, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent) +.BASiCS_MCMCcpp <- function(N, Thin, Burn, Counts, BatchDesign, muSpikes, mu0, delta0, phi0, s0, nu0, theta0, mu_mu, s2mu, adelta, bdelta, s2delta, prior_delta, aphi, as, bs, atheta, btheta, ar, LSmu0, LSdelta0, LSphi0, LSnu0, LStheta0, sumByCellBio, sumByGeneAll, sumByGeneBio, StoreAdapt, EndAdapt, PrintProgress, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent, threads = 1L) { + .Call('_BASiCS_BASiCS_MCMCcpp', PACKAGE = 'BASiCS', N, Thin, Burn, Counts, BatchDesign, muSpikes, mu0, delta0, phi0, s0, nu0, theta0, mu_mu, s2mu, adelta, bdelta, s2delta, prior_delta, aphi, as, bs, atheta, btheta, ar, LSmu0, LSdelta0, LSphi0, LSnu0, LStheta0, sumByCellBio, sumByGeneAll, sumByGeneBio, StoreAdapt, EndAdapt, PrintProgress, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent, threads) } -.BASiCS_DenoisedRates <- function(CountsBio, Mu, TransInvDelta, PhiNu, N, q0, n) { - .Call('_BASiCS_BASiCS_DenoisedRates', PACKAGE = 'BASiCS', CountsBio, Mu, TransInvDelta, PhiNu, N, q0, n) +.BASiCS_MCMCcppNoSpikes <- function(N, Thin, Burn, Counts, BatchDesign, mu0, delta0, s0, nu0, theta0, mu_mu, s2mu, adelta, bdelta, s2delta, prior_delta, as, bs, atheta, btheta, Constrain, Index, RefGene, RefGenes, ConstrainGene, NotConstrainGene, StochasticRef, ar, LSmu0, LSdelta0, LSnu0, LStheta0, sumByCellAll, sumByGeneAll, StoreAdapt, EndAdapt, PrintProgress, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent, threads = 1L) { + .Call('_BASiCS_BASiCS_MCMCcppNoSpikes', PACKAGE = 'BASiCS', N, Thin, Burn, Counts, BatchDesign, mu0, delta0, s0, nu0, theta0, mu_mu, s2mu, adelta, bdelta, s2delta, prior_delta, as, bs, atheta, btheta, Constrain, Index, RefGene, RefGenes, ConstrainGene, NotConstrainGene, StochasticRef, ar, LSmu0, LSdelta0, LSnu0, LStheta0, sumByCellAll, sumByGeneAll, StoreAdapt, EndAdapt, PrintProgress, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent, threads) } -.estimateRBFLocations <- function(log_mu, k, RBFMinMax) { - .Call('_BASiCS_estimateRBFLocations', PACKAGE = 'BASiCS', log_mu, k, RBFMinMax) +.BASiCS_MCMCcppReg <- function(N, Thin, Burn, Counts, BatchDesign, muSpikes, mu0, delta0, phi0, s0, nu0, theta0, mu_mu, s2mu, aphi, as, bs, atheta, btheta, k, m0, V0, sigma2_a0, sigma2_b0, beta0, sigma20, eta0, lambda0, variance, ar, LSmu0, LSdelta0, LSphi0, LSnu0, LStheta0, sumByCellBio, sumByGeneAll, sumByGeneBio, StoreAdapt, EndAdapt, PrintProgress, FixLocations, RBFMinMax, RBFLocations, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent, threads = 1L) { + .Call('_BASiCS_BASiCS_MCMCcppReg', PACKAGE = 'BASiCS', N, Thin, Burn, Counts, BatchDesign, muSpikes, mu0, delta0, phi0, s0, nu0, theta0, mu_mu, s2mu, aphi, as, bs, atheta, btheta, k, m0, V0, sigma2_a0, sigma2_b0, beta0, sigma20, eta0, lambda0, variance, ar, LSmu0, LSdelta0, LSphi0, LSnu0, LStheta0, sumByCellBio, sumByGeneAll, sumByGeneBio, StoreAdapt, EndAdapt, PrintProgress, FixLocations, RBFMinMax, RBFLocations, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent, threads) } -.rDirichlet <- function(alpha) { - .Call('_BASiCS_rDirichlet', PACKAGE = 'BASiCS', alpha) +.BASiCS_MCMCcppRegNoSpikes <- function(N, Thin, Burn, Counts, BatchDesign, mu0, delta0, s0, nu0, theta0, mu_mu, s2mu, as, bs, atheta, btheta, k, m0, V0, sigma2_a0, sigma2_b0, beta0, sigma20, eta0, lambda0, variance, Constrain, Index, RefGene, RefGenes, ConstrainGene, NotConstrainGene, StochasticRef, ar, LSmu0, LSdelta0, LSnu0, LStheta0, sumByCellAll, sumByGeneAll, StoreAdapt, EndAdapt, PrintProgress, RBFMinMax, FixLocations, RBFLocations, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent, threads = 1L) { + .Call('_BASiCS_BASiCS_MCMCcppRegNoSpikes', PACKAGE = 'BASiCS', N, Thin, Burn, Counts, BatchDesign, mu0, delta0, s0, nu0, theta0, mu_mu, s2mu, as, bs, atheta, btheta, k, m0, V0, sigma2_a0, sigma2_b0, beta0, sigma20, eta0, lambda0, variance, Constrain, Index, RefGene, RefGenes, ConstrainGene, NotConstrainGene, StochasticRef, ar, LSmu0, LSdelta0, LSnu0, LStheta0, sumByCellAll, sumByGeneAll, StoreAdapt, EndAdapt, PrintProgress, RBFMinMax, FixLocations, RBFLocations, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent, threads) } .muUpdate <- function(mu0, prop_var, Counts, invdelta, phinu, sum_bycell_bio, mu_mu, s2_mu, q0, n, mu1, u, ind, exponent, mintol) { @@ -53,10 +45,22 @@ .Call('_BASiCS_thetaUpdateBatch', PACKAGE = 'BASiCS', theta0, prop_var, BatchDesign, BatchSizes, s, nu, a_theta, b_theta, n, nBatch, exponent, mintol) } +.muUpdateNoSpikes <- function(mu0, prop_var, Counts, invdelta, nu, sum_bycell_all, mu_mu, s2_mu, q0, n, mu1, u, ind, Constrain, RefGene, ConstrainGene, NotConstrainGene, exponent, mintol) { + .Call('_BASiCS_muUpdateNoSpikes', PACKAGE = 'BASiCS', mu0, prop_var, Counts, invdelta, nu, sum_bycell_all, mu_mu, s2_mu, q0, n, mu1, u, ind, Constrain, RefGene, ConstrainGene, NotConstrainGene, exponent, mintol) +} + +.nuUpdateBatchNoSpikes <- function(nu0, prop_var, Counts, BatchDesign, mu, invdelta, s, thetaBatch, sum_bygene_all, q0, n, nu1, u, ind, exponent, mintol) { + .Call('_BASiCS_nuUpdateBatchNoSpikes', PACKAGE = 'BASiCS', nu0, prop_var, Counts, BatchDesign, mu, invdelta, s, thetaBatch, sum_bygene_all, q0, n, nu1, u, ind, exponent, mintol) +} + .designMatrix <- function(k, RBFLocations, mu, variance) { .Call('_BASiCS_designMatrix', PACKAGE = 'BASiCS', k, RBFLocations, mu, variance) } +.estimateRBFLocations <- function(log_mu, k, RBFMinMax) { + .Call('_BASiCS_estimateRBFLocations', PACKAGE = 'BASiCS', log_mu, k, RBFMinMax) +} + .muUpdateReg <- function(mu0, prop_var, Counts, delta, phinu, sum_bycell_bio, mu_mu, s2_mu, q0, n, mu1, u, ind, k, lambda, beta, X, sigma2, variance, FixLocations, RBFMinMax, RBFLocations, exponent, mintol) { .Call('_BASiCS_muUpdateReg', PACKAGE = 'BASiCS', mu0, prop_var, Counts, delta, phinu, sum_bycell_bio, mu_mu, s2_mu, q0, n, mu1, u, ind, k, lambda, beta, X, sigma2, variance, FixLocations, RBFMinMax, RBFLocations, exponent, mintol) } @@ -77,14 +81,6 @@ .Call('_BASiCS_lambdaUpdateReg', PACKAGE = 'BASiCS', delta, X, beta, sigma2, eta, q0, lambda1, exponent) } -.muUpdateNoSpikes <- function(mu0, prop_var, Counts, invdelta, nu, sum_bycell_all, mu_mu, s2_mu, q0, n, mu1, u, ind, Constrain, RefGene, ConstrainGene, NotConstrainGene, exponent, mintol) { - .Call('_BASiCS_muUpdateNoSpikes', PACKAGE = 'BASiCS', mu0, prop_var, Counts, invdelta, nu, sum_bycell_all, mu_mu, s2_mu, q0, n, mu1, u, ind, Constrain, RefGene, ConstrainGene, NotConstrainGene, exponent, mintol) -} - -.nuUpdateBatchNoSpikes <- function(nu0, prop_var, Counts, BatchDesign, mu, invdelta, s, thetaBatch, sum_bygene_all, q0, n, nu1, u, ind, exponent, mintol) { - .Call('_BASiCS_nuUpdateBatchNoSpikes', PACKAGE = 'BASiCS', nu0, prop_var, Counts, BatchDesign, mu, invdelta, s, thetaBatch, sum_bygene_all, q0, n, nu1, u, ind, exponent, mintol) -} - .muUpdateRegNoSpikes <- function(mu0, prop_var, Counts, delta, invdelta, nu, sum_bycell_all, mu_mu, s2_mu, q0, n, mu1, u, ind, Constrain, RefGene, ConstrainGene, NotConstrainGene, k, lambda, beta, X, sigma2, variance, FixLocations, RBFMinMax, RBFLocations, exponent, mintol) { .Call('_BASiCS_muUpdateRegNoSpikes', PACKAGE = 'BASiCS', mu0, prop_var, Counts, delta, invdelta, nu, sum_bycell_all, mu_mu, s2_mu, q0, n, mu1, u, ind, Constrain, RefGene, ConstrainGene, NotConstrainGene, k, lambda, beta, X, sigma2, variance, FixLocations, RBFMinMax, RBFLocations, exponent, mintol) } @@ -93,3 +89,7 @@ .Call('_BASiCS_deltaUpdateRegNoSpikes', PACKAGE = 'BASiCS', delta0, prop_var, Counts, mu, nu, q0, n, delta1, u, ind, lambda, X, sigma2, beta, exponent, mintol) } +.rDirichlet <- function(alpha) { + .Call('_BASiCS_rDirichlet', PACKAGE = 'BASiCS', alpha) +} + diff --git a/R/utils_MCMC.R b/R/utils_MCMC.R index 82316ef6..29e2fe65 100644 --- a/R/utils_MCMC.R +++ b/R/utils_MCMC.R @@ -3,12 +3,19 @@ Thin, Burn, Regression, - WithSpikes) { + WithSpikes, + Threads) { if (!is(Data, "SingleCellExperiment")) { stop("'Data' is not a SingleCellExperiment class object.\n") } - + if (!is.numeric(Threads) && + length(Threads) == 1 && + round(Threads) == Threads && + Threads > 0) { + stop("Threads must be a positive integer.") + } + # Check if `Data` contains more than one type of types if ((length(SingleCellExperiment::altExpNames(Data)) > 1)) { stop( @@ -24,7 +31,7 @@ "see help(altExp) for details. \n" ) # If SpikeInput slot is missing and WithSpikes == TRUE - if(is.null(rowData(altExp(Data)))) { + if (is.null(rowData(altExp(Data)))) { stop( "'altExp(Data)' does not contain 'rowData' \n" ) @@ -41,17 +48,17 @@ } # If isSpike slot is missing and WithSpikes == TRUE - if((length(altExpNames(Data)) == 0) & WithSpikes) { + if ((length(altExpNames(Data)) == 0) & WithSpikes) { stop( "'Data' does not contain information about spike-in genes \n", - "Please indicate include this information using 'altExp' \n", + "Please include this information using 'altExp' \n", "or set 'WithSpikes = FALSE' \n." ) } # If BatchInfo slot is missing and WithSpikes == FALSE - if(!WithSpikes & is.null(colData(Data)$BatchInfo)) { - stop( + if (!WithSpikes & is.null(colData(Data)$BatchInfo)) { + warning( "'Data' should contain a BatchInfo vector when 'WithSpikes = FALSE'. \n", "Please assign the batch information to: \n 'colData(Data)$BatchInfo = BatchInfo'. \n" @@ -59,7 +66,7 @@ } # Checking how counts are stored - if(!("counts" %in% assayNames(Data))) + if (!("counts" %in% assayNames(Data))) stop( "'Data' does not contain a 'counts' slot. \n", "Please make sure to include the raw data in the \n", diff --git a/src/BASiCS.cpp b/src/BASiCS.cpp new file mode 100644 index 00000000..52a827f1 --- /dev/null +++ b/src/BASiCS.cpp @@ -0,0 +1,9 @@ +#include "updates.h" +#include "updatesReg.h" +#include "updatesNoSpikes.h" +#include "updatesRegNoSpikes.h" +#include "MCMC.h" +#include "MCMCReg.h" +#include "MCMCNoSpikes.h" +#include "MCMCRegNoSpikes.h" +#include "DenoisedRates.h" diff --git a/src/BASiCS_DenoisedRates.cpp b/src/DenoisedRates.h similarity index 90% rename from src/BASiCS_DenoisedRates.cpp rename to src/DenoisedRates.h index df466055..a3d2b9af 100644 --- a/src/BASiCS_DenoisedRates.cpp +++ b/src/DenoisedRates.h @@ -1,5 +1,9 @@ +#ifndef DENOISED_H +#define DENOISED_H + #include "utils.h" +// [[Rcpp::export(".BASiCS_DenoisedRates")]] arma::mat BASiCS_DenoisedRates( NumericMatrix CountsBio, NumericMatrix Mu, @@ -33,8 +37,4 @@ arma::mat BASiCS_DenoisedRates( return(Rho / N); } - - - - - +#endif diff --git a/src/BASiCS_MCMCcpp.cpp b/src/MCMC.h similarity index 98% rename from src/BASiCS_MCMCcpp.cpp rename to src/MCMC.h index 36aca19b..7dc9f1d4 100644 --- a/src/BASiCS_MCMCcpp.cpp +++ b/src/MCMC.h @@ -1,3 +1,6 @@ +#ifndef MCMC_H +#define MCMC_H + #include "utils.h" /* MCMC sampler @@ -36,6 +39,7 @@ * EndAdapt: when to stop the adaptation * PrintProgress: whether to print progress report */ +// [[Rcpp::export(".BASiCS_MCMCcpp")]] Rcpp::List BASiCS_MCMCcpp( int N, int Thin, @@ -77,7 +81,12 @@ Rcpp::List BASiCS_MCMCcpp( double const& mintol_nu, double const& mintol_theta, double const& geneExponent, - double const& cellExponent) { + double const& cellExponent, + int threads = 1) { + + #if defined(_OPENMP) + omp_set_num_threads(threads); + #endif using arma::ones; using arma::zeros; @@ -406,4 +415,4 @@ Rcpp::List BASiCS_MCMCcpp( } } - +#endif diff --git a/src/BASiCS_MCMCcppNoSpikes.cpp b/src/MCMCNoSpikes.h similarity index 98% rename from src/BASiCS_MCMCcppNoSpikes.cpp rename to src/MCMCNoSpikes.h index 54801b58..945091f7 100644 --- a/src/BASiCS_MCMCcppNoSpikes.cpp +++ b/src/MCMCNoSpikes.h @@ -1,3 +1,6 @@ +#ifndef MCMCNOSPIKES_H +#define MCMCNOSPIKES_H + #include "utils.h" /* MCMC sampler for the non-spike case @@ -46,6 +49,7 @@ * NotConstrainGene: * StochasticRef: */ +// [[Rcpp::export(".BASiCS_MCMCcppNoSpikes")]] Rcpp::List BASiCS_MCMCcppNoSpikes( int N, int Thin, @@ -89,7 +93,13 @@ Rcpp::List BASiCS_MCMCcppNoSpikes( double const& mintol_nu, double const& mintol_theta, double const& geneExponent, - double const& cellExponent) { + double const& cellExponent, + int threads = 1) { + + + #if defined(_OPENMP) + omp_set_num_threads(threads); + #endif using arma::ones; using arma::zeros; @@ -400,3 +410,5 @@ Rcpp::List BASiCS_MCMCcppNoSpikes( Rcpp::Named("RefFreq") = RefFreq/(N-Burn))); } } + +#endif diff --git a/src/BASiCS_MCMCcppReg.cpp b/src/MCMCReg.h old mode 100755 new mode 100644 similarity index 98% rename from src/BASiCS_MCMCcppReg.cpp rename to src/MCMCReg.h index ccdc9966..d049c0a1 --- a/src/BASiCS_MCMCcppReg.cpp +++ b/src/MCMCReg.h @@ -1,3 +1,6 @@ +#ifndef MCMCREG_H +#define MCMCREG_H + #include "utils.h" /* MCMC sampler for regression case @@ -43,6 +46,7 @@ * lambda0: Starting values for gene-wise error term * variance: Fixed width (scale) for GRBFs */ +// [[Rcpp::export(".BASiCS_MCMCcppReg")]] Rcpp::List BASiCS_MCMCcppReg( int N, int Thin, @@ -93,7 +97,13 @@ Rcpp::List BASiCS_MCMCcppReg( double const& mintol_nu, double const& mintol_theta, double const& geneExponent, - double const& cellExponent) { + double const& cellExponent, + int threads = 1) { + + + #if defined(_OPENMP) + omp_set_num_threads(threads); + #endif using arma::ones; using arma::zeros; @@ -548,3 +558,5 @@ Rcpp::List BASiCS_MCMCcppReg( } } + +#endif diff --git a/src/BASiCS_MCMCcppRegNoSpikes.cpp b/src/MCMCRegNoSpikes.h old mode 100755 new mode 100644 similarity index 98% rename from src/BASiCS_MCMCcppRegNoSpikes.cpp rename to src/MCMCRegNoSpikes.h index 68393ad3..9e7e0ec6 --- a/src/BASiCS_MCMCcppRegNoSpikes.cpp +++ b/src/MCMCRegNoSpikes.h @@ -1,3 +1,6 @@ +#ifndef MCMCREGNOSPIKES_H +#define MCMCREGNOSPIKES_H + #include "utils.h" /* MCMC sampler for regression and non-spikes case @@ -43,6 +46,7 @@ * lambda0: Starting values for gene-wise error term * variance: Fixed width (scale) for GRBFs */ +// [[Rcpp::export(".BASiCS_MCMCcppRegNoSpikes")]] Rcpp::List BASiCS_MCMCcppRegNoSpikes( int N, int Thin, @@ -95,7 +99,14 @@ Rcpp::List BASiCS_MCMCcppRegNoSpikes( double const& mintol_nu, double const& mintol_theta, double const& geneExponent, - double const& cellExponent) { + double const& cellExponent, + int threads = 1) { + + + #if defined(_OPENMP) + omp_set_num_threads(threads); + // REprintf("Number of threads=%i\n", omp_get_max_threads()); + #endif using arma::ones; using arma::zeros; @@ -524,3 +535,5 @@ Rcpp::List BASiCS_MCMCcppRegNoSpikes( ); } } + +#endif diff --git a/src/MCMCcpp.h b/src/MCMCcpp.h deleted file mode 100644 index 3b6820ea..00000000 --- a/src/MCMCcpp.h +++ /dev/null @@ -1,233 +0,0 @@ -/* C++ implementation of BASiCS - * Author: Catalina A. Vallejos (cnvallej@uc.cl) & Nils Eling - */ - -#ifndef MCMCCPP_H -#define MCMCCPP_H - -#include -#include -#include -#include -#include -//File: matprod_arma.cpp -//[[Rcpp::depends(RcppArmadillo)]] -#include - -using namespace Rcpp; - -// Main MCMC sampler -// [[Rcpp::export(".BASiCS_MCMCcpp")]] -Rcpp::List BASiCS_MCMCcpp( - int N, - int Thin, - int Burn, - arma::mat Counts, - arma::mat BatchDesign, - arma::vec muSpikes, - arma::vec mu0, - arma::vec delta0, - arma::vec phi0, - arma::vec s0, - arma::vec nu0, - arma::vec theta0, - arma::vec mu_mu, - double s2mu, - double adelta, - double bdelta, - double s2delta, - double prior_delta, - arma::vec aphi, - double as, - double bs, - double atheta, - double btheta, - double ar, - arma::vec LSmu0, - arma::vec LSdelta0, - double LSphi0, - arma::vec LSnu0, - arma::vec LStheta0, - arma::vec sumByCellBio, - arma::vec sumByGeneAll, - arma::vec sumByGeneBio, - int StoreAdapt, - int EndAdapt, - int PrintProgress, - double const& mintol_mu, - double const& mintol_delta, - double const& mintol_nu, - double const& mintol_theta, - double const& geneExponent, - double const& cellExponent); - -// MCMC sampler for regression case -// [[Rcpp::export(".BASiCS_MCMCcppReg")]] -Rcpp::List BASiCS_MCMCcppReg( - int N, - int Thin, - int Burn, - arma::mat Counts, - arma::mat BatchDesign, - arma::vec muSpikes, - arma::vec mu0, - arma::vec delta0, - arma::vec phi0, - arma::vec s0, - arma::vec nu0, - arma::vec theta0, - arma::vec mu_mu, - double s2mu, - arma::vec aphi, - double as, - double bs, - double atheta, - double btheta, - int k, - arma::vec m0, - arma::mat V0, - double sigma2_a0, - double sigma2_b0, - arma::vec beta0, - double sigma20, - double eta0, - arma::vec lambda0, - double const& variance, - double ar, - arma::vec LSmu0, - arma::vec LSdelta0, - double LSphi0, - arma::vec LSnu0, - arma::vec LStheta0, - arma::vec sumByCellBio, - arma::vec sumByGeneAll, - arma::vec sumByGeneBio, - int StoreAdapt, - int EndAdapt, - int PrintProgress, - bool FixLocations, - bool RBFMinMax, - arma::vec RBFLocations, - double const& mintol_mu, - double const& mintol_delta, - double const& mintol_nu, - double const& mintol_theta, - double const& geneExponent, - double const& cellExponent); - - -// MCMC sampler for the non-spike case -// [[Rcpp::export(".BASiCS_MCMCcppNoSpikes")]] -Rcpp::List BASiCS_MCMCcppNoSpikes( - int N, - int Thin, - int Burn, - arma::mat Counts, - arma::mat BatchDesign, - arma::vec mu0, - arma::vec delta0, - arma::vec s0, - arma::vec nu0, - arma::vec theta0, - arma::vec mu_mu, - double s2mu, - double adelta, - double bdelta, - double s2delta, - double prior_delta, - double as, - double bs, - double atheta, - double btheta, - double Constrain, - arma::vec Index, - int RefGene, - arma::vec RefGenes, - arma::vec ConstrainGene, - arma::vec NotConstrainGene, - int StochasticRef, - double ar, - arma::vec LSmu0, - arma::vec LSdelta0, - arma::vec LSnu0, - arma::vec LStheta0, - arma::vec sumByCellAll, - arma::vec sumByGeneAll, - int StoreAdapt, - int EndAdapt, - int PrintProgress, - double const& mintol_mu, - double const& mintol_delta, - double const& mintol_nu, - double const& mintol_theta, - double const& geneExponent, - double const& cellExponent); - -// MCMC sampler for regression and non-spikes case -// [[Rcpp::export(".BASiCS_MCMCcppRegNoSpikes")]] -Rcpp::List BASiCS_MCMCcppRegNoSpikes( - int N, - int Thin, - int Burn, - arma::mat Counts, - arma::mat BatchDesign, - arma::vec mu0, - arma::vec delta0, - arma::vec s0, - arma::vec nu0, - arma::vec theta0, - arma::vec mu_mu, - double s2mu, - double as, - double bs, - double atheta, - double btheta, - int k, - arma::vec m0, - arma::mat V0, - double sigma2_a0, - double sigma2_b0, - arma::vec beta0, - double sigma20, - double eta0, - arma::vec lambda0, - double const& variance, - double Constrain, - arma::vec Index, - int RefGene, - arma::vec RefGenes, - arma::vec ConstrainGene, - arma::vec NotConstrainGene, - int StochasticRef, - double ar, - arma::vec LSmu0, - arma::vec LSdelta0, - arma::vec LSnu0, - arma::vec LStheta0, - arma::vec sumByCellAll, - arma::vec sumByGeneAll, - int StoreAdapt, - int EndAdapt, - int PrintProgress, - bool RBFMinMax, - bool FixLocations, - arma::vec RBFLocations, - double const& mintol_mu, - double const& mintol_delta, - double const& mintol_nu, - double const& mintol_theta, - double const& geneExponent, - double const& cellExponent); - -// Function to cumpute denoised rates -// [[Rcpp::export(".BASiCS_DenoisedRates")]] -arma::mat BASiCS_DenoisedRates( - NumericMatrix CountsBio, - NumericMatrix Mu, - NumericMatrix TransInvDelta, - NumericMatrix PhiNu, - int N, - int q0, - int n); - -#endif diff --git a/src/Makevars b/src/Makevars index d198bd1f..b34835a3 100644 --- a/src/Makevars +++ b/src/Makevars @@ -1,5 +1,5 @@ ## This assume that we can call Rscript to ask Rcpp about its locations ## Use the R_HOME indirection to support installations of multiple R version -PKG_LIBS=$(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -#PKG_LIBS = $(shell $(R_HOME)/bin/Rscript.exe -e "Rcpp:::LdFlags()") $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -#$(shell $(R_HOME)/bin/Rscript -e "RcppArmadillo:::CxxFlags()") PKG_CXXFLAGS = $(PKG_CXXFLAGS) $(shell $(R_HOME)/bin/Rscript -e "RcppArmadillo:::CxxFlags()") +CXX_STD = CXX11 +PKG_LIBS=$(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) $(SHLIB_OPENMP_CXXFLAGS) +PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) diff --git a/src/Makevars.win b/src/Makevars.win index dfd835eb..049b221a 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -1,6 +1,4 @@ ## This assume that we can call Rscript to ask Rcpp about its locations ## Use the R_HOME indirection to support installations of multiple R version -PKG_LIBS=$(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -#PKG_LIBS = $(shell $(R_HOME)/bin/Rscript.exe -e "Rcpp:::LdFlags()") $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -#$(shell $(R_HOME)/bin/Rscript -e "RcppArmadillo:::CxxFlags()") PKG_CXXFLAGS = $(PKG_CXXFLAGS) $(shell $(R_HOME)/bin/Rscript -e "RcppArmadillo:::CxxFlags()") - +CXX_STD = CXX11 +PKG_LIBS=$(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) $(SHLIB_OPENMP_CXXFLAGS) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 66b00f9e..49a0868e 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -6,9 +6,26 @@ using namespace Rcpp; +// BASiCS_DenoisedRates +arma::mat BASiCS_DenoisedRates(NumericMatrix CountsBio, NumericMatrix Mu, NumericMatrix TransInvDelta, NumericMatrix PhiNu, int N, int q0, int n); +RcppExport SEXP _BASiCS_BASiCS_DenoisedRates(SEXP CountsBioSEXP, SEXP MuSEXP, SEXP TransInvDeltaSEXP, SEXP PhiNuSEXP, SEXP NSEXP, SEXP q0SEXP, SEXP nSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericMatrix >::type CountsBio(CountsBioSEXP); + Rcpp::traits::input_parameter< NumericMatrix >::type Mu(MuSEXP); + Rcpp::traits::input_parameter< NumericMatrix >::type TransInvDelta(TransInvDeltaSEXP); + Rcpp::traits::input_parameter< NumericMatrix >::type PhiNu(PhiNuSEXP); + Rcpp::traits::input_parameter< int >::type N(NSEXP); + Rcpp::traits::input_parameter< int >::type q0(q0SEXP); + Rcpp::traits::input_parameter< int >::type n(nSEXP); + rcpp_result_gen = Rcpp::wrap(BASiCS_DenoisedRates(CountsBio, Mu, TransInvDelta, PhiNu, N, q0, n)); + return rcpp_result_gen; +END_RCPP +} // BASiCS_MCMCcpp -Rcpp::List BASiCS_MCMCcpp(int N, int Thin, int Burn, arma::mat Counts, arma::mat BatchDesign, arma::vec muSpikes, arma::vec mu0, arma::vec delta0, arma::vec phi0, arma::vec s0, arma::vec nu0, arma::vec theta0, arma::vec mu_mu, double s2mu, double adelta, double bdelta, double s2delta, double prior_delta, arma::vec aphi, double as, double bs, double atheta, double btheta, double ar, arma::vec LSmu0, arma::vec LSdelta0, double LSphi0, arma::vec LSnu0, arma::vec LStheta0, arma::vec sumByCellBio, arma::vec sumByGeneAll, arma::vec sumByGeneBio, int StoreAdapt, int EndAdapt, int PrintProgress, double const& mintol_mu, double const& mintol_delta, double const& mintol_nu, double const& mintol_theta, double const& geneExponent, double const& cellExponent); -RcppExport SEXP _BASiCS_BASiCS_MCMCcpp(SEXP NSEXP, SEXP ThinSEXP, SEXP BurnSEXP, SEXP CountsSEXP, SEXP BatchDesignSEXP, SEXP muSpikesSEXP, SEXP mu0SEXP, SEXP delta0SEXP, SEXP phi0SEXP, SEXP s0SEXP, SEXP nu0SEXP, SEXP theta0SEXP, SEXP mu_muSEXP, SEXP s2muSEXP, SEXP adeltaSEXP, SEXP bdeltaSEXP, SEXP s2deltaSEXP, SEXP prior_deltaSEXP, SEXP aphiSEXP, SEXP asSEXP, SEXP bsSEXP, SEXP athetaSEXP, SEXP bthetaSEXP, SEXP arSEXP, SEXP LSmu0SEXP, SEXP LSdelta0SEXP, SEXP LSphi0SEXP, SEXP LSnu0SEXP, SEXP LStheta0SEXP, SEXP sumByCellBioSEXP, SEXP sumByGeneAllSEXP, SEXP sumByGeneBioSEXP, SEXP StoreAdaptSEXP, SEXP EndAdaptSEXP, SEXP PrintProgressSEXP, SEXP mintol_muSEXP, SEXP mintol_deltaSEXP, SEXP mintol_nuSEXP, SEXP mintol_thetaSEXP, SEXP geneExponentSEXP, SEXP cellExponentSEXP) { +Rcpp::List BASiCS_MCMCcpp(int N, int Thin, int Burn, arma::mat Counts, arma::mat BatchDesign, arma::vec muSpikes, arma::vec mu0, arma::vec delta0, arma::vec phi0, arma::vec s0, arma::vec nu0, arma::vec theta0, arma::vec mu_mu, double s2mu, double adelta, double bdelta, double s2delta, double prior_delta, arma::vec aphi, double as, double bs, double atheta, double btheta, double ar, arma::vec LSmu0, arma::vec LSdelta0, double LSphi0, arma::vec LSnu0, arma::vec LStheta0, arma::vec sumByCellBio, arma::vec sumByGeneAll, arma::vec sumByGeneBio, int StoreAdapt, int EndAdapt, int PrintProgress, double const& mintol_mu, double const& mintol_delta, double const& mintol_nu, double const& mintol_theta, double const& geneExponent, double const& cellExponent, int threads); +RcppExport SEXP _BASiCS_BASiCS_MCMCcpp(SEXP NSEXP, SEXP ThinSEXP, SEXP BurnSEXP, SEXP CountsSEXP, SEXP BatchDesignSEXP, SEXP muSpikesSEXP, SEXP mu0SEXP, SEXP delta0SEXP, SEXP phi0SEXP, SEXP s0SEXP, SEXP nu0SEXP, SEXP theta0SEXP, SEXP mu_muSEXP, SEXP s2muSEXP, SEXP adeltaSEXP, SEXP bdeltaSEXP, SEXP s2deltaSEXP, SEXP prior_deltaSEXP, SEXP aphiSEXP, SEXP asSEXP, SEXP bsSEXP, SEXP athetaSEXP, SEXP bthetaSEXP, SEXP arSEXP, SEXP LSmu0SEXP, SEXP LSdelta0SEXP, SEXP LSphi0SEXP, SEXP LSnu0SEXP, SEXP LStheta0SEXP, SEXP sumByCellBioSEXP, SEXP sumByGeneAllSEXP, SEXP sumByGeneBioSEXP, SEXP StoreAdaptSEXP, SEXP EndAdaptSEXP, SEXP PrintProgressSEXP, SEXP mintol_muSEXP, SEXP mintol_deltaSEXP, SEXP mintol_nuSEXP, SEXP mintol_thetaSEXP, SEXP geneExponentSEXP, SEXP cellExponentSEXP, SEXP threadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -53,13 +70,14 @@ BEGIN_RCPP Rcpp::traits::input_parameter< double const& >::type mintol_theta(mintol_thetaSEXP); Rcpp::traits::input_parameter< double const& >::type geneExponent(geneExponentSEXP); Rcpp::traits::input_parameter< double const& >::type cellExponent(cellExponentSEXP); - rcpp_result_gen = Rcpp::wrap(BASiCS_MCMCcpp(N, Thin, Burn, Counts, BatchDesign, muSpikes, mu0, delta0, phi0, s0, nu0, theta0, mu_mu, s2mu, adelta, bdelta, s2delta, prior_delta, aphi, as, bs, atheta, btheta, ar, LSmu0, LSdelta0, LSphi0, LSnu0, LStheta0, sumByCellBio, sumByGeneAll, sumByGeneBio, StoreAdapt, EndAdapt, PrintProgress, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent)); + Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); + rcpp_result_gen = Rcpp::wrap(BASiCS_MCMCcpp(N, Thin, Burn, Counts, BatchDesign, muSpikes, mu0, delta0, phi0, s0, nu0, theta0, mu_mu, s2mu, adelta, bdelta, s2delta, prior_delta, aphi, as, bs, atheta, btheta, ar, LSmu0, LSdelta0, LSphi0, LSnu0, LStheta0, sumByCellBio, sumByGeneAll, sumByGeneBio, StoreAdapt, EndAdapt, PrintProgress, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent, threads)); return rcpp_result_gen; END_RCPP } -// BASiCS_MCMCcppReg -Rcpp::List BASiCS_MCMCcppReg(int N, int Thin, int Burn, arma::mat Counts, arma::mat BatchDesign, arma::vec muSpikes, arma::vec mu0, arma::vec delta0, arma::vec phi0, arma::vec s0, arma::vec nu0, arma::vec theta0, arma::vec mu_mu, double s2mu, arma::vec aphi, double as, double bs, double atheta, double btheta, int k, arma::vec m0, arma::mat V0, double sigma2_a0, double sigma2_b0, arma::vec beta0, double sigma20, double eta0, arma::vec lambda0, double const& variance, double ar, arma::vec LSmu0, arma::vec LSdelta0, double LSphi0, arma::vec LSnu0, arma::vec LStheta0, arma::vec sumByCellBio, arma::vec sumByGeneAll, arma::vec sumByGeneBio, int StoreAdapt, int EndAdapt, int PrintProgress, bool FixLocations, bool RBFMinMax, arma::vec RBFLocations, double const& mintol_mu, double const& mintol_delta, double const& mintol_nu, double const& mintol_theta, double const& geneExponent, double const& cellExponent); -RcppExport SEXP _BASiCS_BASiCS_MCMCcppReg(SEXP NSEXP, SEXP ThinSEXP, SEXP BurnSEXP, SEXP CountsSEXP, SEXP BatchDesignSEXP, SEXP muSpikesSEXP, SEXP mu0SEXP, SEXP delta0SEXP, SEXP phi0SEXP, SEXP s0SEXP, SEXP nu0SEXP, SEXP theta0SEXP, SEXP mu_muSEXP, SEXP s2muSEXP, SEXP aphiSEXP, SEXP asSEXP, SEXP bsSEXP, SEXP athetaSEXP, SEXP bthetaSEXP, SEXP kSEXP, SEXP m0SEXP, SEXP V0SEXP, SEXP sigma2_a0SEXP, SEXP sigma2_b0SEXP, SEXP beta0SEXP, SEXP sigma20SEXP, SEXP eta0SEXP, SEXP lambda0SEXP, SEXP varianceSEXP, SEXP arSEXP, SEXP LSmu0SEXP, SEXP LSdelta0SEXP, SEXP LSphi0SEXP, SEXP LSnu0SEXP, SEXP LStheta0SEXP, SEXP sumByCellBioSEXP, SEXP sumByGeneAllSEXP, SEXP sumByGeneBioSEXP, SEXP StoreAdaptSEXP, SEXP EndAdaptSEXP, SEXP PrintProgressSEXP, SEXP FixLocationsSEXP, SEXP RBFMinMaxSEXP, SEXP RBFLocationsSEXP, SEXP mintol_muSEXP, SEXP mintol_deltaSEXP, SEXP mintol_nuSEXP, SEXP mintol_thetaSEXP, SEXP geneExponentSEXP, SEXP cellExponentSEXP) { +// BASiCS_MCMCcppNoSpikes +Rcpp::List BASiCS_MCMCcppNoSpikes(int N, int Thin, int Burn, arma::mat Counts, arma::mat BatchDesign, arma::vec mu0, arma::vec delta0, arma::vec s0, arma::vec nu0, arma::vec theta0, arma::vec mu_mu, double s2mu, double adelta, double bdelta, double s2delta, double prior_delta, double as, double bs, double atheta, double btheta, double Constrain, arma::vec Index, int RefGene, arma::vec RefGenes, arma::vec ConstrainGene, arma::vec NotConstrainGene, int StochasticRef, double ar, arma::vec LSmu0, arma::vec LSdelta0, arma::vec LSnu0, arma::vec LStheta0, arma::vec sumByCellAll, arma::vec sumByGeneAll, int StoreAdapt, int EndAdapt, int PrintProgress, double const& mintol_mu, double const& mintol_delta, double const& mintol_nu, double const& mintol_theta, double const& geneExponent, double const& cellExponent, int threads); +RcppExport SEXP _BASiCS_BASiCS_MCMCcppNoSpikes(SEXP NSEXP, SEXP ThinSEXP, SEXP BurnSEXP, SEXP CountsSEXP, SEXP BatchDesignSEXP, SEXP mu0SEXP, SEXP delta0SEXP, SEXP s0SEXP, SEXP nu0SEXP, SEXP theta0SEXP, SEXP mu_muSEXP, SEXP s2muSEXP, SEXP adeltaSEXP, SEXP bdeltaSEXP, SEXP s2deltaSEXP, SEXP prior_deltaSEXP, SEXP asSEXP, SEXP bsSEXP, SEXP athetaSEXP, SEXP bthetaSEXP, SEXP ConstrainSEXP, SEXP IndexSEXP, SEXP RefGeneSEXP, SEXP RefGenesSEXP, SEXP ConstrainGeneSEXP, SEXP NotConstrainGeneSEXP, SEXP StochasticRefSEXP, SEXP arSEXP, SEXP LSmu0SEXP, SEXP LSdelta0SEXP, SEXP LSnu0SEXP, SEXP LStheta0SEXP, SEXP sumByCellAllSEXP, SEXP sumByGeneAllSEXP, SEXP StoreAdaptSEXP, SEXP EndAdaptSEXP, SEXP PrintProgressSEXP, SEXP mintol_muSEXP, SEXP mintol_deltaSEXP, SEXP mintol_nuSEXP, SEXP mintol_thetaSEXP, SEXP geneExponentSEXP, SEXP cellExponentSEXP, SEXP threadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -68,58 +86,52 @@ BEGIN_RCPP Rcpp::traits::input_parameter< int >::type Burn(BurnSEXP); Rcpp::traits::input_parameter< arma::mat >::type Counts(CountsSEXP); Rcpp::traits::input_parameter< arma::mat >::type BatchDesign(BatchDesignSEXP); - Rcpp::traits::input_parameter< arma::vec >::type muSpikes(muSpikesSEXP); Rcpp::traits::input_parameter< arma::vec >::type mu0(mu0SEXP); Rcpp::traits::input_parameter< arma::vec >::type delta0(delta0SEXP); - Rcpp::traits::input_parameter< arma::vec >::type phi0(phi0SEXP); Rcpp::traits::input_parameter< arma::vec >::type s0(s0SEXP); Rcpp::traits::input_parameter< arma::vec >::type nu0(nu0SEXP); Rcpp::traits::input_parameter< arma::vec >::type theta0(theta0SEXP); Rcpp::traits::input_parameter< arma::vec >::type mu_mu(mu_muSEXP); Rcpp::traits::input_parameter< double >::type s2mu(s2muSEXP); - Rcpp::traits::input_parameter< arma::vec >::type aphi(aphiSEXP); + Rcpp::traits::input_parameter< double >::type adelta(adeltaSEXP); + Rcpp::traits::input_parameter< double >::type bdelta(bdeltaSEXP); + Rcpp::traits::input_parameter< double >::type s2delta(s2deltaSEXP); + Rcpp::traits::input_parameter< double >::type prior_delta(prior_deltaSEXP); Rcpp::traits::input_parameter< double >::type as(asSEXP); Rcpp::traits::input_parameter< double >::type bs(bsSEXP); Rcpp::traits::input_parameter< double >::type atheta(athetaSEXP); Rcpp::traits::input_parameter< double >::type btheta(bthetaSEXP); - Rcpp::traits::input_parameter< int >::type k(kSEXP); - Rcpp::traits::input_parameter< arma::vec >::type m0(m0SEXP); - Rcpp::traits::input_parameter< arma::mat >::type V0(V0SEXP); - Rcpp::traits::input_parameter< double >::type sigma2_a0(sigma2_a0SEXP); - Rcpp::traits::input_parameter< double >::type sigma2_b0(sigma2_b0SEXP); - Rcpp::traits::input_parameter< arma::vec >::type beta0(beta0SEXP); - Rcpp::traits::input_parameter< double >::type sigma20(sigma20SEXP); - Rcpp::traits::input_parameter< double >::type eta0(eta0SEXP); - Rcpp::traits::input_parameter< arma::vec >::type lambda0(lambda0SEXP); - Rcpp::traits::input_parameter< double const& >::type variance(varianceSEXP); + Rcpp::traits::input_parameter< double >::type Constrain(ConstrainSEXP); + Rcpp::traits::input_parameter< arma::vec >::type Index(IndexSEXP); + Rcpp::traits::input_parameter< int >::type RefGene(RefGeneSEXP); + Rcpp::traits::input_parameter< arma::vec >::type RefGenes(RefGenesSEXP); + Rcpp::traits::input_parameter< arma::vec >::type ConstrainGene(ConstrainGeneSEXP); + Rcpp::traits::input_parameter< arma::vec >::type NotConstrainGene(NotConstrainGeneSEXP); + Rcpp::traits::input_parameter< int >::type StochasticRef(StochasticRefSEXP); Rcpp::traits::input_parameter< double >::type ar(arSEXP); Rcpp::traits::input_parameter< arma::vec >::type LSmu0(LSmu0SEXP); Rcpp::traits::input_parameter< arma::vec >::type LSdelta0(LSdelta0SEXP); - Rcpp::traits::input_parameter< double >::type LSphi0(LSphi0SEXP); Rcpp::traits::input_parameter< arma::vec >::type LSnu0(LSnu0SEXP); Rcpp::traits::input_parameter< arma::vec >::type LStheta0(LStheta0SEXP); - Rcpp::traits::input_parameter< arma::vec >::type sumByCellBio(sumByCellBioSEXP); + Rcpp::traits::input_parameter< arma::vec >::type sumByCellAll(sumByCellAllSEXP); Rcpp::traits::input_parameter< arma::vec >::type sumByGeneAll(sumByGeneAllSEXP); - Rcpp::traits::input_parameter< arma::vec >::type sumByGeneBio(sumByGeneBioSEXP); Rcpp::traits::input_parameter< int >::type StoreAdapt(StoreAdaptSEXP); Rcpp::traits::input_parameter< int >::type EndAdapt(EndAdaptSEXP); Rcpp::traits::input_parameter< int >::type PrintProgress(PrintProgressSEXP); - Rcpp::traits::input_parameter< bool >::type FixLocations(FixLocationsSEXP); - Rcpp::traits::input_parameter< bool >::type RBFMinMax(RBFMinMaxSEXP); - Rcpp::traits::input_parameter< arma::vec >::type RBFLocations(RBFLocationsSEXP); Rcpp::traits::input_parameter< double const& >::type mintol_mu(mintol_muSEXP); Rcpp::traits::input_parameter< double const& >::type mintol_delta(mintol_deltaSEXP); Rcpp::traits::input_parameter< double const& >::type mintol_nu(mintol_nuSEXP); Rcpp::traits::input_parameter< double const& >::type mintol_theta(mintol_thetaSEXP); Rcpp::traits::input_parameter< double const& >::type geneExponent(geneExponentSEXP); Rcpp::traits::input_parameter< double const& >::type cellExponent(cellExponentSEXP); - rcpp_result_gen = Rcpp::wrap(BASiCS_MCMCcppReg(N, Thin, Burn, Counts, BatchDesign, muSpikes, mu0, delta0, phi0, s0, nu0, theta0, mu_mu, s2mu, aphi, as, bs, atheta, btheta, k, m0, V0, sigma2_a0, sigma2_b0, beta0, sigma20, eta0, lambda0, variance, ar, LSmu0, LSdelta0, LSphi0, LSnu0, LStheta0, sumByCellBio, sumByGeneAll, sumByGeneBio, StoreAdapt, EndAdapt, PrintProgress, FixLocations, RBFMinMax, RBFLocations, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent)); + Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); + rcpp_result_gen = Rcpp::wrap(BASiCS_MCMCcppNoSpikes(N, Thin, Burn, Counts, BatchDesign, mu0, delta0, s0, nu0, theta0, mu_mu, s2mu, adelta, bdelta, s2delta, prior_delta, as, bs, atheta, btheta, Constrain, Index, RefGene, RefGenes, ConstrainGene, NotConstrainGene, StochasticRef, ar, LSmu0, LSdelta0, LSnu0, LStheta0, sumByCellAll, sumByGeneAll, StoreAdapt, EndAdapt, PrintProgress, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent, threads)); return rcpp_result_gen; END_RCPP } -// BASiCS_MCMCcppNoSpikes -Rcpp::List BASiCS_MCMCcppNoSpikes(int N, int Thin, int Burn, arma::mat Counts, arma::mat BatchDesign, arma::vec mu0, arma::vec delta0, arma::vec s0, arma::vec nu0, arma::vec theta0, arma::vec mu_mu, double s2mu, double adelta, double bdelta, double s2delta, double prior_delta, double as, double bs, double atheta, double btheta, double Constrain, arma::vec Index, int RefGene, arma::vec RefGenes, arma::vec ConstrainGene, arma::vec NotConstrainGene, int StochasticRef, double ar, arma::vec LSmu0, arma::vec LSdelta0, arma::vec LSnu0, arma::vec LStheta0, arma::vec sumByCellAll, arma::vec sumByGeneAll, int StoreAdapt, int EndAdapt, int PrintProgress, double const& mintol_mu, double const& mintol_delta, double const& mintol_nu, double const& mintol_theta, double const& geneExponent, double const& cellExponent); -RcppExport SEXP _BASiCS_BASiCS_MCMCcppNoSpikes(SEXP NSEXP, SEXP ThinSEXP, SEXP BurnSEXP, SEXP CountsSEXP, SEXP BatchDesignSEXP, SEXP mu0SEXP, SEXP delta0SEXP, SEXP s0SEXP, SEXP nu0SEXP, SEXP theta0SEXP, SEXP mu_muSEXP, SEXP s2muSEXP, SEXP adeltaSEXP, SEXP bdeltaSEXP, SEXP s2deltaSEXP, SEXP prior_deltaSEXP, SEXP asSEXP, SEXP bsSEXP, SEXP athetaSEXP, SEXP bthetaSEXP, SEXP ConstrainSEXP, SEXP IndexSEXP, SEXP RefGeneSEXP, SEXP RefGenesSEXP, SEXP ConstrainGeneSEXP, SEXP NotConstrainGeneSEXP, SEXP StochasticRefSEXP, SEXP arSEXP, SEXP LSmu0SEXP, SEXP LSdelta0SEXP, SEXP LSnu0SEXP, SEXP LStheta0SEXP, SEXP sumByCellAllSEXP, SEXP sumByGeneAllSEXP, SEXP StoreAdaptSEXP, SEXP EndAdaptSEXP, SEXP PrintProgressSEXP, SEXP mintol_muSEXP, SEXP mintol_deltaSEXP, SEXP mintol_nuSEXP, SEXP mintol_thetaSEXP, SEXP geneExponentSEXP, SEXP cellExponentSEXP) { +// BASiCS_MCMCcppReg +Rcpp::List BASiCS_MCMCcppReg(int N, int Thin, int Burn, arma::mat Counts, arma::mat BatchDesign, arma::vec muSpikes, arma::vec mu0, arma::vec delta0, arma::vec phi0, arma::vec s0, arma::vec nu0, arma::vec theta0, arma::vec mu_mu, double s2mu, arma::vec aphi, double as, double bs, double atheta, double btheta, int k, arma::vec m0, arma::mat V0, double sigma2_a0, double sigma2_b0, arma::vec beta0, double sigma20, double eta0, arma::vec lambda0, double const& variance, double ar, arma::vec LSmu0, arma::vec LSdelta0, double LSphi0, arma::vec LSnu0, arma::vec LStheta0, arma::vec sumByCellBio, arma::vec sumByGeneAll, arma::vec sumByGeneBio, int StoreAdapt, int EndAdapt, int PrintProgress, bool FixLocations, bool RBFMinMax, arma::vec RBFLocations, double const& mintol_mu, double const& mintol_delta, double const& mintol_nu, double const& mintol_theta, double const& geneExponent, double const& cellExponent, int threads); +RcppExport SEXP _BASiCS_BASiCS_MCMCcppReg(SEXP NSEXP, SEXP ThinSEXP, SEXP BurnSEXP, SEXP CountsSEXP, SEXP BatchDesignSEXP, SEXP muSpikesSEXP, SEXP mu0SEXP, SEXP delta0SEXP, SEXP phi0SEXP, SEXP s0SEXP, SEXP nu0SEXP, SEXP theta0SEXP, SEXP mu_muSEXP, SEXP s2muSEXP, SEXP aphiSEXP, SEXP asSEXP, SEXP bsSEXP, SEXP athetaSEXP, SEXP bthetaSEXP, SEXP kSEXP, SEXP m0SEXP, SEXP V0SEXP, SEXP sigma2_a0SEXP, SEXP sigma2_b0SEXP, SEXP beta0SEXP, SEXP sigma20SEXP, SEXP eta0SEXP, SEXP lambda0SEXP, SEXP varianceSEXP, SEXP arSEXP, SEXP LSmu0SEXP, SEXP LSdelta0SEXP, SEXP LSphi0SEXP, SEXP LSnu0SEXP, SEXP LStheta0SEXP, SEXP sumByCellBioSEXP, SEXP sumByGeneAllSEXP, SEXP sumByGeneBioSEXP, SEXP StoreAdaptSEXP, SEXP EndAdaptSEXP, SEXP PrintProgressSEXP, SEXP FixLocationsSEXP, SEXP RBFMinMaxSEXP, SEXP RBFLocationsSEXP, SEXP mintol_muSEXP, SEXP mintol_deltaSEXP, SEXP mintol_nuSEXP, SEXP mintol_thetaSEXP, SEXP geneExponentSEXP, SEXP cellExponentSEXP, SEXP threadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -128,51 +140,59 @@ BEGIN_RCPP Rcpp::traits::input_parameter< int >::type Burn(BurnSEXP); Rcpp::traits::input_parameter< arma::mat >::type Counts(CountsSEXP); Rcpp::traits::input_parameter< arma::mat >::type BatchDesign(BatchDesignSEXP); + Rcpp::traits::input_parameter< arma::vec >::type muSpikes(muSpikesSEXP); Rcpp::traits::input_parameter< arma::vec >::type mu0(mu0SEXP); Rcpp::traits::input_parameter< arma::vec >::type delta0(delta0SEXP); + Rcpp::traits::input_parameter< arma::vec >::type phi0(phi0SEXP); Rcpp::traits::input_parameter< arma::vec >::type s0(s0SEXP); Rcpp::traits::input_parameter< arma::vec >::type nu0(nu0SEXP); Rcpp::traits::input_parameter< arma::vec >::type theta0(theta0SEXP); Rcpp::traits::input_parameter< arma::vec >::type mu_mu(mu_muSEXP); Rcpp::traits::input_parameter< double >::type s2mu(s2muSEXP); - Rcpp::traits::input_parameter< double >::type adelta(adeltaSEXP); - Rcpp::traits::input_parameter< double >::type bdelta(bdeltaSEXP); - Rcpp::traits::input_parameter< double >::type s2delta(s2deltaSEXP); - Rcpp::traits::input_parameter< double >::type prior_delta(prior_deltaSEXP); + Rcpp::traits::input_parameter< arma::vec >::type aphi(aphiSEXP); Rcpp::traits::input_parameter< double >::type as(asSEXP); Rcpp::traits::input_parameter< double >::type bs(bsSEXP); Rcpp::traits::input_parameter< double >::type atheta(athetaSEXP); Rcpp::traits::input_parameter< double >::type btheta(bthetaSEXP); - Rcpp::traits::input_parameter< double >::type Constrain(ConstrainSEXP); - Rcpp::traits::input_parameter< arma::vec >::type Index(IndexSEXP); - Rcpp::traits::input_parameter< int >::type RefGene(RefGeneSEXP); - Rcpp::traits::input_parameter< arma::vec >::type RefGenes(RefGenesSEXP); - Rcpp::traits::input_parameter< arma::vec >::type ConstrainGene(ConstrainGeneSEXP); - Rcpp::traits::input_parameter< arma::vec >::type NotConstrainGene(NotConstrainGeneSEXP); - Rcpp::traits::input_parameter< int >::type StochasticRef(StochasticRefSEXP); + Rcpp::traits::input_parameter< int >::type k(kSEXP); + Rcpp::traits::input_parameter< arma::vec >::type m0(m0SEXP); + Rcpp::traits::input_parameter< arma::mat >::type V0(V0SEXP); + Rcpp::traits::input_parameter< double >::type sigma2_a0(sigma2_a0SEXP); + Rcpp::traits::input_parameter< double >::type sigma2_b0(sigma2_b0SEXP); + Rcpp::traits::input_parameter< arma::vec >::type beta0(beta0SEXP); + Rcpp::traits::input_parameter< double >::type sigma20(sigma20SEXP); + Rcpp::traits::input_parameter< double >::type eta0(eta0SEXP); + Rcpp::traits::input_parameter< arma::vec >::type lambda0(lambda0SEXP); + Rcpp::traits::input_parameter< double const& >::type variance(varianceSEXP); Rcpp::traits::input_parameter< double >::type ar(arSEXP); Rcpp::traits::input_parameter< arma::vec >::type LSmu0(LSmu0SEXP); Rcpp::traits::input_parameter< arma::vec >::type LSdelta0(LSdelta0SEXP); + Rcpp::traits::input_parameter< double >::type LSphi0(LSphi0SEXP); Rcpp::traits::input_parameter< arma::vec >::type LSnu0(LSnu0SEXP); Rcpp::traits::input_parameter< arma::vec >::type LStheta0(LStheta0SEXP); - Rcpp::traits::input_parameter< arma::vec >::type sumByCellAll(sumByCellAllSEXP); + Rcpp::traits::input_parameter< arma::vec >::type sumByCellBio(sumByCellBioSEXP); Rcpp::traits::input_parameter< arma::vec >::type sumByGeneAll(sumByGeneAllSEXP); + Rcpp::traits::input_parameter< arma::vec >::type sumByGeneBio(sumByGeneBioSEXP); Rcpp::traits::input_parameter< int >::type StoreAdapt(StoreAdaptSEXP); Rcpp::traits::input_parameter< int >::type EndAdapt(EndAdaptSEXP); Rcpp::traits::input_parameter< int >::type PrintProgress(PrintProgressSEXP); + Rcpp::traits::input_parameter< bool >::type FixLocations(FixLocationsSEXP); + Rcpp::traits::input_parameter< bool >::type RBFMinMax(RBFMinMaxSEXP); + Rcpp::traits::input_parameter< arma::vec >::type RBFLocations(RBFLocationsSEXP); Rcpp::traits::input_parameter< double const& >::type mintol_mu(mintol_muSEXP); Rcpp::traits::input_parameter< double const& >::type mintol_delta(mintol_deltaSEXP); Rcpp::traits::input_parameter< double const& >::type mintol_nu(mintol_nuSEXP); Rcpp::traits::input_parameter< double const& >::type mintol_theta(mintol_thetaSEXP); Rcpp::traits::input_parameter< double const& >::type geneExponent(geneExponentSEXP); Rcpp::traits::input_parameter< double const& >::type cellExponent(cellExponentSEXP); - rcpp_result_gen = Rcpp::wrap(BASiCS_MCMCcppNoSpikes(N, Thin, Burn, Counts, BatchDesign, mu0, delta0, s0, nu0, theta0, mu_mu, s2mu, adelta, bdelta, s2delta, prior_delta, as, bs, atheta, btheta, Constrain, Index, RefGene, RefGenes, ConstrainGene, NotConstrainGene, StochasticRef, ar, LSmu0, LSdelta0, LSnu0, LStheta0, sumByCellAll, sumByGeneAll, StoreAdapt, EndAdapt, PrintProgress, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent)); + Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); + rcpp_result_gen = Rcpp::wrap(BASiCS_MCMCcppReg(N, Thin, Burn, Counts, BatchDesign, muSpikes, mu0, delta0, phi0, s0, nu0, theta0, mu_mu, s2mu, aphi, as, bs, atheta, btheta, k, m0, V0, sigma2_a0, sigma2_b0, beta0, sigma20, eta0, lambda0, variance, ar, LSmu0, LSdelta0, LSphi0, LSnu0, LStheta0, sumByCellBio, sumByGeneAll, sumByGeneBio, StoreAdapt, EndAdapt, PrintProgress, FixLocations, RBFMinMax, RBFLocations, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent, threads)); return rcpp_result_gen; END_RCPP } // BASiCS_MCMCcppRegNoSpikes -Rcpp::List BASiCS_MCMCcppRegNoSpikes(int N, int Thin, int Burn, arma::mat Counts, arma::mat BatchDesign, arma::vec mu0, arma::vec delta0, arma::vec s0, arma::vec nu0, arma::vec theta0, arma::vec mu_mu, double s2mu, double as, double bs, double atheta, double btheta, int k, arma::vec m0, arma::mat V0, double sigma2_a0, double sigma2_b0, arma::vec beta0, double sigma20, double eta0, arma::vec lambda0, double const& variance, double Constrain, arma::vec Index, int RefGene, arma::vec RefGenes, arma::vec ConstrainGene, arma::vec NotConstrainGene, int StochasticRef, double ar, arma::vec LSmu0, arma::vec LSdelta0, arma::vec LSnu0, arma::vec LStheta0, arma::vec sumByCellAll, arma::vec sumByGeneAll, int StoreAdapt, int EndAdapt, int PrintProgress, bool RBFMinMax, bool FixLocations, arma::vec RBFLocations, double const& mintol_mu, double const& mintol_delta, double const& mintol_nu, double const& mintol_theta, double const& geneExponent, double const& cellExponent); -RcppExport SEXP _BASiCS_BASiCS_MCMCcppRegNoSpikes(SEXP NSEXP, SEXP ThinSEXP, SEXP BurnSEXP, SEXP CountsSEXP, SEXP BatchDesignSEXP, SEXP mu0SEXP, SEXP delta0SEXP, SEXP s0SEXP, SEXP nu0SEXP, SEXP theta0SEXP, SEXP mu_muSEXP, SEXP s2muSEXP, SEXP asSEXP, SEXP bsSEXP, SEXP athetaSEXP, SEXP bthetaSEXP, SEXP kSEXP, SEXP m0SEXP, SEXP V0SEXP, SEXP sigma2_a0SEXP, SEXP sigma2_b0SEXP, SEXP beta0SEXP, SEXP sigma20SEXP, SEXP eta0SEXP, SEXP lambda0SEXP, SEXP varianceSEXP, SEXP ConstrainSEXP, SEXP IndexSEXP, SEXP RefGeneSEXP, SEXP RefGenesSEXP, SEXP ConstrainGeneSEXP, SEXP NotConstrainGeneSEXP, SEXP StochasticRefSEXP, SEXP arSEXP, SEXP LSmu0SEXP, SEXP LSdelta0SEXP, SEXP LSnu0SEXP, SEXP LStheta0SEXP, SEXP sumByCellAllSEXP, SEXP sumByGeneAllSEXP, SEXP StoreAdaptSEXP, SEXP EndAdaptSEXP, SEXP PrintProgressSEXP, SEXP RBFMinMaxSEXP, SEXP FixLocationsSEXP, SEXP RBFLocationsSEXP, SEXP mintol_muSEXP, SEXP mintol_deltaSEXP, SEXP mintol_nuSEXP, SEXP mintol_thetaSEXP, SEXP geneExponentSEXP, SEXP cellExponentSEXP) { +Rcpp::List BASiCS_MCMCcppRegNoSpikes(int N, int Thin, int Burn, arma::mat Counts, arma::mat BatchDesign, arma::vec mu0, arma::vec delta0, arma::vec s0, arma::vec nu0, arma::vec theta0, arma::vec mu_mu, double s2mu, double as, double bs, double atheta, double btheta, int k, arma::vec m0, arma::mat V0, double sigma2_a0, double sigma2_b0, arma::vec beta0, double sigma20, double eta0, arma::vec lambda0, double const& variance, double Constrain, arma::vec Index, int RefGene, arma::vec RefGenes, arma::vec ConstrainGene, arma::vec NotConstrainGene, int StochasticRef, double ar, arma::vec LSmu0, arma::vec LSdelta0, arma::vec LSnu0, arma::vec LStheta0, arma::vec sumByCellAll, arma::vec sumByGeneAll, int StoreAdapt, int EndAdapt, int PrintProgress, bool RBFMinMax, bool FixLocations, arma::vec RBFLocations, double const& mintol_mu, double const& mintol_delta, double const& mintol_nu, double const& mintol_theta, double const& geneExponent, double const& cellExponent, int threads); +RcppExport SEXP _BASiCS_BASiCS_MCMCcppRegNoSpikes(SEXP NSEXP, SEXP ThinSEXP, SEXP BurnSEXP, SEXP CountsSEXP, SEXP BatchDesignSEXP, SEXP mu0SEXP, SEXP delta0SEXP, SEXP s0SEXP, SEXP nu0SEXP, SEXP theta0SEXP, SEXP mu_muSEXP, SEXP s2muSEXP, SEXP asSEXP, SEXP bsSEXP, SEXP athetaSEXP, SEXP bthetaSEXP, SEXP kSEXP, SEXP m0SEXP, SEXP V0SEXP, SEXP sigma2_a0SEXP, SEXP sigma2_b0SEXP, SEXP beta0SEXP, SEXP sigma20SEXP, SEXP eta0SEXP, SEXP lambda0SEXP, SEXP varianceSEXP, SEXP ConstrainSEXP, SEXP IndexSEXP, SEXP RefGeneSEXP, SEXP RefGenesSEXP, SEXP ConstrainGeneSEXP, SEXP NotConstrainGeneSEXP, SEXP StochasticRefSEXP, SEXP arSEXP, SEXP LSmu0SEXP, SEXP LSdelta0SEXP, SEXP LSnu0SEXP, SEXP LStheta0SEXP, SEXP sumByCellAllSEXP, SEXP sumByGeneAllSEXP, SEXP StoreAdaptSEXP, SEXP EndAdaptSEXP, SEXP PrintProgressSEXP, SEXP RBFMinMaxSEXP, SEXP FixLocationsSEXP, SEXP RBFLocationsSEXP, SEXP mintol_muSEXP, SEXP mintol_deltaSEXP, SEXP mintol_nuSEXP, SEXP mintol_thetaSEXP, SEXP geneExponentSEXP, SEXP cellExponentSEXP, SEXP threadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -228,48 +248,8 @@ BEGIN_RCPP Rcpp::traits::input_parameter< double const& >::type mintol_theta(mintol_thetaSEXP); Rcpp::traits::input_parameter< double const& >::type geneExponent(geneExponentSEXP); Rcpp::traits::input_parameter< double const& >::type cellExponent(cellExponentSEXP); - rcpp_result_gen = Rcpp::wrap(BASiCS_MCMCcppRegNoSpikes(N, Thin, Burn, Counts, BatchDesign, mu0, delta0, s0, nu0, theta0, mu_mu, s2mu, as, bs, atheta, btheta, k, m0, V0, sigma2_a0, sigma2_b0, beta0, sigma20, eta0, lambda0, variance, Constrain, Index, RefGene, RefGenes, ConstrainGene, NotConstrainGene, StochasticRef, ar, LSmu0, LSdelta0, LSnu0, LStheta0, sumByCellAll, sumByGeneAll, StoreAdapt, EndAdapt, PrintProgress, RBFMinMax, FixLocations, RBFLocations, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent)); - return rcpp_result_gen; -END_RCPP -} -// BASiCS_DenoisedRates -arma::mat BASiCS_DenoisedRates(NumericMatrix CountsBio, NumericMatrix Mu, NumericMatrix TransInvDelta, NumericMatrix PhiNu, int N, int q0, int n); -RcppExport SEXP _BASiCS_BASiCS_DenoisedRates(SEXP CountsBioSEXP, SEXP MuSEXP, SEXP TransInvDeltaSEXP, SEXP PhiNuSEXP, SEXP NSEXP, SEXP q0SEXP, SEXP nSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericMatrix >::type CountsBio(CountsBioSEXP); - Rcpp::traits::input_parameter< NumericMatrix >::type Mu(MuSEXP); - Rcpp::traits::input_parameter< NumericMatrix >::type TransInvDelta(TransInvDeltaSEXP); - Rcpp::traits::input_parameter< NumericMatrix >::type PhiNu(PhiNuSEXP); - Rcpp::traits::input_parameter< int >::type N(NSEXP); - Rcpp::traits::input_parameter< int >::type q0(q0SEXP); - Rcpp::traits::input_parameter< int >::type n(nSEXP); - rcpp_result_gen = Rcpp::wrap(BASiCS_DenoisedRates(CountsBio, Mu, TransInvDelta, PhiNu, N, q0, n)); - return rcpp_result_gen; -END_RCPP -} -// estimateRBFLocations -arma::vec estimateRBFLocations(arma::vec const& log_mu, int const& k, bool RBFMinMax); -RcppExport SEXP _BASiCS_estimateRBFLocations(SEXP log_muSEXP, SEXP kSEXP, SEXP RBFMinMaxSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< arma::vec const& >::type log_mu(log_muSEXP); - Rcpp::traits::input_parameter< int const& >::type k(kSEXP); - Rcpp::traits::input_parameter< bool >::type RBFMinMax(RBFMinMaxSEXP); - rcpp_result_gen = Rcpp::wrap(estimateRBFLocations(log_mu, k, RBFMinMax)); - return rcpp_result_gen; -END_RCPP -} -// rDirichlet -arma::vec rDirichlet(arma::vec alpha); -RcppExport SEXP _BASiCS_rDirichlet(SEXP alphaSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< arma::vec >::type alpha(alphaSEXP); - rcpp_result_gen = Rcpp::wrap(rDirichlet(alpha)); + Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); + rcpp_result_gen = Rcpp::wrap(BASiCS_MCMCcppRegNoSpikes(N, Thin, Burn, Counts, BatchDesign, mu0, delta0, s0, nu0, theta0, mu_mu, s2mu, as, bs, atheta, btheta, k, m0, V0, sigma2_a0, sigma2_b0, beta0, sigma20, eta0, lambda0, variance, Constrain, Index, RefGene, RefGenes, ConstrainGene, NotConstrainGene, StochasticRef, ar, LSmu0, LSdelta0, LSnu0, LStheta0, sumByCellAll, sumByGeneAll, StoreAdapt, EndAdapt, PrintProgress, RBFMinMax, FixLocations, RBFLocations, mintol_mu, mintol_delta, mintol_nu, mintol_theta, geneExponent, cellExponent, threads)); return rcpp_result_gen; END_RCPP } @@ -415,20 +395,88 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// muUpdateNoSpikes +arma::mat muUpdateNoSpikes(arma::vec const& mu0, arma::vec const& prop_var, arma::mat const& Counts, arma::vec const& invdelta, arma::vec const& nu, arma::vec const& sum_bycell_all, arma::vec const& mu_mu, double const& s2_mu, int const& q0, int const& n, arma::vec& mu1, arma::vec& u, arma::vec& ind, double const& Constrain, int const& RefGene, arma::uvec const& ConstrainGene, arma::uvec const& NotConstrainGene, double const& exponent, double const& mintol); +RcppExport SEXP _BASiCS_muUpdateNoSpikes(SEXP mu0SEXP, SEXP prop_varSEXP, SEXP CountsSEXP, SEXP invdeltaSEXP, SEXP nuSEXP, SEXP sum_bycell_allSEXP, SEXP mu_muSEXP, SEXP s2_muSEXP, SEXP q0SEXP, SEXP nSEXP, SEXP mu1SEXP, SEXP uSEXP, SEXP indSEXP, SEXP ConstrainSEXP, SEXP RefGeneSEXP, SEXP ConstrainGeneSEXP, SEXP NotConstrainGeneSEXP, SEXP exponentSEXP, SEXP mintolSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::vec const& >::type mu0(mu0SEXP); + Rcpp::traits::input_parameter< arma::vec const& >::type prop_var(prop_varSEXP); + Rcpp::traits::input_parameter< arma::mat const& >::type Counts(CountsSEXP); + Rcpp::traits::input_parameter< arma::vec const& >::type invdelta(invdeltaSEXP); + Rcpp::traits::input_parameter< arma::vec const& >::type nu(nuSEXP); + Rcpp::traits::input_parameter< arma::vec const& >::type sum_bycell_all(sum_bycell_allSEXP); + Rcpp::traits::input_parameter< arma::vec const& >::type mu_mu(mu_muSEXP); + Rcpp::traits::input_parameter< double const& >::type s2_mu(s2_muSEXP); + Rcpp::traits::input_parameter< int const& >::type q0(q0SEXP); + Rcpp::traits::input_parameter< int const& >::type n(nSEXP); + Rcpp::traits::input_parameter< arma::vec& >::type mu1(mu1SEXP); + Rcpp::traits::input_parameter< arma::vec& >::type u(uSEXP); + Rcpp::traits::input_parameter< arma::vec& >::type ind(indSEXP); + Rcpp::traits::input_parameter< double const& >::type Constrain(ConstrainSEXP); + Rcpp::traits::input_parameter< int const& >::type RefGene(RefGeneSEXP); + Rcpp::traits::input_parameter< arma::uvec const& >::type ConstrainGene(ConstrainGeneSEXP); + Rcpp::traits::input_parameter< arma::uvec const& >::type NotConstrainGene(NotConstrainGeneSEXP); + Rcpp::traits::input_parameter< double const& >::type exponent(exponentSEXP); + Rcpp::traits::input_parameter< double const& >::type mintol(mintolSEXP); + rcpp_result_gen = Rcpp::wrap(muUpdateNoSpikes(mu0, prop_var, Counts, invdelta, nu, sum_bycell_all, mu_mu, s2_mu, q0, n, mu1, u, ind, Constrain, RefGene, ConstrainGene, NotConstrainGene, exponent, mintol)); + return rcpp_result_gen; +END_RCPP +} +// nuUpdateBatchNoSpikes +arma::mat nuUpdateBatchNoSpikes(arma::vec const& nu0, arma::vec const& prop_var, arma::mat const& Counts, arma::mat const& BatchDesign, arma::vec const& mu, arma::vec const& invdelta, arma::vec const& s, arma::vec const& thetaBatch, arma::vec const& sum_bygene_all, int const& q0, int const& n, arma::vec& nu1, arma::vec& u, arma::vec& ind, double const& exponent, double const& mintol); +RcppExport SEXP _BASiCS_nuUpdateBatchNoSpikes(SEXP nu0SEXP, SEXP prop_varSEXP, SEXP CountsSEXP, SEXP BatchDesignSEXP, SEXP muSEXP, SEXP invdeltaSEXP, SEXP sSEXP, SEXP thetaBatchSEXP, SEXP sum_bygene_allSEXP, SEXP q0SEXP, SEXP nSEXP, SEXP nu1SEXP, SEXP uSEXP, SEXP indSEXP, SEXP exponentSEXP, SEXP mintolSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::vec const& >::type nu0(nu0SEXP); + Rcpp::traits::input_parameter< arma::vec const& >::type prop_var(prop_varSEXP); + Rcpp::traits::input_parameter< arma::mat const& >::type Counts(CountsSEXP); + Rcpp::traits::input_parameter< arma::mat const& >::type BatchDesign(BatchDesignSEXP); + Rcpp::traits::input_parameter< arma::vec const& >::type mu(muSEXP); + Rcpp::traits::input_parameter< arma::vec const& >::type invdelta(invdeltaSEXP); + Rcpp::traits::input_parameter< arma::vec const& >::type s(sSEXP); + Rcpp::traits::input_parameter< arma::vec const& >::type thetaBatch(thetaBatchSEXP); + Rcpp::traits::input_parameter< arma::vec const& >::type sum_bygene_all(sum_bygene_allSEXP); + Rcpp::traits::input_parameter< int const& >::type q0(q0SEXP); + Rcpp::traits::input_parameter< int const& >::type n(nSEXP); + Rcpp::traits::input_parameter< arma::vec& >::type nu1(nu1SEXP); + Rcpp::traits::input_parameter< arma::vec& >::type u(uSEXP); + Rcpp::traits::input_parameter< arma::vec& >::type ind(indSEXP); + Rcpp::traits::input_parameter< double const& >::type exponent(exponentSEXP); + Rcpp::traits::input_parameter< double const& >::type mintol(mintolSEXP); + rcpp_result_gen = Rcpp::wrap(nuUpdateBatchNoSpikes(nu0, prop_var, Counts, BatchDesign, mu, invdelta, s, thetaBatch, sum_bygene_all, q0, n, nu1, u, ind, exponent, mintol)); + return rcpp_result_gen; +END_RCPP +} // designMatrix -arma::mat designMatrix(int const& k, arma::vec RBFLocations, arma::vec const& mu, double const& variance); +arma::mat designMatrix(int const& k, /* Number of Gaussian radial basis functions to use for regression */ arma::vec RBFLocations, arma::vec const& mu, double const& variance); RcppExport SEXP _BASiCS_designMatrix(SEXP kSEXP, SEXP RBFLocationsSEXP, SEXP muSEXP, SEXP varianceSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< int const& >::type k(kSEXP); - Rcpp::traits::input_parameter< arma::vec >::type RBFLocations(RBFLocationsSEXP); + Rcpp::traits::input_parameter< /* Number of Gaussian radial basis functions to use for regression */ arma::vec >::type RBFLocations(RBFLocationsSEXP); Rcpp::traits::input_parameter< arma::vec const& >::type mu(muSEXP); Rcpp::traits::input_parameter< double const& >::type variance(varianceSEXP); rcpp_result_gen = Rcpp::wrap(designMatrix(k, RBFLocations, mu, variance)); return rcpp_result_gen; END_RCPP } +// estimateRBFLocations +arma::vec estimateRBFLocations(arma::vec const& log_mu, int const& k, bool RBFMinMax); +RcppExport SEXP _BASiCS_estimateRBFLocations(SEXP log_muSEXP, SEXP kSEXP, SEXP RBFMinMaxSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::vec const& >::type log_mu(log_muSEXP); + Rcpp::traits::input_parameter< int const& >::type k(kSEXP); + Rcpp::traits::input_parameter< bool >::type RBFMinMax(RBFMinMaxSEXP); + rcpp_result_gen = Rcpp::wrap(estimateRBFLocations(log_mu, k, RBFMinMax)); + return rcpp_result_gen; +END_RCPP +} // muUpdateReg arma::mat muUpdateReg(arma::vec const& mu0, arma::vec const& prop_var, arma::mat const& Counts, arma::vec const& delta, arma::vec const& phinu, arma::vec const& sum_bycell_bio, arma::vec const& mu_mu, double const& s2_mu, int const& q0, int const& n, arma::vec& mu1, arma::vec& u, arma::vec& ind, int const& k, arma::vec const& lambda, arma::vec const& beta, arma::mat const& X, double const& sigma2, double variance, bool FixLocations, bool RBFMinMax, arma::vec RBFLocations, double const& exponent, double const& mintol); RcppExport SEXP _BASiCS_muUpdateReg(SEXP mu0SEXP, SEXP prop_varSEXP, SEXP CountsSEXP, SEXP deltaSEXP, SEXP phinuSEXP, SEXP sum_bycell_bioSEXP, SEXP mu_muSEXP, SEXP s2_muSEXP, SEXP q0SEXP, SEXP nSEXP, SEXP mu1SEXP, SEXP uSEXP, SEXP indSEXP, SEXP kSEXP, SEXP lambdaSEXP, SEXP betaSEXP, SEXP XSEXP, SEXP sigma2SEXP, SEXP varianceSEXP, SEXP FixLocationsSEXP, SEXP RBFMinMaxSEXP, SEXP RBFLocationsSEXP, SEXP exponentSEXP, SEXP mintolSEXP) { @@ -540,63 +588,8 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// muUpdateNoSpikes -arma::mat muUpdateNoSpikes(arma::vec const& mu0, arma::vec const& prop_var, arma::mat const& Counts, arma::vec const& invdelta, arma::vec const& nu, arma::vec const& sum_bycell_all, arma::vec const& mu_mu, double const& s2_mu, int const& q0, int const& n, arma::vec& mu1, arma::vec& u, arma::vec& ind, double const& Constrain, int const& RefGene, arma::uvec const& ConstrainGene, arma::uvec const& NotConstrainGene, double const& exponent, double const& mintol); -RcppExport SEXP _BASiCS_muUpdateNoSpikes(SEXP mu0SEXP, SEXP prop_varSEXP, SEXP CountsSEXP, SEXP invdeltaSEXP, SEXP nuSEXP, SEXP sum_bycell_allSEXP, SEXP mu_muSEXP, SEXP s2_muSEXP, SEXP q0SEXP, SEXP nSEXP, SEXP mu1SEXP, SEXP uSEXP, SEXP indSEXP, SEXP ConstrainSEXP, SEXP RefGeneSEXP, SEXP ConstrainGeneSEXP, SEXP NotConstrainGeneSEXP, SEXP exponentSEXP, SEXP mintolSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< arma::vec const& >::type mu0(mu0SEXP); - Rcpp::traits::input_parameter< arma::vec const& >::type prop_var(prop_varSEXP); - Rcpp::traits::input_parameter< arma::mat const& >::type Counts(CountsSEXP); - Rcpp::traits::input_parameter< arma::vec const& >::type invdelta(invdeltaSEXP); - Rcpp::traits::input_parameter< arma::vec const& >::type nu(nuSEXP); - Rcpp::traits::input_parameter< arma::vec const& >::type sum_bycell_all(sum_bycell_allSEXP); - Rcpp::traits::input_parameter< arma::vec const& >::type mu_mu(mu_muSEXP); - Rcpp::traits::input_parameter< double const& >::type s2_mu(s2_muSEXP); - Rcpp::traits::input_parameter< int const& >::type q0(q0SEXP); - Rcpp::traits::input_parameter< int const& >::type n(nSEXP); - Rcpp::traits::input_parameter< arma::vec& >::type mu1(mu1SEXP); - Rcpp::traits::input_parameter< arma::vec& >::type u(uSEXP); - Rcpp::traits::input_parameter< arma::vec& >::type ind(indSEXP); - Rcpp::traits::input_parameter< double const& >::type Constrain(ConstrainSEXP); - Rcpp::traits::input_parameter< int const& >::type RefGene(RefGeneSEXP); - Rcpp::traits::input_parameter< arma::uvec const& >::type ConstrainGene(ConstrainGeneSEXP); - Rcpp::traits::input_parameter< arma::uvec const& >::type NotConstrainGene(NotConstrainGeneSEXP); - Rcpp::traits::input_parameter< double const& >::type exponent(exponentSEXP); - Rcpp::traits::input_parameter< double const& >::type mintol(mintolSEXP); - rcpp_result_gen = Rcpp::wrap(muUpdateNoSpikes(mu0, prop_var, Counts, invdelta, nu, sum_bycell_all, mu_mu, s2_mu, q0, n, mu1, u, ind, Constrain, RefGene, ConstrainGene, NotConstrainGene, exponent, mintol)); - return rcpp_result_gen; -END_RCPP -} -// nuUpdateBatchNoSpikes -arma::mat nuUpdateBatchNoSpikes(arma::vec const& nu0, arma::vec const& prop_var, arma::mat const& Counts, arma::mat const& BatchDesign, arma::vec const& mu, arma::vec const& invdelta, arma::vec const& s, arma::vec const& thetaBatch, arma::vec const& sum_bygene_all, int const& q0, int const& n, arma::vec& nu1, arma::vec& u, arma::vec& ind, double const& exponent, double const& mintol); -RcppExport SEXP _BASiCS_nuUpdateBatchNoSpikes(SEXP nu0SEXP, SEXP prop_varSEXP, SEXP CountsSEXP, SEXP BatchDesignSEXP, SEXP muSEXP, SEXP invdeltaSEXP, SEXP sSEXP, SEXP thetaBatchSEXP, SEXP sum_bygene_allSEXP, SEXP q0SEXP, SEXP nSEXP, SEXP nu1SEXP, SEXP uSEXP, SEXP indSEXP, SEXP exponentSEXP, SEXP mintolSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< arma::vec const& >::type nu0(nu0SEXP); - Rcpp::traits::input_parameter< arma::vec const& >::type prop_var(prop_varSEXP); - Rcpp::traits::input_parameter< arma::mat const& >::type Counts(CountsSEXP); - Rcpp::traits::input_parameter< arma::mat const& >::type BatchDesign(BatchDesignSEXP); - Rcpp::traits::input_parameter< arma::vec const& >::type mu(muSEXP); - Rcpp::traits::input_parameter< arma::vec const& >::type invdelta(invdeltaSEXP); - Rcpp::traits::input_parameter< arma::vec const& >::type s(sSEXP); - Rcpp::traits::input_parameter< arma::vec const& >::type thetaBatch(thetaBatchSEXP); - Rcpp::traits::input_parameter< arma::vec const& >::type sum_bygene_all(sum_bygene_allSEXP); - Rcpp::traits::input_parameter< int const& >::type q0(q0SEXP); - Rcpp::traits::input_parameter< int const& >::type n(nSEXP); - Rcpp::traits::input_parameter< arma::vec& >::type nu1(nu1SEXP); - Rcpp::traits::input_parameter< arma::vec& >::type u(uSEXP); - Rcpp::traits::input_parameter< arma::vec& >::type ind(indSEXP); - Rcpp::traits::input_parameter< double const& >::type exponent(exponentSEXP); - Rcpp::traits::input_parameter< double const& >::type mintol(mintolSEXP); - rcpp_result_gen = Rcpp::wrap(nuUpdateBatchNoSpikes(nu0, prop_var, Counts, BatchDesign, mu, invdelta, s, thetaBatch, sum_bygene_all, q0, n, nu1, u, ind, exponent, mintol)); - return rcpp_result_gen; -END_RCPP -} // muUpdateRegNoSpikes -arma::mat muUpdateRegNoSpikes(arma::vec const& mu0, arma::vec const& prop_var, arma::mat const& Counts, arma::vec const& delta, arma::vec const& invdelta, arma::vec const& nu, arma::vec const& sum_bycell_all, arma::vec const& mu_mu, double const& s2_mu, int const& q0, int const& n, arma::vec& mu1, arma::vec& u, arma::vec& ind, double const& Constrain, int const& RefGene, arma::uvec const& ConstrainGene, arma::uvec const& NotConstrainGene, int const& k, arma::vec const& lambda, arma::vec const& beta, arma::mat const& X, double const& sigma2, double variance, bool FixLocations, bool RBFMinMax, arma::vec RBFLocations, double const& exponent, double const& mintol); +arma::mat muUpdateRegNoSpikes(arma::vec const& mu0, arma::vec const& prop_var, arma::mat const& Counts, arma::vec const& delta, arma::vec const& invdelta, arma::vec const& nu, arma::vec const& sum_bycell_all, arma::vec const& mu_mu, double const& s2_mu, int const& q0, int const& n, arma::vec& mu1, arma::vec& u, arma::vec& ind, double const& Constrain, /* No-spikes arguments from here */ int const& RefGene, arma::uvec const& ConstrainGene, arma::uvec const& NotConstrainGene, int const& k, /* Regression arguments from here */ arma::vec const& lambda, arma::vec const& beta, arma::mat const& X, double const& sigma2, double variance, bool FixLocations, bool RBFMinMax, arma::vec RBFLocations, double const& exponent, double const& mintol); RcppExport SEXP _BASiCS_muUpdateRegNoSpikes(SEXP mu0SEXP, SEXP prop_varSEXP, SEXP CountsSEXP, SEXP deltaSEXP, SEXP invdeltaSEXP, SEXP nuSEXP, SEXP sum_bycell_allSEXP, SEXP mu_muSEXP, SEXP s2_muSEXP, SEXP q0SEXP, SEXP nSEXP, SEXP mu1SEXP, SEXP uSEXP, SEXP indSEXP, SEXP ConstrainSEXP, SEXP RefGeneSEXP, SEXP ConstrainGeneSEXP, SEXP NotConstrainGeneSEXP, SEXP kSEXP, SEXP lambdaSEXP, SEXP betaSEXP, SEXP XSEXP, SEXP sigma2SEXP, SEXP varianceSEXP, SEXP FixLocationsSEXP, SEXP RBFMinMaxSEXP, SEXP RBFLocationsSEXP, SEXP exponentSEXP, SEXP mintolSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; @@ -616,11 +609,11 @@ BEGIN_RCPP Rcpp::traits::input_parameter< arma::vec& >::type u(uSEXP); Rcpp::traits::input_parameter< arma::vec& >::type ind(indSEXP); Rcpp::traits::input_parameter< double const& >::type Constrain(ConstrainSEXP); - Rcpp::traits::input_parameter< int const& >::type RefGene(RefGeneSEXP); + Rcpp::traits::input_parameter< /* No-spikes arguments from here */ int const& >::type RefGene(RefGeneSEXP); Rcpp::traits::input_parameter< arma::uvec const& >::type ConstrainGene(ConstrainGeneSEXP); Rcpp::traits::input_parameter< arma::uvec const& >::type NotConstrainGene(NotConstrainGeneSEXP); Rcpp::traits::input_parameter< int const& >::type k(kSEXP); - Rcpp::traits::input_parameter< arma::vec const& >::type lambda(lambdaSEXP); + Rcpp::traits::input_parameter< /* Regression arguments from here */ arma::vec const& >::type lambda(lambdaSEXP); Rcpp::traits::input_parameter< arma::vec const& >::type beta(betaSEXP); Rcpp::traits::input_parameter< arma::mat const& >::type X(XSEXP); Rcpp::traits::input_parameter< double const& >::type sigma2(sigma2SEXP); @@ -660,31 +653,42 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// rDirichlet +arma::vec rDirichlet(arma::vec alpha); +RcppExport SEXP _BASiCS_rDirichlet(SEXP alphaSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::vec >::type alpha(alphaSEXP); + rcpp_result_gen = Rcpp::wrap(rDirichlet(alpha)); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { - {"_BASiCS_BASiCS_MCMCcpp", (DL_FUNC) &_BASiCS_BASiCS_MCMCcpp, 41}, - {"_BASiCS_BASiCS_MCMCcppReg", (DL_FUNC) &_BASiCS_BASiCS_MCMCcppReg, 50}, - {"_BASiCS_BASiCS_MCMCcppNoSpikes", (DL_FUNC) &_BASiCS_BASiCS_MCMCcppNoSpikes, 43}, - {"_BASiCS_BASiCS_MCMCcppRegNoSpikes", (DL_FUNC) &_BASiCS_BASiCS_MCMCcppRegNoSpikes, 52}, {"_BASiCS_BASiCS_DenoisedRates", (DL_FUNC) &_BASiCS_BASiCS_DenoisedRates, 7}, - {"_BASiCS_estimateRBFLocations", (DL_FUNC) &_BASiCS_estimateRBFLocations, 3}, - {"_BASiCS_rDirichlet", (DL_FUNC) &_BASiCS_rDirichlet, 1}, + {"_BASiCS_BASiCS_MCMCcpp", (DL_FUNC) &_BASiCS_BASiCS_MCMCcpp, 42}, + {"_BASiCS_BASiCS_MCMCcppNoSpikes", (DL_FUNC) &_BASiCS_BASiCS_MCMCcppNoSpikes, 44}, + {"_BASiCS_BASiCS_MCMCcppReg", (DL_FUNC) &_BASiCS_BASiCS_MCMCcppReg, 51}, + {"_BASiCS_BASiCS_MCMCcppRegNoSpikes", (DL_FUNC) &_BASiCS_BASiCS_MCMCcppRegNoSpikes, 53}, {"_BASiCS_muUpdate", (DL_FUNC) &_BASiCS_muUpdate, 15}, {"_BASiCS_deltaUpdate", (DL_FUNC) &_BASiCS_deltaUpdate, 16}, {"_BASiCS_phiUpdate", (DL_FUNC) &_BASiCS_phiUpdate, 12}, {"_BASiCS_sUpdateBatch", (DL_FUNC) &_BASiCS_sUpdateBatch, 9}, {"_BASiCS_nuUpdateBatch", (DL_FUNC) &_BASiCS_nuUpdateBatch, 18}, {"_BASiCS_thetaUpdateBatch", (DL_FUNC) &_BASiCS_thetaUpdateBatch, 12}, + {"_BASiCS_muUpdateNoSpikes", (DL_FUNC) &_BASiCS_muUpdateNoSpikes, 19}, + {"_BASiCS_nuUpdateBatchNoSpikes", (DL_FUNC) &_BASiCS_nuUpdateBatchNoSpikes, 16}, {"_BASiCS_designMatrix", (DL_FUNC) &_BASiCS_designMatrix, 4}, + {"_BASiCS_estimateRBFLocations", (DL_FUNC) &_BASiCS_estimateRBFLocations, 3}, {"_BASiCS_muUpdateReg", (DL_FUNC) &_BASiCS_muUpdateReg, 24}, {"_BASiCS_deltaUpdateReg", (DL_FUNC) &_BASiCS_deltaUpdateReg, 16}, {"_BASiCS_betaUpdateReg", (DL_FUNC) &_BASiCS_betaUpdateReg, 3}, {"_BASiCS_sigma2UpdateReg", (DL_FUNC) &_BASiCS_sigma2UpdateReg, 10}, {"_BASiCS_lambdaUpdateReg", (DL_FUNC) &_BASiCS_lambdaUpdateReg, 8}, - {"_BASiCS_muUpdateNoSpikes", (DL_FUNC) &_BASiCS_muUpdateNoSpikes, 19}, - {"_BASiCS_nuUpdateBatchNoSpikes", (DL_FUNC) &_BASiCS_nuUpdateBatchNoSpikes, 16}, {"_BASiCS_muUpdateRegNoSpikes", (DL_FUNC) &_BASiCS_muUpdateRegNoSpikes, 29}, {"_BASiCS_deltaUpdateRegNoSpikes", (DL_FUNC) &_BASiCS_deltaUpdateRegNoSpikes, 16}, + {"_BASiCS_rDirichlet", (DL_FUNC) &_BASiCS_rDirichlet, 1}, {NULL, NULL, 0} }; diff --git a/src/general_utils.cpp b/src/general_utils.cpp deleted file mode 100644 index e168a2f7..00000000 --- a/src/general_utils.cpp +++ /dev/null @@ -1,455 +0,0 @@ -#include "utils.h" - -/* -* Rgig is an adaptation of the rgig.c function implemented by -* Ester Pantaleo and Robert B. Gramacy, 2010 -* (which in turn, was adapted from the C code in the ghyp -* v_1.5.2 package for R, rgig function) -* Our adaptation makes the function compatible with the Rcpp classes -*/ - -/* Auxiliary function taken from rgig.c */ -double gig_y_gfn(double y, double m, double beta, double lambda) -{ - /* - * gig_y_gfn: - * - * evaluate the function that we need to find the root - * of in order to construct the optimal rejection - * sampler to obtain GIG samples - */ - double y2, g; - y2 = y * y; - g = 0.5 * beta * y2 * y; - g -= y2 * (0.5 * beta * m + lambda + 1.0); - g += y * ((lambda - 1.0) * m - 0.5 * beta) + 0.5 * beta * m; - return(g); -} - -/* Auxiliary function taken from rgig.c */ -double rinvgauss(const double mu, const double lambda) -{ - /* - * rinvgauss: - * - * Michael/Schucany/Haas Method for generating Inverse Gaussian - * random variable with mean mu and shape parameter lambda, - * as given in Gentle's book on page 193 - */ - double u, y, x1, mu2, l2; - - y = norm_rand(); /// ****** - y *= y; - mu2 = mu * mu; - l2 = 2.0*lambda; - x1 = mu + mu2*y/l2 - (mu/l2)* sqrt(4.0*mu*lambda*y + mu2*y*y); - - u = unif_rand(); // ****** - if(u <= mu/(mu + x1)) return x1; - else return mu2/x1; -} - -/* Auxiliary function taken from rgig.c */ -double zeroin_gig(double ax,double bx, - double f(double x, double m, double beta, double lambda), - double tol,double m,double beta,double lambda) - /* An estimate to the root */ - //double ax; /* Left border | of the range */ - //double bx; /* Right border| the root is seeked*/ - /* Function under investigation */ - //double (*f)(double x, double m, double beta, double lambda); - //double tol; /* Acceptable tolerance */ - //double m; /* specific to gig_y_gfn */ - //double beta; /* specific to gig_y_gfn */ - //double lambda; /* specific to gig_y_gfn */ -{ - - double a,b,c; /* Abscissae, descr. see above */ - double fa; /* f(a) */ - double fb; /* f(b) */ - double fc; /* f(c) */ - - a = ax; b = bx; fa = (f)(a, m, beta, lambda); fb = (f)(b, m, beta, lambda); - c = a; fc = fa; - - for(;;) /* Main iteration loop */ - { - double prev_step = b-a; /* Distance from the last but one*/ - /* to the last approximation */ - double tol_act; /* Actual tolerance */ - double p; /* Interpolation step is calcu- */ - double q; /* lated in the form p/q; divi- */ - /* sion operations is delayed */ - /* until the last moment */ - double new_step; /* Step at this iteration */ - - if( fabs(fc) < fabs(fb) ) - { /* Swap data for b to be the */ - a = b; b = c; c = a; /* best approximation */ - fa=fb; fb=fc; fc=fa; - } - tol_act = 2.0*DOUBLE_EPS*fabs(b) + tol/2.0; - new_step = (c-b)/2.0; - - if( fabs(new_step) <= tol_act || fb == (double)0 ) - return b; /* Acceptable approx. is found */ - - /* Decide if the interpolation can be tried */ - if( fabs(prev_step) >= tol_act /* If prev_step was large enough*/ - && fabs(fa) > fabs(fb) ) /* and was in true direction, */ - { /* Interpolatiom may be tried */ - //register - double t1,cb,t2; - cb = c-b; - if( a==c ) /* If we have only two distinct */ - { /* points linear interpolation */ - t1 = fb/fa; /* can only be applied */ - p = cb*t1; - q = 1.0 - t1; - } - else /* Quadric inverse interpolation*/ - { - q = fa/fc; t1 = fb/fc; t2 = fb/fa; - p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) ); - q = (q-1.0) * (t1-1.0) * (t2-1.0); - } - - if( p>(double)0 ) /* p was calculated with the op-*/ - q = -q; /* posite sign; make p positive */ - else /* and assign possible minus to */ - p = -p; /* q */ - - if( p < (0.75*cb*q-fabs(tol_act*q)/2.0) /* If b+p/q falls in [b,c]*/ - && p < fabs(prev_step*q/2.0) ) /* and isn't too large */ - new_step = p/q; /* it is accepted */ - /* If p/q is too large then the */ - /* bissection procedure can */ - /* reduce [b,c] range to more */ - /* extent */ - } - - if( fabs(new_step) < tol_act ) - { /* Adjust the step to be not less*/ - if( new_step > (double)0 ) /* than tolerance */ - new_step = tol_act; - else - new_step = -tol_act; - } - - a = b; fa = fb; /* Save the previous approx. */ - b += new_step; fb = (*f)(b, m, beta, lambda); /* Do step to a new approxim. */ - if( (fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0) ) - { /* Adjust c for it to have a sign*/ - c = a; fc = fa; /* opposite to that of b */ - } - } -} - -/* Random sample generator from a Generalized Inverse Gaussian distribution -* (modified version of rgig.c using Rcpp classes) -*/ -NumericVector Rgig(const int n, - const double lambda, - const double chi, - const double psi) -{ - NumericVector samps(n); - /* special case which is basically a gamma distribution */ - if ((chi < ZTOL) & (lambda > 0.0)) { - int i; - for (i = 0; i < n; i++) { - samps(i) = R::rgamma(lambda, 2.0/psi); - } - return samps; - } - - /* special cases which is basically an inverse gamma distribution */ - if ((psi < ZTOL) & (lambda < 0.0)) { - int i; - for (i = 0; i < n; i++) { - samps(i) = 1.0/R::rgamma(0.0-lambda, 2.0/chi); - } - return samps; - } - - /* special case which is basically an inverse gaussian distribution */ - if (lambda == -0.5) { - double alpha; - int i; - alpha = sqrt(chi/psi); - for (i = 0; i < n; i++) { - samps(i) = rinvgauss(alpha, chi); - } - return samps; - } - - /* - * begin general purpose rgig code, which was basically - * translated from the R function rgig in the ghyp package v_1.5.2 - */ - - double alpha, beta, beta2, m, m1, lm1, lm12, upper, yM, yP, a, b, c, R1, R2, Y; - int i, need; - - alpha = sqrt(chi/psi); - beta2 = psi*chi; - beta = sqrt(psi*chi); - lm1 = lambda - 1.0; - lm12 = lm1*lm1; - m = (lm1 + sqrt(lm12 + beta2))/beta; - m1 = m + 1.0/m; - - upper = m; - while (gig_y_gfn(upper, m, beta, lambda) <= 0) { upper *= 2.0; } - - yM = zeroin_gig(0.0, m, gig_y_gfn, ZTOL, m, beta, lambda); - yP = zeroin_gig(m, upper, gig_y_gfn, ZTOL, m, beta, lambda); - - a = (yP - m) * pow(yP/m, 0.5 * lm1); - a *= exp(-0.25 * beta * (yP + 1.0/yP - m1)); - b = (yM - m) * pow(yM/m, 0.5 * lm1); - b *= exp(-0.25 * beta * (yM + 1/yM - m1)); - c = -0.25 * beta * m1 + 0.5 * lm1 * log(m); - - for (i=0; i 0.0) - { - if (-log(R1) >= - 0.5 * lm1 * log(Y) + 0.25 * beta * (Y + 1.0/Y) + c) - { - need = 0; - } - } - } - samps[i] = Y*alpha; - } - return(samps); -} - -/* Random sample generator from a Generalized Inverse Gaussian distribution -* (modified version of rgig.c using Rcpp classes) -*/ -// double RgigDouble(const double lambda, const double chi, const double psi) -// { -// double samps; -// /* special case which is basically a gamma distribution */ -// if((chi < ZTOL) & (lambda > 0.0)) { -// samps = R::rgamma(lambda, 2.0/psi); -// return samps; -// } - -// /* special cases which is basically an inverse gamma distribution */ -// if((psi < ZTOL) & (lambda < 0.0)) { -// samps = 1.0/R::rgamma(0.0-lambda, 2.0/chi); -// return samps; -// } - -// /* special case which is basically an inverse gaussian distribution */ -// if(lambda == -0.5) { -// double alpha; -// alpha = sqrt(chi/psi); -// samps = rinvgauss(alpha, chi); -// return samps; -// } - -// /* -// * begin general purpose rgig code, which was basically -// * translated from the R function rgig in the ghyp package v_1.5.2 -// */ - -// double alpha, beta, beta2, m, m1, lm1, lm12, upper, yM, yP, a, b, c, R1, R2, Y; -// int need; - -// alpha = sqrt(chi/psi); -// beta2 = psi*chi; -// beta = sqrt(psi*chi); -// lm1 = lambda - 1.0; -// lm12 = lm1*lm1; -// m = (lm1 + sqrt(lm12 + beta2))/beta; -// m1 = m + 1.0/m; - -// upper = m; -// while (gig_y_gfn(upper, m, beta, lambda) <= 0) { upper *= 2.0; } - -// yM = zeroin_gig(0.0, m, gig_y_gfn, ZTOL, m, beta, lambda); -// yP = zeroin_gig(m, upper, gig_y_gfn, ZTOL, m, beta, lambda); - -// a = (yP - m) * pow(yP/m, 0.5 * lm1); -// a *= exp(-0.25 * beta * (yP + 1.0/yP - m1)); -// b = (yM - m) * pow(yM/m, 0.5 * lm1); -// b *= exp(-0.25 * beta * (yM + 1/yM - m1)); -// c = -0.25 * beta * m1 + 0.5 * lm1 * log(m); - -// need = 1; -// while (need) -// { -// R1 = unif_rand(); -// R2 = unif_rand(); - -// Y = m + a * R2/R1 + b * (1.0 - R2)/R1; -// if (Y > 0.0) -// { -// if (-log(R1) >= - 0.5 * lm1 * log(Y) + 0.25 * beta * (Y + 1.0/Y) + c) -// { -// need = 0; -// } -// } -// } -// samps = Y*alpha; -// return(samps); -// } - - - -/* END OF ADAPTED CODE*/ - -void StartSampler(int const& N) -{ - Rcout << "-----------------------------------------------------" << std::endl; - Rcout << "MCMC sampler has been started: " << N << " iterations to go." << std::endl; - Rcout << "-----------------------------------------------------" << std::endl; -} - -void EndBurn() -{ - Rcout << "-----------------------------------------------------" << std::endl; - Rcout << "End of Burn-in period."<< std::endl; - Rcout << "-----------------------------------------------------" << std::endl; -} - -void CurrentIter(int const& i, int const& N) -{ - Rcout << "-----------------------------------------------------" << std::endl; - Rcout << "Iteration " << i << " out of " << N << " has been completed." << std::endl; - Rcout << "-----------------------------------------------------" << std::endl; - Rcout << "Current draws for selected parameters are displayed below." << std::endl; -} - -void EndSampler(int const& N) -{ - Rcout << " " << std::endl; - Rcout << "-----------------------------------------------------" << std::endl; - Rcout << "-----------------------------------------------------" << std::endl; - Rcout << "All " << N << " MCMC iterations have been completed." << std::endl; - Rcout << "-----------------------------------------------------" << std::endl; - Rcout << "-----------------------------------------------------" << std::endl; - Rcout << " " << std::endl; - Rcout << "-----------------------------------------------------" << std::endl; - Rcout << "Please see below a summary of the overall acceptance rates." << std::endl; - Rcout << "-----------------------------------------------------" << std::endl; -} - -void ReportAR(arma::vec const& AR, - std::string const& Param) -{ - Rcout << " " << std::endl; - Rcout << "Minimum acceptance rate among " << - Param << ": " << min(AR) << std::endl; - Rcout << "Average acceptance rate among " << - Param << ": " << mean(AR) << std::endl; - Rcout << "Maximum acceptance rate among " << - Param << ": " << max(AR) << std::endl; - Rcout << " " << std::endl; -} - -/* Auxiliary function that converts Rcpp::NumericVector - * objects into arma::vec objects - */ -arma::vec as_arma(NumericVector& x) -{ - return arma::vec(x.begin(), x.length(), false); -} - -/* Auxiliary function that converts Rcpp::NumericMatrix -* objects into arma::mac objects -*/ -arma::mat as_arma(NumericMatrix& x) -{ - return arma::mat(x.begin(), x.nrow(), x.ncol(), false); -} - -/* Auxiliary function: logarithm of a Gamma function applied -* element-wise to an arma::mat object -*/ -arma::mat lgamma_cpp(arma::mat const& x) -{ - arma::mat output = x; - for (unsigned int i=0; i max(x_arma) ) { offset = min(x_arma);} - else { offset = max(x_arma);} - return log(sum(exp(x_arma - offset))) + offset; -} - -arma::vec DegubInd(arma::vec ind, - int const q, - arma::vec const& u, - arma::vec const& log_aux, - arma::vec const& y, - double const& threshold, - std::string const& param) -{ - for (int i=0; i < q; i++) { - if( arma::is_finite(log_aux(i)) & arma::is_finite(y(i)) ) { - if((log(u(i)) < log_aux(i)) & (y(i) > threshold)) { ind(i) = 1; } - else{ind(i) = 0;} - } - else { - ind(i) = 0; - Rcpp::Rcout << "Error when updating " << param << " " << i+1 << std::endl; - Rcpp::Rcout << "Consider applying additional quality control" << std::endl; - Rcpp::Rcout << "to remove genes/cells when low total counts." << std::endl; - } - } - return ind; -} - -/*Dirichlet sampler*/ -arma::vec rDirichlet(arma::vec alpha) { - arma::vec aux = arma::ones(alpha.size()); - unsigned int i; - for(i=0; i +#include +#include +#include +#include +//File: matprod_arma.cpp +//[[Rcpp::depends(RcppArmadillo)]] +#include + +#ifdef _OPENMP + #include +#endif + +// [[Rcpp::plugins(openmp)]] + + +using namespace Rcpp; + +#endif diff --git a/src/utils_MCMCcpp.cpp b/src/updates.h similarity index 92% rename from src/utils_MCMCcpp.cpp rename to src/updates.h index 28ca5ba5..74764dba 100644 --- a/src/utils_MCMCcpp.cpp +++ b/src/updates.h @@ -1,8 +1,12 @@ +#ifndef UPDATES_H +#define UPDATES_H + #include "utils.h" /* Metropolis-Hastings updates of mu * Updates are implemented simulateaneously for all biological genes */ +// [[Rcpp::export(".muUpdate")]] arma::mat muUpdate( arma::vec const& mu0, arma::vec const& prop_var, @@ -33,6 +37,7 @@ arma::mat muUpdate( log_aux -= (0.5 / s2_mu) * (pow(log(mu1) - mu_mu, 2) - pow(log(mu0) - mu_mu, 2)) * exponent; + #pragma omp parallel for for (int i=0; i < q0; i++) { for (int j=0; j < n; j++) { log_aux(i) -= ( Counts(i,j) + invdelta(i) ) * @@ -59,6 +64,7 @@ arma::mat muUpdate( /* Metropolis-Hastings updates of delta * Updates are implemented simulateaneously for all biological genes. */ +// [[Rcpp::export(".deltaUpdate")]] arma::mat deltaUpdate( arma::vec const& delta0, arma::vec const& prop_var, @@ -88,6 +94,7 @@ arma::mat deltaUpdate( */ arma::vec log_aux = - n * (lgamma_cpp(1/delta1) - lgamma_cpp(1/delta0)); log_aux -= n * ( (log(delta1)/delta1) - (log(delta0)/delta0) ); + #pragma omp parallel for for (int i=0; i < q0; i++) { for (int j=0; j < n; j++) { log_aux(i) += std::lgamma(Counts(i,j) + (1/delta1(i))); @@ -131,6 +138,7 @@ arma::mat deltaUpdate( /* Metropolis-Hastings updates of phi * Joint updates using Dirichlet proposals */ +// [[Rcpp::export(".phiUpdate")]] Rcpp::List phiUpdate( arma::vec const& phi0, double const& prop_var, @@ -218,6 +226,7 @@ Rcpp::List phiUpdate( * have a closed form (Generalized Inverse Gaussian) * Updates are implemented simulateaneously for all cells. */ +// [[Rcpp::export(".sUpdateBatch")]] arma::vec sUpdateBatch( arma::vec const& s0, arma::vec const& nu, @@ -266,6 +275,7 @@ arma::vec sUpdateBatch( /* Metropolis-Hastings updates of nu (batch case) * Updates are implemented simulateaneously for all cells. */ +// [[Rcpp::export(".nuUpdateBatch")]] arma::mat nuUpdateBatch( arma::vec const& nu0, arma::vec const& prop_var, @@ -294,7 +304,7 @@ arma::mat nuUpdateBatch( // ACCEPT/REJECT STEP arma::vec log_aux = arma::zeros(n); - + #pragma omp parallel for for (int j=0; j < n; j++) { for (int i=0; i < q0; i++) { log_aux(j) -= (Counts(i, j) + invdelta(i)) * @@ -328,15 +338,22 @@ arma::mat nuUpdateBatch( /* Metropolis-Hastings updates of theta */ +/* theta0: Current value of $\theta$ */ +/* prop_var: Current value of the proposal variances for $\theta$ */ +/* s: Current value of $s=(s_1,...,s_n)$' */ +/* nu: Current value of $\nu=(\nu_1,...,\nu_n)'$ */ +/* a_theta: Shape hyper-parameter of the Gamma($a_{\theta}$,$b_{\theta}$) prior assigned to $\theta$ */ +/* b_theta: Rate hyper-parameter of the Gamma($a_{\theta}$,$b_{\theta}$) prior assigned to $\theta$ */ +// [[Rcpp::export(".thetaUpdateBatch")]] arma::mat thetaUpdateBatch( - arma::vec const& theta0, /* Current value of $\theta$ */ - arma::vec const& prop_var, /* Current value of the proposal variances for $\theta$ */ + arma::vec const& theta0, + arma::vec const& prop_var, arma::mat const& BatchDesign, arma::vec const& BatchSizes, - arma::vec const& s, /* Current value of $s=(s_1,...,s_n)$' */ - arma::vec const& nu, /* Current value of $\nu=(\nu_1,...,\nu_n)'$ */ - double const& a_theta, /* Shape hyper-parameter of the Gamma($a_{\theta}$,$b_{\theta}$) prior assigned to $\theta$ */ - double const& b_theta, /* Rate hyper-parameter of the Gamma($a_{\theta}$,$b_{\theta}$) prior assigned to $\theta$ */ + arma::vec const& s, + arma::vec const& nu, + double const& a_theta, + double const& b_theta, int const& n, int const& nBatch, double const& exponent, @@ -382,3 +399,5 @@ arma::mat thetaUpdateBatch( // OUTPUT return join_rows(theta, arma::conv_to::from(ind)); } + +#endif diff --git a/src/utils_MCMCcppNoSpikes.cpp b/src/updatesNoSpikes.h similarity index 95% rename from src/utils_MCMCcppNoSpikes.cpp rename to src/updatesNoSpikes.h index 00b8aa03..deae3c1f 100644 --- a/src/utils_MCMCcppNoSpikes.cpp +++ b/src/updatesNoSpikes.h @@ -1,8 +1,12 @@ +#ifndef UPDATESNOSPIKES_H +#define UPDATESNOSPIKES_H + #include "utils.h" /* Metropolis-Hastings updates of mu * Updates are implemented simulateaneously for all biological genes */ +// [[Rcpp::export(".muUpdateNoSpikes")]] arma::mat muUpdateNoSpikes( arma::vec const& mu0, arma::vec const& prop_var, @@ -42,6 +46,7 @@ arma::mat muUpdateNoSpikes( // Calculated in the same way for all genes, // but the reference one (no need to be sequential) arma::vec log_aux = (log(mu1) - log(mu0)) % sum_bycell_all; + #pragma omp parallel for for (int i = 0; i < q0; i++) { if (i != RefGene) { for (int j = 0; j < n; j++) { @@ -84,6 +89,7 @@ arma::mat muUpdateNoSpikes( // Step 2.3: For genes that are *not* under the constrain // Only relevant for a trimmed constrain + #pragma omp parallel for for (int i=0; i < nNotConstrainGene; i++) { iAux = NotConstrainGene(i); log_aux(iAux) -= (0.5 / s2_mu) * @@ -108,6 +114,7 @@ arma::mat muUpdateNoSpikes( /* Metropolis-Hastings updates of nu (batch case) * Updates are implemented simulateaneously for all cells. */ +// [[Rcpp::export(".nuUpdateBatchNoSpikes")]] arma::mat nuUpdateBatchNoSpikes( arma::vec const& nu0, arma::vec const& prop_var, @@ -137,6 +144,7 @@ arma::mat nuUpdateBatchNoSpikes( ((sum_bygene_all + 1 / thetaBatch) * exponent); log_aux -= (nu1 - nu0) % (1 / (thetaBatch % s * exponent)); + #pragma omp parallel for for (int j = 0; j < n; j++) { for (int i = 0; i < q0; i++) { log_aux(j) -= (Counts(i,j) + invdelta(i)) * @@ -162,3 +170,5 @@ arma::mat nuUpdateBatchNoSpikes( // OUTPUT return join_rows(nu1, ind); } + +#endif diff --git a/src/utils_MCMCcppReg.cpp b/src/updatesReg.h old mode 100755 new mode 100644 similarity index 95% rename from src/utils_MCMCcppReg.cpp rename to src/updatesReg.h index c68d2617..a25199ac --- a/src/utils_MCMCcppReg.cpp +++ b/src/updatesReg.h @@ -1,8 +1,12 @@ +#ifndef UPDATESREG_H +#define UPDATESREG_H + #include "utils.h" // Functions for regression case of BASiCS // Model matrix generation for regression +// [[Rcpp::export(".designMatrix")]] arma::mat designMatrix( int const& k, /* Number of Gaussian radial basis functions to use for regression */ arma::vec RBFLocations, @@ -23,6 +27,7 @@ arma::mat designMatrix( return X; } +// [[Rcpp::export(".estimateRBFLocations")]] arma::vec estimateRBFLocations( arma::vec const& log_mu, int const& k, @@ -51,6 +56,7 @@ arma::vec estimateRBFLocations( /* Metropolis-Hastings updates of mu * Updates are implemented simulateaneously for all biological genes */ +// [[Rcpp::export(".muUpdateReg")]] arma::mat muUpdateReg( arma::vec const& mu0, arma::vec const& prop_var, @@ -88,6 +94,7 @@ arma::mat muUpdateReg( arma::vec log_aux = (log(mu1) - log(mu0)) % sum_bycell_bio; log_aux -= (0.5 / s2_mu) * (pow(log(mu1) - mu_mu, 2) - pow(log(mu0) - mu_mu, 2)) * exponent; + #pragma omp parallel for for (int i = 0; i < q0; i++) { for (int j = 0; j < n; j++) { log_aux(i) -= ( Counts(i,j) + 1/delta(i) ) * @@ -150,6 +157,7 @@ arma::mat muUpdateReg( * sigma2: current value of sigma2 * beta: current value of beta */ +// [[Rcpp::export(".deltaUpdateReg")]] arma::mat deltaUpdateReg( arma::vec const& delta0, arma::vec const& prop_var, @@ -178,6 +186,7 @@ arma::mat deltaUpdateReg( */ arma::vec log_aux = - n * (lgamma_cpp(1/delta1) - lgamma_cpp(1/delta0)); log_aux -= n * ( (log(delta1)/delta1) - (log(delta0)/delta0) ); + #pragma omp parallel for for (int i=0; i < q0; i++) { for (int j=0; j < n; j++) { log_aux(i) += std::lgamma(Counts(i,j) + (1/delta1(i))); @@ -214,7 +223,7 @@ arma::mat deltaUpdateReg( // OUTPUT return join_rows(delta1, ind); } - +// [[Rcpp::export(".betaUpdateReg")]] arma::vec betaUpdateReg(double const& sigma2, arma::mat const& VAux, arma::vec const& mAux) { @@ -223,7 +232,7 @@ arma::vec betaUpdateReg(double const& sigma2, arma::vec beta = MVRNORM.row(0).t(); return beta; } - +// [[Rcpp::export(".sigma2UpdateReg")]] double sigma2UpdateReg(arma::vec const& delta, arma::vec const& beta, arma::vec const& lambda, @@ -243,7 +252,7 @@ double sigma2UpdateReg(arma::vec const& delta, double sigma2 = pow(R::rgamma(a, 1.0 / b), -1); return sigma2; } - +// [[Rcpp::export(".lambdaUpdateReg")]] arma::vec lambdaUpdateReg(arma::vec const& delta, arma::mat const& X, arma::vec const& beta, @@ -270,3 +279,5 @@ arma::vec lambdaUpdateReg(arma::vec const& delta, } return lambda1; } + +#endif diff --git a/src/utils_MCMCcppRegNoSpikes.cpp b/src/updatesRegNoSpikes.h old mode 100755 new mode 100644 similarity index 96% rename from src/utils_MCMCcppRegNoSpikes.cpp rename to src/updatesRegNoSpikes.h index bb05ada2..cbc0389a --- a/src/utils_MCMCcppRegNoSpikes.cpp +++ b/src/updatesRegNoSpikes.h @@ -1,8 +1,13 @@ +#ifndef UPDATESREGNOSPIKES_H +#define UPDATESREGNOSPIKES_H + + #include "utils.h" /* Metropolis-Hastings updates of mu * Updates are implemented simulateaneously for all biological genes */ +// [[Rcpp::export(".muUpdateRegNoSpikes")]] arma::mat muUpdateRegNoSpikes( arma::vec const& mu0, arma::vec const& prop_var, @@ -53,6 +58,7 @@ arma::mat muUpdateRegNoSpikes( // Calculated in the same way for all genes, // but the reference one (no need to be sequential) arma::vec log_aux = (log(mu1) - log(mu0)) % sum_bycell_all; + #pragma omp parallel for for (int i = 0; i < q0; i++) { if(i != RefGene) { for (int j=0; j < n; j++) { @@ -109,6 +115,7 @@ arma::mat muUpdateRegNoSpikes( // Step 2.3: For genes that are *not* under the constrain // Only relevant for a trimmed constrain + #pragma omp parallel for for (int i=0; i < nNotConstrainGene; i++) { iAux = NotConstrainGene(i); log_aux(iAux) -= (0.5 / s2_mu) * @@ -132,6 +139,7 @@ arma::mat muUpdateRegNoSpikes( /* Metropolis-Hastings updates of delta * Updates are implemented simulateaneously for all biological genes */ +// [[Rcpp::export(".deltaUpdateRegNoSpikes")]] arma::mat deltaUpdateRegNoSpikes( arma::vec const& delta0, arma::vec const& prop_var, @@ -163,6 +171,7 @@ arma::mat deltaUpdateRegNoSpikes( log_aux -= n * ((log(delta1) / delta1) - (log(delta0) / delta0)); // Loop to replace matrix operations, through genes and cells + #pragma omp parallel for for (int i = 0; i < q0; i++) { for (int j = 0; j < n; j++) { log_aux(i) += std::lgamma(Counts(i, j) + (1 / delta1(i))); @@ -194,5 +203,4 @@ arma::mat deltaUpdateRegNoSpikes( return join_rows(delta1, ind); } - - +#endif diff --git a/src/utils.h b/src/utils.h old mode 100755 new mode 100644 index 77ee40cf..3caa673f --- a/src/utils.h +++ b/src/utils.h @@ -1,338 +1,390 @@ -#include "MCMCcpp.h" - #ifndef UTILS_H #define UTILS_H -/* Declarations for general utility functions - * Stored in: general_utils.cpp - */ +#include "libraries.h" + #define ZTOL sqrt(DOUBLE_EPS) +/* +* Rgig is an adaptation of the rgig.c function implemented by +* Ester Pantaleo and Robert B. Gramacy, 2010 +* (which in turn, was adapted from the C code in the ghyp +* v_1.5.2 package for R, rgig function) +* Our adaptation makes the function compatible with the Rcpp classes +*/ + +/* Auxiliary function taken from rgig.c */ +double gig_y_gfn(double y, double m, double beta, double lambda) +{ + /* + * gig_y_gfn: + * + * evaluate the function that we need to find the root + * of in order to construct the optimal rejection + * sampler to obtain GIG samples + */ + double y2, g; + y2 = y * y; + g = 0.5 * beta * y2 * y; + g -= y2 * (0.5 * beta * m + lambda + 1.0); + g += y * ((lambda - 1.0) * m - 0.5 * beta) + 0.5 * beta * m; + return(g); +} + +/* Auxiliary function taken from rgig.c */ +double rinvgauss(const double mu, const double lambda) +{ + /* + * rinvgauss: + * + * Michael/Schucany/Haas Method for generating Inverse Gaussian + * random variable with mean mu and shape parameter lambda, + * as given in Gentle's book on page 193 + */ + double u, y, x1, mu2, l2; + + y = norm_rand(); /// ****** + y *= y; + mu2 = mu * mu; + l2 = 2.0*lambda; + x1 = mu + mu2*y/l2 - (mu/l2)* sqrt(4.0*mu*lambda*y + mu2*y*y); + + u = unif_rand(); // ****** + if(u <= mu/(mu + x1)) return x1; + else return mu2/x1; +} -double gig_y_gfn(double y, double m, double beta, double lambda); -double rinvgauss(const double mu, const double lambda); +/* Auxiliary function taken from rgig.c */ double zeroin_gig(double ax,double bx, double f(double x, double m, double beta, double lambda), - double tol,double m,double beta,double lambda); + double tol,double m,double beta,double lambda) + /* An estimate to the root */ + //double ax; /* Left border | of the range */ + //double bx; /* Right border| the root is seeked*/ + /* Function under investigation */ + //double (*f)(double x, double m, double beta, double lambda); + //double tol; /* Acceptable tolerance */ + //double m; /* specific to gig_y_gfn */ + //double beta; /* specific to gig_y_gfn */ + //double lambda; /* specific to gig_y_gfn */ +{ + + double a,b,c; /* Abscissae, descr. see above */ + double fa; /* f(a) */ + double fb; /* f(b) */ + double fc; /* f(c) */ + + a = ax; b = bx; fa = (f)(a, m, beta, lambda); fb = (f)(b, m, beta, lambda); + c = a; fc = fa; + + for(;;) /* Main iteration loop */ + { + double prev_step = b-a; /* Distance from the last but one*/ + /* to the last approximation */ + double tol_act; /* Actual tolerance */ + double p; /* Interpolation step is calcu- */ + double q; /* lated in the form p/q; divi- */ + /* sion operations is delayed */ + /* until the last moment */ + double new_step; /* Step at this iteration */ + + if( fabs(fc) < fabs(fb) ) + { /* Swap data for b to be the */ + a = b; b = c; c = a; /* best approximation */ + fa=fb; fb=fc; fc=fa; + } + tol_act = 2.0*DOUBLE_EPS*fabs(b) + tol/2.0; + new_step = (c-b)/2.0; + + if( fabs(new_step) <= tol_act || fb == (double)0 ) + return b; /* Acceptable approx. is found */ + + /* Decide if the interpolation can be tried */ + if( fabs(prev_step) >= tol_act /* If prev_step was large enough*/ + && fabs(fa) > fabs(fb) ) /* and was in true direction, */ + { /* Interpolatiom may be tried */ + //register + double t1,cb,t2; + cb = c-b; + if( a==c ) /* If we have only two distinct */ + { /* points linear interpolation */ + t1 = fb/fa; /* can only be applied */ + p = cb*t1; + q = 1.0 - t1; + } + else /* Quadric inverse interpolation*/ + { + q = fa/fc; t1 = fb/fc; t2 = fb/fa; + p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) ); + q = (q-1.0) * (t1-1.0) * (t2-1.0); + } + + if( p>(double)0 ) /* p was calculated with the op-*/ + q = -q; /* posite sign; make p positive */ + else /* and assign possible minus to */ + p = -p; /* q */ + + if( p < (0.75*cb*q-fabs(tol_act*q)/2.0) /* If b+p/q falls in [b,c]*/ + && p < fabs(prev_step*q/2.0) ) /* and isn't too large */ + new_step = p/q; /* it is accepted */ + /* If p/q is too large then the */ + /* bissection procedure can */ + /* reduce [b,c] range to more */ + /* extent */ + } + + if( fabs(new_step) < tol_act ) + { /* Adjust the step to be not less*/ + if( new_step > (double)0 ) /* than tolerance */ + new_step = tol_act; + else + new_step = -tol_act; + } + + a = b; fa = fb; /* Save the previous approx. */ + b += new_step; fb = (*f)(b, m, beta, lambda); /* Do step to a new approxim. */ + if( (fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0) ) + { /* Adjust c for it to have a sign*/ + c = a; fc = fa; /* opposite to that of b */ + } + } +} + +/* Random sample generator from a Generalized Inverse Gaussian distribution +* (modified version of rgig.c using Rcpp classes) +*/ NumericVector Rgig(const int n, const double lambda, const double chi, - const double psi); - -// double RgigDouble(const double lambda, const double chi, const double psi); - -void StartSampler(int const& N); -void EndBurn(); -void CurrentIter(int const& i, int const& N); -void EndSampler(int const& N); -void ReportAR(arma::vec const& AR, - std::string const& Param); -arma::vec as_arma(NumericVector& x); -arma::mat as_arma(NumericMatrix& x); -arma::mat lgamma_cpp(arma::mat const& x); -arma::vec lgamma_cpp_vec(arma::vec const& x); -double log_sum_exp_cpp(arma::vec const& x_arma); -arma::vec DegubInd(arma::vec ind, - int const q, - arma::vec const& u, - arma::vec const& log_aux, - arma::vec const& y, - double const& threshold, - std::string const& param); + const double psi) +{ + NumericVector samps(n); + /* special case which is basically a gamma distribution */ + if ((chi < ZTOL) & (lambda > 0.0)) { + int i; + for (i = 0; i < n; i++) { + samps(i) = R::rgamma(lambda, 2.0/psi); + } + return samps; + } + + /* special cases which is basically an inverse gamma distribution */ + if ((psi < ZTOL) & (lambda < 0.0)) { + int i; + for (i = 0; i < n; i++) { + samps(i) = 1.0/R::rgamma(0.0-lambda, 2.0/chi); + } + return samps; + } + + /* special case which is basically an inverse gaussian distribution */ + if (lambda == -0.5) { + double alpha; + int i; + alpha = sqrt(chi/psi); + for (i = 0; i < n; i++) { + samps(i) = rinvgauss(alpha, chi); + } + return samps; + } + + /* + * begin general purpose rgig code, which was basically + * translated from the R function rgig in the ghyp package v_1.5.2 + */ + + double alpha, beta, beta2, m, m1, lm1, lm12, upper, yM, yP, a, b, c, R1, R2, Y; + int i, need; + + alpha = sqrt(chi/psi); + beta2 = psi*chi; + beta = sqrt(psi*chi); + lm1 = lambda - 1.0; + lm12 = lm1*lm1; + m = (lm1 + sqrt(lm12 + beta2))/beta; + m1 = m + 1.0/m; + + upper = m; + while (gig_y_gfn(upper, m, beta, lambda) <= 0) { upper *= 2.0; } + + yM = zeroin_gig(0.0, m, gig_y_gfn, ZTOL, m, beta, lambda); + yP = zeroin_gig(m, upper, gig_y_gfn, ZTOL, m, beta, lambda); + + a = (yP - m) * pow(yP/m, 0.5 * lm1); + a *= exp(-0.25 * beta * (yP + 1.0/yP - m1)); + b = (yM - m) * pow(yM/m, 0.5 * lm1); + b *= exp(-0.25 * beta * (yM + 1/yM - m1)); + c = -0.25 * beta * m1 + 0.5 * lm1 * log(m); + + for (i=0; i 0.0) + { + if (-log(R1) >= - 0.5 * lm1 * log(Y) + 0.25 * beta * (Y + 1.0/Y) + c) + { + need = 0; + } + } + } + samps[i] = Y*alpha; + } + return(samps); +} -// [[Rcpp::export(".estimateRBFLocations")]] -arma::vec estimateRBFLocations( - arma::vec const& log_mu, - int const& k, - bool RBFMinMax); -// [[Rcpp::export(.rDirichlet)]] -arma::vec rDirichlet(arma::vec alpha); -arma::mat mvrnormArma(int n, arma::vec mu, arma::mat sigma); -/* Declarations for functions used by the main MCMC sampler - * Stored in: utils_MCMCcpp.cpp - */ -// [[Rcpp::export(".muUpdate")]] -arma::mat muUpdate( - arma::vec const& mu0, - arma::vec const& prop_var, - arma::mat const& Counts, - arma::vec const& invdelta, - arma::vec const& phinu, - arma::vec const& sum_bycell_bio, - arma::vec const& mu_mu, - double const& s2_mu, - int const& q0, - int const& n, - arma::vec & mu1, - arma::vec & u, - arma::vec & ind, - double const& exponent, - double const& mintol); -// [[Rcpp::export(".deltaUpdate")]] -arma::mat deltaUpdate( - arma::vec const& delta0, - arma::vec const& prop_var, - arma::mat const& Counts, - arma::vec const& mu, - arma::vec const& phinu, - double const& a_delta, - double const& b_delta, - double const& s2delta, - double const& prior_delta, - int const& q0, - int const& n, - arma::vec & delta1, - arma::vec & u, - arma::vec & ind, - double const& exponent, - double const& mintol); -// [[Rcpp::export(".phiUpdate")]] -Rcpp::List phiUpdate( - arma::vec const& phi0, - double const& prop_var, - arma::mat const& Counts, - arma::vec const& mu, - arma::vec const& invdelta, - arma::vec const& nu, - arma::vec const& aphi, - arma::vec const& sum_bygene_bio, - int const& q0, - int const& n, - arma::vec & phi1, - double const& exponent); +void StartSampler(int const& N) +{ + Rcout << "-----------------------------------------------------" << std::endl; + Rcout << "MCMC sampler has been started: " << N << " iterations to go." << std::endl; + Rcout << "-----------------------------------------------------" << std::endl; +} -// [[Rcpp::export(".sUpdateBatch")]] -arma::vec sUpdateBatch( - arma::vec const& s0, - arma::vec const& nu, - arma::vec const& thetaBatch, - double const& as, - double const& bs, - arma::mat const& BatchDesign, - int const& n, - arma::vec & s1, - double const& exponent); +void EndBurn() +{ + Rcout << "-----------------------------------------------------" << std::endl; + Rcout << "End of Burn-in period."<< std::endl; + Rcout << "-----------------------------------------------------" << std::endl; +} -// [[Rcpp::export(".nuUpdateBatch")]] -arma::mat nuUpdateBatch( - arma::vec const& nu0, - arma::vec const& prop_var, - arma::mat const& Counts, - double const& SumSpikeInput, - arma::mat const& BatchDesign, - arma::vec const& mu, - arma::vec const& invdelta, - arma::vec const& phi, - arma::vec const& s, - arma::vec const& thetaBatch, - arma::vec const& sum_bygene_all, - int const& q0, - int const& n, - arma::vec & nu1, - arma::vec & u, - arma::vec & ind, - double const& exponent, - double const& mintol); +void CurrentIter(int const& i, int const& N) +{ + Rcout << "-----------------------------------------------------" << std::endl; + Rcout << "Iteration " << i << " out of " << N << " has been completed." << std::endl; + Rcout << "-----------------------------------------------------" << std::endl; + Rcout << "Current draws for selected parameters are displayed below." << std::endl; +} -// [[Rcpp::export(".thetaUpdateBatch")]] -arma::mat thetaUpdateBatch( - arma::vec const& theta0, - arma::vec const& prop_var, - arma::mat const& BatchDesign, - arma::vec const& BatchSizes, - arma::vec const& s, - arma::vec const& nu, - double const& a_theta, - double const& b_theta, - int const& n, - int const& nBatch, - double const& exponent, - double const& mintol); - - /* Declarations for functions used by the MCMC sampler for the regression case - * Stored in: utils_MCMCcppReg.cpp - */ +void EndSampler(int const& N) +{ + Rcout << " " << std::endl; + Rcout << "-----------------------------------------------------" << std::endl; + Rcout << "-----------------------------------------------------" << std::endl; + Rcout << "All " << N << " MCMC iterations have been completed." << std::endl; + Rcout << "-----------------------------------------------------" << std::endl; + Rcout << "-----------------------------------------------------" << std::endl; + Rcout << " " << std::endl; + Rcout << "-----------------------------------------------------" << std::endl; + Rcout << "Please see below a summary of the overall acceptance rates." << std::endl; + Rcout << "-----------------------------------------------------" << std::endl; +} -// [[Rcpp::export(".designMatrix")]] -arma::mat designMatrix( - int const& k, - arma::vec RBFLocations, - arma::vec const& mu, - double const& variance); +void ReportAR(arma::vec const& AR, + std::string const& Param) +{ + Rcout << " " << std::endl; + Rcout << "Minimum acceptance rate among " << + Param << ": " << min(AR) << std::endl; + Rcout << "Average acceptance rate among " << + Param << ": " << mean(AR) << std::endl; + Rcout << "Maximum acceptance rate among " << + Param << ": " << max(AR) << std::endl; + Rcout << " " << std::endl; +} -// [[Rcpp::export(".muUpdateReg")]] -arma::mat muUpdateReg( - arma::vec const& mu0, - arma::vec const& prop_var, - arma::mat const& Counts, - arma::vec const& delta, - arma::vec const& phinu, - arma::vec const& sum_bycell_bio, - arma::vec const& mu_mu, - double const& s2_mu, - int const& q0, - int const& n, - arma::vec & mu1, - arma::vec & u, - arma::vec & ind, - int const& k, - arma::vec const& lambda, - arma::vec const& beta, - arma::mat const& X, - double const& sigma2, - double variance, - bool FixLocations, - bool RBFMinMax, - arma::vec RBFLocations, - double const& exponent, - double const& mintol); +/* Auxiliary function that converts Rcpp::NumericVector + * objects into arma::vec objects + */ +arma::vec as_arma(NumericVector& x) +{ + return arma::vec(x.begin(), x.length(), false); +} -// [[Rcpp::export(".deltaUpdateReg")]] -arma::mat deltaUpdateReg( - arma::vec const& delta0, - arma::vec const& prop_var, - arma::mat const& Counts, - arma::vec const& mu, - arma::vec const& phinu, - int const& q0, - int const& n, - arma::vec & delta1, - arma::vec & u, - arma::vec & ind, - arma::vec const& lambda, - arma::mat const& X, - double const& sigma2, - arma::vec const& beta, - double const& exponent, - double const& mintol); +/* Auxiliary function that converts Rcpp::NumericMatrix +* objects into arma::mac objects +*/ +arma::mat as_arma(NumericMatrix& x) +{ + return arma::mat(x.begin(), x.nrow(), x.ncol(), false); +} -// [[Rcpp::export(".betaUpdateReg")]] -arma::vec betaUpdateReg(double const& sigma2, - arma::mat const& VAux, - arma::vec const& mAux); +/* Auxiliary function: logarithm of a Gamma function applied +* element-wise to an arma::mat object +*/ +arma::mat lgamma_cpp(arma::mat const& x) +{ + arma::mat output = x; + for (unsigned int i=0; i max(x_arma) ) { offset = min(x_arma);} + else { offset = max(x_arma);} + return log(sum(exp(x_arma - offset))) + offset; +} -/* Declarations for functions used by the MCMC sampler for the non-spikes case - * Stored in: utils_MCMCcppNoSpikes.cpp - */ -// [[Rcpp::export(".muUpdateNoSpikes")]] -arma::mat muUpdateNoSpikes( - arma::vec const& mu0, - arma::vec const& prop_var, - arma::mat const& Counts, - arma::vec const& invdelta, - arma::vec const& nu, - arma::vec const& sum_bycell_all, - arma::vec const& mu_mu, - double const& s2_mu, - int const& q0, - int const& n, - arma::vec & mu1, - arma::vec & u, - arma::vec & ind, - double const& Constrain, - int const& RefGene, - arma::uvec const& ConstrainGene, - arma::uvec const& NotConstrainGene, - double const& exponent, - double const& mintol); +arma::vec DegubInd(arma::vec ind, + int const q, + arma::vec const& u, + arma::vec const& log_aux, + arma::vec const& y, + double const& threshold, + std::string const& param) +{ + for (int i=0; i < q; i++) { + if( arma::is_finite(log_aux(i)) & arma::is_finite(y(i)) ) { + if((log(u(i)) < log_aux(i)) & (y(i) > threshold)) { ind(i) = 1; } + else{ind(i) = 0;} + } + else { + ind(i) = 0; + Rcpp::Rcout << "Error when updating " << param << " " << i+1 << std::endl; + Rcpp::Rcout << "Consider applying additional quality control" << std::endl; + Rcpp::Rcout << "to remove genes/cells when low total counts." << std::endl; + } + } + return ind; +} -// [[Rcpp::export(".nuUpdateBatchNoSpikes")]] -arma::mat nuUpdateBatchNoSpikes( - arma::vec const& nu0, - arma::vec const& prop_var, - arma::mat const& Counts, - arma::mat const& BatchDesign, - arma::vec const& mu, - arma::vec const& invdelta, - arma::vec const& s, - arma::vec const& thetaBatch, - arma::vec const& sum_bygene_all, - int const& q0, - int const& n, - arma::vec & nu1, - arma::vec & u, - arma::vec & ind, - double const& exponent, - double const& mintol); +/*Dirichlet sampler*/ +// [[Rcpp::export(".rDirichlet")]] +arma::vec rDirichlet(arma::vec alpha) { + arma::vec aux = arma::ones(alpha.size()); + unsigned int i; + for(i=0; i