Fixed a bug in perl scripts for printing error messages

[rsem.git] / calcCI.cpp

diff --git a/calcCI.cpp b/calcCI.cpp
index 79a6bb2f924ddd0746e0cd816430541ad00e0370..99ce14858825c989216dbd6ab309562218a8ae62 100644 (file)
--- 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);