Add ability to set seed

R/RcppExports.R

@@ -37,6 +37,10 @@
     .Call(`_ToxicR_polyk`, dose, tumor, daysOnStudy)
 }
 
+setseedGSL <- function(s) {
+    invisible(.Call(`_ToxicR_setseedGSL`, s))
+}
+
 .set_threads <- function(num_threads) {
     invisible(.Call(`_ToxicR_set_threads`, num_threads))
 }

R/continuous_wrappers.R

@@ -73,6 +73,7 @@
 #' @param transform Transforms doses using \eqn{\log(dose+\sqrt{dose^2+1})}. Note: this is a log transform that has a derivative defined when dose =0.
 #' @param BMD_TYPE Deprecated version of BMR_TYPE that specifies the type of benchmark dose analysis to be performed
 #' @param threads specify the number of OpenMP threads to use for the calculations. Default = 2
+#' @param seed sets the GSL seed. Default = 12331
 #' @return Returns a model object class with the following structure:
 #' \itemize{
 #'    \item \code{full_model}:  The model along with the likelihood distribution. 
@@ -112,7 +113,8 @@ single_continuous_fit <- function(D,Y,model_type="hill", fit_type = "laplace",
                                    BMR = 0.1, point_p = 0.01, distribution = "normal-ncv",
                                    alpha = 0.05, samples = 25000, degree=2,
                                    burnin = 1000, BMD_priors = FALSE, ewald = FALSE,
-                                   transform = FALSE, BMD_TYPE = NA, threads = 2){
+                                   transform = FALSE, BMD_TYPE = NA, threads = 2, seed = 12331){
+    setseedGSL(seed)
     Y <- as.matrix(Y) 
     D <- as.matrix(D) 
 

R/dichotomous_wrappers.R

@@ -14,6 +14,7 @@
 #' @param samples the number of samples to take (MCMC only)
 #' @param burnin the number of burnin samples to take (MCMC only)
 #' @param threads specify the number of OpenMP threads to use for the calculations. Default = 2
+#' @param seed set the GSL seed. Default = 12331
 #'
 #' @return Returns a model object class with the following structure:
 #' \itemize{
@@ -55,7 +56,8 @@
 single_dichotomous_fit <- function(D, Y, N, model_type, fit_type = "laplace",
                                    prior = NULL, BMR = 0.1,
                                    alpha = 0.05, degree = 2, samples = 21000,
-                                   burnin = 1000, threads=2) {
+                                   burnin = 1000, threads=2, seed = 12331) {
+  setseedGSL(seed)
   Y <- as.matrix(Y)
   D <- as.matrix(D)
   N <- as.matrix(N)

src/RcppExports.cpp

@@ -149,6 +149,16 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// setseedGSL
+void setseedGSL(const int s);
+RcppExport SEXP _ToxicR_setseedGSL(SEXP sSEXP) {
+BEGIN_RCPP
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const int >::type s(sSEXP);
+    setseedGSL(s);
+    return R_NilValue;
+END_RCPP
+}
 // set_threads
 void set_threads(int num_threads);
 RcppExport SEXP _ToxicR_set_threads(SEXP num_threadsSEXP) {
@@ -170,6 +180,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_ToxicR_run_dichotomous_single_mcmc", (DL_FUNC) &_ToxicR_run_dichotomous_single_mcmc, 5},
     {"_ToxicR_run_continuous_single_mcmc", (DL_FUNC) &_ToxicR_run_continuous_single_mcmc, 7},
     {"_ToxicR_polyk", (DL_FUNC) &_ToxicR_polyk, 3},
+    {"_ToxicR_setseedGSL", (DL_FUNC) &_ToxicR_setseedGSL, 1},
     {"_ToxicR_set_threads", (DL_FUNC) &_ToxicR_set_threads, 1},
     {NULL, NULL, 0}
 };

src/include/seeder.h

@@ -0,0 +1,19 @@
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+
+class Seeder {
+    public:
+    Seeder(int seed=12345678) {
+        gsl_rng_env_setup();
+        r = gsl_rng_alloc(gsl_rng_default);
+        gsl_rng_set(r, seed);
+    }
+    ~Seeder() {
+        gsl_rng_free(r);
+    }
+    void setSeed(const int seed) {
+        gsl_rng_set(r, seed);
+    }
+    private:
+    gsl_rng *r;
+};

src/main_seed.cpp

@@ -0,0 +1,8 @@
+#include "seeder.h"
+static Seeder seeder;
+
+// [[Rcpp::export]]
+void setseedGSL(const int s) {
+    seeder.setSeed(s);
+    return;
+}