Add documentation for seeder function. Refactor seeder class for thre…

…ad safety and better functionality

.Rbuildignore

@@ -10,4 +10,5 @@
 \.so$
 \.git/
 ^\.dll$
-^\.github$
+^\.github$
+^\.log$

.gitignore

@@ -4,4 +4,5 @@
 .Rhistory
 .vscode/*
 .RData
-*.tar.gz
+*.tar.gz
+*.dll

DESCRIPTION

@@ -78,7 +78,7 @@ Imports: Rcpp (>= 1.0.0), ggplot2 (>= 3.3.2), shiny (>= 1.5.0), coda
         (>= 0.19-4), scales (>= 1.1.1), forcats, ggridges (>= 0.5.3),
         doBy (>= 4.6.11), multcomp (>= 1.4), dplyr (>= 1.0.7), remotes
 LinkingTo: Rcpp, RcppEigen, RcppGSL
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.0
 Encoding: UTF-8
 VignetteBuilder: knitr
 Suggests: plotly (>= 4.9.2.1), rmarkdown, actuar (>= 3.2-0),ggpubr (>=

R/RcppExports.R

@@ -37,6 +37,15 @@
     .Call(`_ToxicR_polyk`, dose, tumor, daysOnStudy)
 }
 
+#' Sets a GSL seed
+#' 
+#' @title setseedGSL - Sets a GSL seed
+#' @param s positive integer used as the seeds
+#' @return None
+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 specify 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. 
@@ -104,15 +105,16 @@
 #' M2[,3] <- c(20,20,19,20,20)
 #' M2[,4] <- c(1.2,1.1,0.81,0.74,0.66)
 #' model = single_continuous_fit(M2[,1,drop=FALSE], M2[,2:4], BMR_TYPE="sd", BMR=1, ewald = TRUE,
-#'                              distribution = "normal",fit_type="laplace",model_type = "hill",threads = 2)
+#'                              distribution = "normal",fit_type="laplace",model_type = "hill",threads = 2, seed = 12331)
 #' 
 #' summary(model)
 single_continuous_fit <- function(D,Y,model_type="hill", fit_type = "laplace",
                                    prior=NA, BMR_TYPE = "sd", 
                                    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 specify 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)

R/model_averaging_fits.R

@@ -415,6 +415,7 @@ ma_continuous_fit <- function(D, Y, model_list = NA, fit_type = "laplace",
 #' @param samples the number of samples to take (MCMC only)
 #' @param burnin the number of burnin samples to take (MCMC only)
 #' @param BMD_TYPE Deprecated version of BMR_TYPE that specifies the type of benchmark dose analysis to be performed
+#' @param threads number of threads to use. Default = 2
 #' @return a model object containing a list of single models
 #' \itemize{
 #'  \item \code{Individual_Model_X}: Here \code{X} is a number \eqn{1\leq X \leq n,} where \eqn{n}

src/Makevars.in

@@ -1,4 +1,4 @@
-PKG_CXXFLAGS=  -I./code_base -I./include @OPENMP@ @NLOPT_CPPFLAGS@ @GSL_CPPFLAGS@ -DR_COMPILATION -DEIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS
+PKG_CXXFLAGS=  -I./code_base -I./include @OPENMP@ @NLOPT_CPPFLAGS@ @GSL_CPPFLAGS@ -DR_COMPILATION -DEIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS -DEIGEN_WARNINGS_DISABLED
 PKG_LIBS= $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) @NLOPT_LIBS@  @GSL_LIBS@ @OPENMP@
 CXX14FLAGS= -Wno-ignored-attributes -Wno-sign-compare -Wno-unused -Wno-builtin-macro-redefined
 CXX11FLAGS= -Wno-ignored-attributes -Wno-sign-compare -Wno-unused -Wno-builtin-macro-redefined

src/Makevars.ucrt

@@ -1,4 +1,4 @@
-PKG_CXXFLAGS=  -I./code_base -I./include $(SHLIB_OPENMP_CXXFLAGS) -I$(R_TOOLS_SOFT)/include/nlopt -I$(R_TOOLS_SOFT)/include/gsl -DR_COMPILATION -ftree-vectorize -Os -O3 -Wno-ignored-attributes -Wno-sign-compare -Wno-unused -Wa,-mbig-obj -Wno-builtin-macro-redefined -Wno-comment -DEIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS
+PKG_CXXFLAGS=  -I./code_base -I./include $(SHLIB_OPENMP_CXXFLAGS) -I$(R_TOOLS_SOFT)/include/nlopt -I$(R_TOOLS_SOFT)/include/gsl -DR_COMPILATION -ftree-vectorize -Os -Wno-ignored-attributes -Wno-sign-compare -Wno-unused -Wa,-mbig-obj -Wno-builtin-macro-redefined -Wno-comment -DEIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS -DEIGEN_WARNINGS_DISABLED
 PKG_LIBS= $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) $(SHLIB_OPENMP_CXXFLAGS) -lgsl -lgslcblas -lnlopt
 CXX14FLAGS= -Wno-ignored-attributes -Wno-sign-compare -Wno-unused -Wno-builtin-macro-redefined
 CXX11FLAGS= -Wno-ignored-attributes -Wno-sign-compare -Wno-unused -Wno-builtin-macro-redefined

src/Makevars.win

@@ -1,4 +1,4 @@
-PKG_CPPFLAGS= -I../windows/gsl-2.7/include -I../windows/nlopt-2.7.1/include  -I./ -I./code_base  -I./include -DR_COMPILATION -ftree-vectorize -Os -O3 -Wno-ignored-attributes -Wno-sign-compare -Wno-unused -Wa,-mbig-obj -Wno-builtin-macro-redefined -Wno-comment -DEIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS
+PKG_CPPFLAGS= -I../windows/gsl-2.7/include -I../windows/nlopt-2.7.1/include  -I./ -I./code_base  -I./include -DR_COMPILATION -ftree-vectorize -Os -DEIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS -DEIGEN_WARNINGS_DISABLED
 PKG_LIBS= -L../windows/nlopt-2.7.1/lib${R_ARCH}${CRT}  -L../windows/gsl-2.7/lib${R_ARCH}${CRT} -lgsl -lgslcblas -lnlopt
 
 CXX14FLAGS= -Wno-ignored-attributes -Wno-sign-compare -Wno-unused -Wno-builtin-macro-redefined

src/RcppExports.cpp

@@ -149,6 +149,16 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// setseedGSL
+void setseedGSL(int s);
+RcppExport SEXP _ToxicR_setseedGSL(SEXP sSEXP) {
+BEGIN_RCPP
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< 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/code_base/IDPrior.cpp

@@ -83,7 +83,7 @@ double IDPrior::neg_log_prior(Eigen::MatrixXd theta, bool bound_check = false)
 		}
 	}
 	return -1.0 * returnV;
-};
+}
 
 ///////////////////////////////////////////////////////////
 // Function: neg_log_prior(Eigen::MatrixXd theta)
@@ -136,7 +136,7 @@ Eigen::MatrixXd IDPrior::log_prior(Eigen::MatrixXd theta)
 		}
 	}
 	return -1 * returnV.asDiagonal();
-};
+}
 
 //////////////////////////////////////////////////////
 // Function: prior_mean()
@@ -162,4 +162,4 @@ Eigen::MatrixXd IDPrior::prior_mean()
 		}
 	}
 	return pmean;
-};
+}

src/include/bmd_calculate.h

@@ -630,7 +630,7 @@ bmd_analysis bmd_fast_BMD_cont(LL likelihood, PR prior,
 	 */
 	std::vector<double> x(500);
 	std::vector<double> y(500);
