X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=calcCI.cpp;h=7cb58ffa07c24f5c6eb77f871e847aeb430ae7bf;hb=b5ca7d97e07526eb3c66205e25fc643510a31fcf;hp=97eba64e1ea9965c3983409090dd252abb78cbad;hpb=4f7502168c3816ba3283385f093e599527e2b144;p=rsem.git

diff --git a/calcCI.cpp b/calcCI.cpp
index 97eba64..7cb58ff 100644
--- a/calcCI.cpp
+++ b/calcCI.cpp
@@ -80,6 +80,9 @@ pthread_t *threads;
 pthread_attr_t attr;
 int rc;
 
+bool hasSeed;
+seedType seed;
+
 CIParams *ciParamsArray;
 
 void* sample_theta_from_c(void* arg) {
@@ -103,18 +106,22 @@ void* sample_theta_from_c(void* arg) {
 	int cnt = 0;
 	while (fscanf(fi, "%d", &cvec[0]) == 1) {
 		for (int j = 1; j <= M; j++) assert(fscanf(fi, "%d", &cvec[j]) == 1);
+		assert(cvec[0] > 0);
 
 		++cnt;
 
 		for (int j = 0; j <= M; j++) {
-			gammas[j] = new gamma_dist(cvec[j]);
-			rgs[j] = new gamma_generator(*(params->engine), *gammas[j]);
+		  gammas[j] = NULL; rgs[j] = NULL;
+		  if (cvec[j] > 0) {
+		    gammas[j] = new gamma_dist(cvec[j]);
+		    rgs[j] = new gamma_generator(*(params->engine), *gammas[j]);
+		  }
 		}
-
+		  
 		for (int i = 0; i < nSpC; i++) {
 			double sum = 0.0;
 			for (int j = 0; j <= M; j++) {
-				theta[j] = ((j == 0 || (eel[j] >= EPSILON && mw[j] >= EPSILON)) ? (*rgs[j])() / mw[j] : 0.0);
+				theta[j] = ((j == 0 || (cvec[j] > 0 && eel[j] >= EPSILON && mw[j] >= EPSILON)) ? (*rgs[j])() / mw[j] : 0.0);
 				sum += theta[j];
 			}
 			assert(sum >= EPSILON);
@@ -135,8 +142,8 @@ void* sample_theta_from_c(void* arg) {
 		}
 
 		for (int j = 0; j <= M; j++) {
-			delete gammas[j];
-			delete rgs[j];
+		  if (gammas[j] != NULL) delete gammas[j];
+		  if (rgs[j] != NULL) delete rgs[j];
 		}
 
 		if (verbose && cnt % 100 == 0) { printf("Thread %d, %d count vectors are processed!\n", params->no, cnt); }
@@ -165,6 +172,7 @@ void sample_theta_vectors_from_count_vectors() {
 	threads = new pthread_t[num_threads];
 
 	char inpF[STRLEN];
+	hasSeed ? engineFactory::init(seed) : engineFactory::init();
 	for (int i = 0; i < num_threads; i++) {
 		paramsArray[i].no = i;
 		sprintf(inpF, "%s%d", cvsF, i);
@@ -172,6 +180,7 @@ void sample_theta_vectors_from_count_vectors() {
 		paramsArray[i].engine = engineFactory::new_engine();
 		paramsArray[i].mw = model.getMW();
 	}
+	engineFactory::finish();
 
 	/* set thread attribute to be joinable */
 	pthread_attr_init(&attr);
@@ -458,7 +467,7 @@ void calculate_credibility_intervals(char* imdName) {
 
 int main(int argc, char* argv[]) {
 	if (argc < 8) {
-		printf("Usage: rsem-calculate-credibility-intervals reference_name imdName statName confidence nCV nSpC nMB [-p #Threads] [-q]\n");
+		printf("Usage: rsem-calculate-credibility-intervals reference_name imdName statName confidence nCV nSpC nMB [-p #Threads] [--seed seed] [-q]\n");
 		exit(-1);
 	}
 
@@ -473,8 +482,15 @@ int main(int argc, char* argv[]) {
 
 	nThreads = 1;
 	quiet = false;
+	hasSeed = false;
 	for (int i = 8; i < argc; i++) {
 		if (!strcmp(argv[i], "-p")) nThreads = atoi(argv[i + 1]);
+		if (!strcmp(argv[i], "--seed")) {
+		  hasSeed = true;
+		  int len = strlen(argv[i + 1]);
+		  seed = 0;
+		  for (int k = 0; k < len; k++) seed = seed * 10 + (argv[i + 1][k] - '0');
+		}
 		if (!strcmp(argv[i], "-q")) quiet = true;
 	}
 	verbose = !quiet;