X-Git-Url: https://git.donarmstrong.com/?p=rsem.git;a=blobdiff_plain;f=calcCI.cpp;h=99ce14858825c989216dbd6ab309562218a8ae62;hp=79a6bb2f924ddd0746e0cd816430541ad00e0370;hb=683863b75f8d8bef2461039a6911b0e9619cc113;hpb=635ca2939cfb1f519f19e9dec072ddd05e9fb450

diff --git a/calcCI.cpp b/calcCI.cpp
index 79a6bb2..99ce148 100644
--- a/calcCI.cpp
+++ b/calcCI.cpp
@@ -193,13 +193,16 @@ void sample_theta_vectors_from_count_vectors() {
 	model.read(modelF);
 	calcExpectedEffectiveLengths<ModelType>(model);
 
+
+	int num_threads = min(nThreads, nCV);
+
 	buffer = new Buffer(nMB, nSamples, cvlen, tmpF);
 
-	paramsArray = new Params[nThreads];
-	threads = new pthread_t[nThreads];
+	paramsArray = new Params[num_threads];
+	threads = new pthread_t[num_threads];
 
 	char inpF[STRLEN];
-	for (int i = 0; i < nThreads; i++) {
+	for (int i = 0; i < num_threads; i++) {
 		paramsArray[i].no = i;
 		sprintf(inpF, "%s%d", cvsF, i);
 		paramsArray[i].fi = fopen(inpF, "r");
@@ -211,11 +214,11 @@ void sample_theta_vectors_from_count_vectors() {
 	pthread_attr_init(&attr);
 	pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
 
-	for (int i = 0; i < nThreads; i++) {
+	for (int i = 0; i < num_threads; i++) {
 		rc = pthread_create(&threads[i], &attr, &sample_theta_from_c, (void*)(&paramsArray[i]));
 		pthread_assert(rc, "pthread_create", "Cannot create thread " + itos(i) + " (numbered from 0) in sample_theta_vectors_from_count_vectors!");
 	}
-	for (int i = 0; i < nThreads; i++) {
+	for (int i = 0; i < num_threads; i++) {
 		rc = pthread_join(threads[i], &status);
 		pthread_assert(rc, "pthread_join", "Cannot join thread " + itos(i) + " (numbered from 0) in sample_theta_vectors_from_count_vectors!");
 	}
@@ -224,7 +227,7 @@ void sample_theta_vectors_from_count_vectors() {
 	pthread_attr_destroy(&attr);
 	delete[] threads;
 
-	for (int i = 0; i < nThreads; i++) {
+	for (int i = 0; i < num_threads; i++) {
 		fclose(paramsArray[i].fi);
 		delete paramsArray[i].engine;
 	}
@@ -318,24 +321,25 @@ void* calcCI_batch(void* arg) {
 void calculate_credibility_intervals(char* imdName) {
 	FILE *fo;
 	char outF[STRLEN];
+	int num_threads = nThreads;
 
 	iso_tau = new CIType[M + 1];
 	gene_tau = new CIType[m];
 
 	assert(M > 0);
-	int quotient = M / nThreads;
-	if (quotient < 1) { nThreads = M; quotient = 1; }
+	int quotient = M / num_threads;
+	if (quotient < 1) { num_threads = M; quotient = 1; }
 	int cur_gene_id = 0;
 	int num_isoforms = 0;
 
 	// A just so so strategy for paralleling
-	ciParamsArray = new CIParams[nThreads];
-	for (int i = 0; i < nThreads; i++) {
+	ciParamsArray = new CIParams[num_threads];
+	for (int i = 0; i < num_threads; i++) {
 		ciParamsArray[i].no = i;
 		ciParamsArray[i].start_gene_id = cur_gene_id;
 		num_isoforms = 0;
 
-		while ((m - cur_gene_id > nThreads - i - 1) && (i == nThreads - 1 || num_isoforms < quotient)) {
+		while ((m - cur_gene_id > num_threads - i - 1) && (i == num_threads - 1 || num_isoforms < quotient)) {
 			num_isoforms += gi.spAt(cur_gene_id + 1) - gi.spAt(cur_gene_id);
 			++cur_gene_id;
 		}
@@ -343,18 +347,18 @@ void calculate_credibility_intervals(char* imdName) {
 		ciParamsArray[i].end_gene_id = cur_gene_id;
 	}
 
-	threads = new pthread_t[nThreads];
+	threads = new pthread_t[num_threads];
 
 	/* set thread attribute to be joinable */
 	pthread_attr_init(&attr);
 	pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
 
 	// paralleling
-	for (int i = 0; i < nThreads; i++) {
+	for (int i = 0; i < num_threads; i++) {
 		rc = pthread_create(&threads[i], &attr, &calcCI_batch, (void*)(&ciParamsArray[i]));
 		pthread_assert(rc, "pthread_create", "Cannot create thread " + itos(i) + " (numbered from 0) in calculate_credibility_intervals!");
 	}
-	for (int i = 0; i < nThreads; i++) {
+	for (int i = 0; i < num_threads; i++) {
 		rc = pthread_join(threads[i], &status);
 		pthread_assert(rc, "pthread_join", "Cannot join thread " + itos(i) + " (numbered from 0) in calculate_credibility_intervals!");
 	}
@@ -393,10 +397,13 @@ void calculate_credibility_intervals(char* imdName) {
 
 int main(int argc, char* argv[]) {
 	if (argc < 8) {
-		printf("Usage: rsem-calculate-credibility-intervals reference_name sample_name sampleToken confidence nCV nSpC nMB [-p #Threads] [-q]\n");
+		printf("Usage: rsem-calculate-credibility-intervals reference_name imdName statName confidence nCV nSpC nMB [-p #Threads] [-q]\n");
 		exit(-1);
 	}
 
+	strcpy(imdName, argv[2]);
+	strcpy(statName, argv[3]);
+
 	confidence = atof(argv[4]);
 	nCV = atoi(argv[5]);
 	nSpC = atoi(argv[6]);
@@ -410,11 +417,6 @@ int main(int argc, char* argv[]) {
 	}
 	verbose = !quiet;
 
-	if (nThreads > nCV) {
-		nThreads = nCV;
-		printf("Warning: Number of count vectors is less than number of threads! Change the number of threads to %d!\n", nThreads);
-	}
-
 	sprintf(refF, "%s.seq", argv[1]);
 	refs.loadRefs(refF, 1);
 	M = refs.getM();
@@ -426,8 +428,6 @@ int main(int argc, char* argv[]) {
 	cvlen = M + 1;
 	assert(nSamples > 0 && cvlen > 1); // for Buffter.h: (bufsize_type)nSamples
 
-	sprintf(imdName, "%s.temp/%s", argv[2], argv[3]);
-	sprintf(statName, "%s.stat/%s", argv[2], argv[3]);
 	sprintf(tmpF, "%s.tmp", imdName);
 	sprintf(cvsF, "%s.countvectors", imdName);