-	if (isnormal(var(0, 0)) && (var(0.0) > 1e-7) && isnormal(log(BMD)))
+	if (isnormal(var(0, 0)) && (var(0, 0) > 1e-7) && isnormal(log(BMD)))
 	{
 
 		for (size_t i = 0; i < x.size(); i++)

src/include/bmds_entry.h

@@ -14,7 +14,9 @@
 typedef BSTR BMDS_MODEL_ID;
 #else
 #define BMDS_ENTRY_API
+#ifndef _stdcall
 #define _stdcall
+#endif
 typedef char *BMDS_MODEL_ID;
 #endif // defined WIN32 || defined _WINDOWS
 

src/include/mcmc_analysis.h

@@ -41,12 +41,11 @@
 #else
 #include <Eigen/Dense>
 #endif
-
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_spline.h>
-#include <gsl/gsl_randist.h>
 #include <gsl/gsl_cdf.h>
 
+#include "seeder.h"
 #include "dBMDstatmod.h"
 #include "cBMDstatmod.h"
 #include "bmd_calculate.h"
@@ -89,14 +88,12 @@ mcmcSamples MCMC_bmd_analysis_DNC(Eigen::MatrixXd Y, Eigen::MatrixXd D, Eigen::M
   //  ziggurat is used as it is the fastest sampler algorithm gsl has
   Eigen::MatrixXd rNormal(n, samples);
   rNormal.setZero(); 
-  gsl_rng *r;
-  // gsl_rng_env_setup();
-  r = gsl_rng_alloc(gsl_rng_mt19937);
+  Seeder* seeder = Seeder::getInstance();
   for (int i = 0; i < samples; i++)
   {
     for (int j = 0; j < n; j++)
     {
-      rNormal(j, i) = gsl_ran_gaussian_ziggurat(r, 1.0);
+      rNormal(j, i) = seeder->get_gaussian_ziggurat();
     }
   }
 
@@ -138,7 +135,7 @@ mcmcSamples MCMC_bmd_analysis_DNC(Eigen::MatrixXd Y, Eigen::MatrixXd D, Eigen::M
       double denom = penLike(0, i - 1) + b(0, 0);
 
       double test = exp(numer - denom);
-      double rr = gsl_rng_uniform(r);
+      double rr = seeder->get_uniform();
 
       if (isnan(test))
       {
@@ -167,7 +164,6 @@ mcmcSamples MCMC_bmd_analysis_DNC(Eigen::MatrixXd Y, Eigen::MatrixXd D, Eigen::M
                         : model.log_likelihood.compute_BMD_ADDED_NC(nSamples.col(i), BMR);
   }
 
-  gsl_rng_free(r);
   /////////////////////////////////////////////////////////////////
   rVal.BMD = BMD;          // vector of burn in BMD samples
   rVal.samples = nSamples; // vector of sample parameters
@@ -228,15 +224,12 @@ mcmcSamples mcmc_continuous(cBMDModel<LL, PR> *model, int samples,
   Eigen::MatrixXd rNormal(n, samples);
   rNormal.setZero(); 
 
-  gsl_rng *r;
-  // gsl_rng_env_setup();
-  r = gsl_rng_alloc(gsl_rng_mt19937);
-  gsl_rng_set(r, mySeed);
+  Seeder* seeder = Seeder::getInstance();
   for (int i = 0; i < samples; i++)
   {
     for (int j = 0; j < n; j++)
     {
-      rNormal(j, i) = gsl_ran_gaussian_ziggurat(r, 1.0);
+      rNormal(j, i) = seeder->get_gaussian_ziggurat();
     }
   }
 
@@ -284,7 +277,7 @@ mcmcSamples mcmc_continuous(cBMDModel<LL, PR> *model, int samples,
       double numer = -model->negPenLike(value) + a(0, 0);
       double denom = penLike(0, i - 1) + b(0, 0);
       double test = exp(numer - denom);
-      double rr = gsl_rng_uniform(r);
+      double rr = seeder->get_uniform();
 
       if (isnan(test))
       {
@@ -314,7 +307,6 @@ mcmcSamples mcmc_continuous(cBMDModel<LL, PR> *model, int samples,
                                  BMR, pointP);
   }
 
-  gsl_rng_free(r);
   /////////////////////////////////////////////////////////////////
   rVal.BMD = BMD;          // vector of BMD samples
   rVal.samples = nSamples; // vector of sample parameters

src/include/seeder.h

@@ -0,0 +1,72 @@
+#define STRICT_R_HEADERS
+#pragma once
+#ifndef SEEDER
+#define SEEDER
+
+#include <Rcpp.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+
+class Seeder {
+    private:
+    gsl_rng* r = nullptr;
+    static Seeder* instance;
+    int currentSeed = -1;
+    Seeder() {}
+
+    public:
+    Seeder(Seeder const&) = delete;
+    Seeder& operator=(Seeder const&) = delete;
+    static Seeder* getInstance() {
+        if(!instance) {
+            instance = new Seeder();
+        }
+        return instance;
+    }
+
+    ~Seeder() {
+        if(r) {
+            Rcpp::Rcout << "Singleton destroyed" << std::endl;
+            gsl_rng_free(r);
+        }
+    }
+
+    void setSeed(int seed) {
+        if (seed < 0) {
+            Rcpp::stop("Error: Seed must be a positive integer.");
+        }
+
+        // if(currentSeed > 0) {
+        //     Rcpp::Rcout << "Previous GSL seed: " << currentSeed << std::endl;
+        // }
+
+        if(r) {
+            // Rcpp::Rcout << "Freeing r before setting it to a new value" << std::endl;
+            gsl_rng_free(r);
+        }
+        gsl_rng_env_setup();
+        r = gsl_rng_alloc(gsl_rng_mt19937);
+        gsl_rng_set(r, seed);
+        if (currentSeed != seed) {  
+            Rcpp::Rcout << "Updated GSL seed: " << currentSeed << std::endl;
+        }
+        currentSeed = seed;
+    } 
+
+    double get_uniform() {
+        return gsl_rng_uniform(r);
+    }
+
+    double get_gaussian_ziggurat() {
+        return gsl_ran_gaussian_ziggurat(r, 1.0);
+    }
+
+    double get_ran_flat() {
+        return gsl_ran_flat(r, -1, 1);
+    }
+
+    gsl_rng* get_gsl_rng() {
+        return r;
+    }
+};
+#endif