File: TagAlignments.cpp

package info (click to toggle)
bitseq 0.7.5+dfsg-5
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 676 kB
  • sloc: cpp: 7,043; python: 562; makefile: 150; sh: 52
file content (132 lines) | stat: -rw-r--r-- 4,486 bytes parent folder | download | duplicates (4)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
#include<cmath>

#include "TagAlignments.h"

#include "misc.h"

#include "common.h"

//#define MEM_USAGE

TagAlignments::TagAlignments(bool storeL){//{{{
   knowNtotal=false;
   knowNreads=false;
   Ntotal=0;
   Nreads=0;
   storeLog = storeL;
}//}}}
void TagAlignments::init(long Nreads,long Ntotal, long M){//{{{
   currentRead = 0;
   reservedN = 0;
   if(Nreads>0){
      this->Nreads=Nreads;
      knowNreads=true;
      readIndex.reserve(Nreads+2);
   }      
   readIndex.push_back(0);
   
   if(Ntotal>0){
      this->Ntotal=Ntotal;
      knowNtotal=true;
      reservedN = Ntotal+1;
      trIds.reserve(reservedN);
      probs.reserve(reservedN);
   }
   if(M>0){
      this->M=M;
      readsInIsoform.assign(M,-1);
   }else{
      readsInIsoform.clear();
      this->M=0;
   }
}//}}}
void TagAlignments::pushAlignment(long trId, double prob){//{{{
   if(prob<=0)pushAlignmentL(trId, ns_misc::LOG_ZERO);
   else pushAlignmentL(trId, log(prob));
}//}}}
void TagAlignments::pushAlignmentL(long trId, double lProb){//{{{
   if(trId>=M){
      M=trId+1;
      readsInIsoform.resize(M,-1);
   }
   if(readsInIsoform[trId] == currentRead){
      // The read has already one alignment to this transcript.
     for(long i=readIndex[currentRead];i<(long)trIds.size();i++)
        if(trIds[i] == trId){
           probs[i] = ns_math::logAddExp(probs[i], lProb);
           break;
        }
   }else{
      if(! knowNtotal){
         // the size of arrays is unknown try to reserve sensible amount of space if we know Nreads
         if(knowNreads && reservedN && ((long)probs.size() == reservedN)){
            // we reached the size of reserved space
            double dens = (double)probs.size() / currentRead; 
            dens *= 1.05; //increase it by 5%
            reservedN =(long)( reservedN + (dens) * (Nreads - currentRead + 1000.0) );
         #ifdef MEM_USAGE
            message("TagAlignments:\n   size: %ld  reserving: %ld  capacity before: %ld\n",probs.size(),reservedN,probs.capacity());
         #endif
            trIds.reserve(reservedN);
            probs.reserve(reservedN);
         #ifdef MEM_USAGE
            message("   capacity after: %ld\n",probs.capacity());
         #endif
         }else if(knowNreads && (! reservedN) && (currentRead == Nreads / 4 )){
            // one quarter in, try to reserve sensible amount of space
            double dens = (double)probs.size() / currentRead; 
            dens *= 1.05; //increase it by 5%
            reservedN =(long)((dens) * (Nreads));
         #ifdef MEM_USAGE
            message("TagAlignments:\n   size: %ld  reserving: %ld  capacity before: %ld\n",probs.size(),reservedN,probs.capacity());
         #endif
            trIds.reserve(reservedN);
            probs.reserve(reservedN);
         #ifdef MEM_USAGE
            message("   capacity after: %ld\n",probs.capacity());
         #endif
         }
      }
      trIds.push_back(trId);
      probs.push_back(lProb);
      // Mark that transcript trId already has alignment from this read.
      readsInIsoform[trId] = currentRead;
   }
}//}}}
void TagAlignments::pushRead(){//{{{
   // Check whether there were any valid alignments added for this read:
   if(readIndex[currentRead] == (int_least32_t) probs.size()){
      // If no new alignments, do nothing.
      return;
   }
   // If there are alignments transform from log space if necessary and move to next read.
   if(!storeLog){
      double logSum = ns_math::logSumExp(probs, readIndex[currentRead], probs.size());
      for(long i = readIndex[currentRead]; i<(long)probs.size(); i++)
         probs[i] = exp(probs[i]-logSum);
   }
   // Move to the next read.
   currentRead++;
   readIndex.push_back(probs.size());
}//}}}
void TagAlignments::finalizeRead(long *M, long *Nreads, long *Ntotal){//{{{
   *M = this->M = readsInIsoform.size();
   *Nreads = this->Nreads = readIndex.size() - 1;
   *Ntotal = this->Ntotal = probs.size();
#ifdef MEM_USAGE
   message("TagAlignments: readIndex size: %ld  capacity %ld\n",readIndex.size(),readIndex.capacity());
   message("TagAlignments: probs size: %ld  capacity %ld\n",probs.size(),probs.capacity());
#endif
}//}}}
int_least32_t TagAlignments::getTrId(long i) const {//{{{
   if(i<Ntotal)return trIds[i];
   return 0;
}//}}}
double TagAlignments::getProb(long i) const {//{{{
   if(i<Ntotal)return probs[i];
   return 0;
}//}}}
int_least32_t TagAlignments::getReadsI(long i) const {//{{{
   if(i<=Nreads)return readIndex[i];
   return 0;
}//}}}