rsem v1.1.14

[rsem.git] / EM.cpp

diff --git a/EM.cpp b/EM.cpp
index 2921ad33cd837a84241a5c6f8c410e1d5bacc5d7..cf53fc3881b557631dca0ff6fdc4d0bfdd209b5d 100644 (file)
--- a/EM.cpp
+++ b/EM.cpp
@@ -102,7 +102,7 @@ void init(ReadReader<ReadType> **&readers, HitContainer<HitType> **&hitvs, doubl
                indices[i] = new ReadIndex(readFs[i]);
        }
        for (int i = 0; i < nThreads; i++) {
-               readers[i] = new ReadReader<ReadType>(s, readFs);
+               readers[i] = new ReadReader<ReadType>(s, readFs, refs.hasPolyA(), mparams.seedLen); // allow calculation of calc_lq() function
                readers[i]->setIndices(indices);
        }
 
@@ -340,7 +340,7 @@ void writeResults(ModelType& model, double* counts) {
                fprintf(fo, "%.15g%c", tau[i], (i < M ? '\t' : '\n'));
        for (int i = 1; i <= M; i++) {
                const Transcript& transcript = transcripts.getTranscriptAt(i);
-               fprintf(fo, "%s%c", transcript.getLeft().c_str(), (i < M ? '\t' : '\n'));
+               fprintf(fo, "%s%c", transcript.getGeneID().c_str(), (i < M ? '\t' : '\n'));
        }
        fclose(fo);
 
@@ -465,14 +465,18 @@ void EM() {
                //E step
                for (int i = 0; i < nThreads; i++) {
                        rc = pthread_create(&threads[i], &attr, E_STEP<ReadType, HitType, ModelType>, (void*)(&fparams[i]));
-                       if (rc != 0) { fprintf(stderr, "Cannot create thread %d at ROUND %d! (numbered from 0)\n", i, ROUND); exit(-1); }
-                       //assert(rc == 0);
+                       if (rc != 0) {
+                               fprintf(stderr, "Cannot create thread %d at ROUND %d! (numbered from 0)", i, ROUND);
+                               pthread_exception(rc);
+                       }
                }
 
                for (int i = 0; i < nThreads; i++) {
                        rc = pthread_join(threads[i], &status);
-                       if (rc != 0) { fprintf(stderr, "Cannot join thread %d at ROUND %d! (numbered from 0)\n", i, ROUND); exit(-1); }
-                       //assert(rc == 0);
+                       if (rc != 0) {
+                               fprintf(stderr, "Cannot join thread %d at ROUND %d! (numbered from 0)\n", i, ROUND);
+                               pthread_exception(rc);
+                       }
                }
 
                model.setNeedCalcConPrb(false);
@@ -518,11 +522,17 @@ void EM() {
                if (model.getNeedCalcConPrb()) {
                        for (int i = 0; i < nThreads; i++) {
                                rc = pthread_create(&threads[i], &attr, calcConProbs<ReadType, HitType, ModelType>, (void*)(&fparams[i]));
-                               if (rc != 0) { fprintf(stderr, "Cannot create thread %d when generate files for Gibbs sampler! (numbered from 0)\n", i); exit(-1); }
+                               if (rc != 0) {
+                                       fprintf(stderr, "Cannot create thread %d when generate files for Gibbs sampler! (numbered from 0)\n", i);
+                                       pthread_exception(rc);
+                               }
                        }
                        for (int i = 0; i < nThreads; i++) {
                                rc = pthread_join(threads[i], &status);
-                               if (rc != 0) { fprintf(stderr, "Cannot join thread %d when generate files for Gibbs sampler! (numbered from 0)\n", i); exit(-1); }
+                               if (rc != 0) {
+                                       fprintf(stderr, "Cannot join thread %d when generate files for Gibbs sampler! (numbered from 0)\n", i);
+                                       pthread_exception(rc);
+                               }
                        }
                }
                model.setNeedCalcConPrb(false);
@@ -576,13 +586,17 @@ void EM() {
        for (int i = 0; i <= M; i++) probv[i] = theta[i];
        for (int i = 0; i < nThreads; i++) {
                rc = pthread_create(&threads[i], &attr, E_STEP<ReadType, HitType, ModelType>, (void*)(&fparams[i]));
-               if (rc != 0) { fprintf(stderr, "Cannot create thread %d when calculate expected weights! (numbered from 0)\n", i); exit(-1); }
-               //assert(rc == 0);
+               if (rc != 0) {
+                       fprintf(stderr, "Cannot create thread %d when calculate expected weights! (numbered from 0)\n", i);
+                       pthread_exception(rc);
+               }
        }
        for (int i = 0; i < nThreads; i++) {
                rc = pthread_join(threads[i], &status);
-               if (rc != 0) { fprintf(stderr, "Cannot join thread %d! (numbered from 0) when calculate expected weights!\n", i); exit(-1); }
-               //assert(rc == 0);
+               if (rc != 0) {
+                       fprintf(stderr, "Cannot join thread %d! (numbered from 0) when calculate expected weights!\n", i);
+                       pthread_exception(rc);
+               }
        }
        model.setNeedCalcConPrb(false);
        for (int i = 1; i < nThreads; i++) {
@@ -596,7 +610,7 @@ void EM() {
        pthread_attr_destroy(&attr);
 
        //convert theta' to theta
-               double *mw = model.getMW();
+    double *mw = model.getMW();
        sum = 0.0;
        for (int i = 0; i <= M; i++) {
          theta[i] = (mw[i] < EPSILON ? 0.0 : theta[i] / mw[i]);