Improve C++ speed, remove batch requirement (#187)

version 2.1.8 (2020-09-25)

- 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.

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

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`. 
+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.

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
   )
 
 

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

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

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,25 +48,25 @@
   }
 
   # 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"
     )
   }
 
   # 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", 

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"

src/BASiCS_DenoisedRates.cpp → 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

src/BASiCS_MCMCcpp.cpp → 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

src/BASiCS_MCMCcppNoSpikes.cpp → 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

src/BASiCS_MCMCcppReg.cpp → 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

src/BASiCS_MCMCcppRegNoSpikes.cpp → 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