aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/bslmm.cpp33
-rw-r--r--src/bslmm.h7
-rw-r--r--src/bslmmdap.h4
-rw-r--r--src/param.cpp42
-rw-r--r--src/param.h7
5 files changed, 40 insertions, 53 deletions
diff --git a/src/bslmm.cpp b/src/bslmm.cpp
index 3305639..eb961e9 100644
--- a/src/bslmm.cpp
+++ b/src/bslmm.cpp
@@ -84,7 +84,8 @@ void BSLMM::CopyFromParam(PARAM &cPar) {
w_pace = cPar.w_pace;
n_mh = cPar.n_mh;
geo_mean = cPar.geo_mean;
- randseed = cPar.randseed;
+ // randseed = cPar.randseed;
+ gsl_r = cPar.gsl_r;
trace_G = cPar.trace_G;
ni_total = cPar.ni_total;
@@ -107,7 +108,7 @@ void BSLMM::CopyToParam(PARAM &cPar) {
cPar.cHyp_initial = cHyp_initial;
cPar.n_accept = n_accept;
cPar.pheno_mean = pheno_mean;
- cPar.randseed = randseed;
+ // cPar.randseed = randseed;
return;
}
@@ -938,19 +939,6 @@ void BSLMM::MCMC(const gsl_matrix *U, const gsl_matrix *UtX,
// Calculate proposal distribution for gamma (unnormalized),
// and set up gsl_r and gsl_t.
- gsl_rng_env_setup();
- const gsl_rng_type *gslType;
- gslType = gsl_rng_default;
- if (randseed < 0) {
- time_t rawtime;
- time(&rawtime);
- tm *ptm = gmtime(&rawtime);
-
- randseed =
- (unsigned)(ptm->tm_hour % 24 * 3600 + ptm->tm_min * 60 + ptm->tm_sec);
- }
- gsl_r = gsl_rng_alloc(gslType);
- gsl_rng_set(gsl_r, randseed);
double *p_gamma = new double[ns_test];
CalcPgamma(p_gamma);
@@ -1643,21 +1631,6 @@ void BSLMM::MCMC(const gsl_matrix *X, const gsl_vector *y) {
}
// Calculate proposal distribution for gamma (unnormalized),
- // and set up gsl_r and gsl_t.
- gsl_rng_env_setup();
- const gsl_rng_type *gslType;
- gslType = gsl_rng_default;
- if (randseed < 0) {
- time_t rawtime;
- time(&rawtime);
- tm *ptm = gmtime(&rawtime);
-
- randseed =
- (unsigned)(ptm->tm_hour % 24 * 3600 + ptm->tm_min * 60 + ptm->tm_sec);
- }
- gsl_r = gsl_rng_alloc(gslType);
- gsl_rng_set(gsl_r, randseed);
-
double *p_gamma = new double[ns_test];
CalcPgamma(p_gamma);
diff --git a/src/bslmm.h b/src/bslmm.h
index d2dadbf..93dadf9 100644
--- a/src/bslmm.h
+++ b/src/bslmm.h
@@ -20,7 +20,7 @@
#define __BSLMM_H__
#include <gsl/gsl_randist.h>
-#include <gsl/gsl_rng.h>
+// #include <gsl/gsl_rng.h>
#include <map>
#include <vector>
@@ -60,7 +60,8 @@ public:
size_t n_accept; // Number of acceptances.
size_t n_mh; // Number of MH steps per iter.
double geo_mean; // Mean of geometric dist.
- long int randseed;
+ // long int randseed;
+ gsl_rng *gsl_r; // Track randomizer state
double trace_G;
HYPBSLMM cHyp_initial;
@@ -88,7 +89,7 @@ public:
vector<SNPINFO> snpInfo;
// Not included in PARAM.
- gsl_rng *gsl_r;
+ // gsl_rng *gsl_r;
gsl_ran_discrete_t *gsl_t;
map<size_t, size_t> mapRank2pos;
diff --git a/src/bslmmdap.h b/src/bslmmdap.h
index dc05e34..0f560c2 100644
--- a/src/bslmmdap.h
+++ b/src/bslmmdap.h
@@ -21,7 +21,7 @@
#include "param.h"
#include <gsl/gsl_randist.h>
-#include <gsl/gsl_rng.h>
+// #include <gsl/gsl_rng.h>
#include <map>
#include <vector>
@@ -44,7 +44,7 @@ public:
double pheno_mean;
// BSLMM MCMC related parameters
- long int randseed;
+ // long int randseed;
double trace_G;
HYPBSLMM cHyp_initial;
diff --git a/src/param.cpp b/src/param.cpp
index e4781c4..db6d8d5 100644
--- a/src/param.cpp
+++ b/src/param.cpp
@@ -106,6 +106,10 @@ PARAM::PARAM(void)
time_total(0.0), time_G(0.0), time_eigen(0.0), time_UtX(0.0),
time_UtZ(0.0), time_opt(0.0), time_Omega(0.0) {}
+PARAM::~PARAM() {
+ gsl_rng_free(gsl_r);
+}
+
// Read files: obtain ns_total, ng_total, ns_test, ni_test.
void PARAM::ReadFiles(void) {
string file_str;
@@ -817,6 +821,28 @@ void PARAM::CheckParam(void) {
flag++;
}
+ // Always set up random environment.
+ gsl_rng_env_setup(); // sets gsl_rng_default_seed
+ const gsl_rng_type *T = gsl_rng_default; // pick up environment GSL_RNG_SEED
+
+ if (randseed >= 0)
+ gsl_rng_default_seed = randseed; // CLI option used
+ else if (gsl_rng_default_seed == 0) { // by default we will randomize the seed
+ time_t rawtime;
+ time(&rawtime);
+ tm *ptm = gmtime(&rawtime);
+
+ gsl_rng_default_seed =
+ (unsigned)(ptm->tm_hour % 24 * 3600 + ptm->tm_min * 60 + ptm->tm_sec);
+ }
+ gsl_r = gsl_rng_alloc(T);
+
+ if (is_debug_mode()) {
+ printf ("GSL random generator type: %s; ", gsl_rng_name (gsl_r));
+ printf ("seed = %lu (option %li); ", gsl_rng_default_seed, randseed);
+ printf ("first value = %lu\n", gsl_rng_get (gsl_r));
+ }
+
if (flag != 1 && a_mode != 15 && a_mode != 27 && a_mode != 28 &&
a_mode != 43 && a_mode != 5 && a_mode != 61 && a_mode != 62 &&
a_mode != 63 && a_mode != 66 && a_mode != 67) {
@@ -2015,22 +2041,6 @@ void PARAM::ProcessCvtPhen() {
<< "analyzed individuals. " << endl;
} else {
- // Set up random environment.
- gsl_rng_env_setup();
- gsl_rng *gsl_r;
- const gsl_rng_type *gslType;
- gslType = gsl_rng_default;
- if (randseed < 0) {
- time_t rawtime;
- time(&rawtime);
- tm *ptm = gmtime(&rawtime);
-
- randseed = (unsigned)(ptm->tm_hour % 24 * 3600 + ptm->tm_min * 60 +
- ptm->tm_sec);
- }
- gsl_r = gsl_rng_alloc(gslType);
- gsl_rng_set(gsl_r, randseed);
-
// From ni_test, sub-sample ni_subsample.
vector<size_t> a, b;
for (size_t i = 0; i < ni_subsample; i++) {
diff --git a/src/param.h b/src/param.h
index 9ad14b2..eb2cef7 100644
--- a/src/param.h
+++ b/src/param.h
@@ -23,6 +23,7 @@
#include "debug.h"
#include "gsl/gsl_matrix.h"
+#include <gsl/gsl_rng.h>
#include "gsl/gsl_vector.h"
#include <map>
#include <set>
@@ -211,7 +212,8 @@ public:
size_t n_accept; // Number of acceptance.
size_t n_mh; // # MH steps in each iter.
double geo_mean; // Mean of geometric dist.
- long int randseed;
+ long int randseed; // holds -seed parameter
+ gsl_rng *gsl_r; // Track the randomizer
double trace_G;
HYPBSLMM cHyp_initial;
@@ -324,8 +326,9 @@ public:
set<string> setKSnps; // Set of snps for K (-ksnps and LOCO)
set<string> setGWASnps; // Set of snps for GWA (-gwasnps and LOCO)
- // Constructor.
+ // Constructor and destructor
PARAM();
+ ~PARAM();
// Functions.
void ReadFiles();
generated by cgit v1.2.3 (git 2.44.0) at 2025-01-19 03:23:30 +0000