Mapping speed improvement

source/FastResetVector.h

@@ -0,0 +1,40 @@
+#ifndef H_FastResetVector
+#define H_FastResetVector
+
+#include <vector>
+
+// Implementation of a vector that can be reset to a default value in ~O(1).
+// It is more efficient than an ordinary vector if:
+// - the values are sparse, such that only a few elements need to be deleted.
+// - there are not too many accesses to elements, since upon each access
+//   it must be checked if the element is stale.
+
+template <class T> class FastResetVector {
+
+    private:
+        vector<T> data; // contains actual data to be stored
+        T defaultValue; // all elements of `data` are initialized with this value
+        unsigned int incarnation = 0; // increasing the incarnation invalidates all elements of `data` (=reset)
+        vector<unsigned int> lastUpdate; // for each element in `data`, keep track of the incarnation that it was last updated
+
+    public:
+        FastResetVector(const size_t s, const T& d): data(s,d), defaultValue(d), lastUpdate(s) {}
+
+        inline T& operator[](const size_t i) { // whenever an element is accessed, check if it's stale
+            if (incarnation != lastUpdate[i]) { // is it stale?
+                data[i] = defaultValue; // reset to defaut value
+                lastUpdate[i] = incarnation; // mark as fresh for this incarnation
+            };
+            return data[i];
+        }
+
+        void reset() {
+            incarnation++; // we can invalidate `data` simply through "reincarnation"
+            if (incarnation == 0) // only when there is an integer overflow, we need to reinitialize
+                fill(data.begin(), data.end(), defaultValue);
+        }
+
+};
+
+#endif
+

source/PackedArray.h

@@ -27,7 +27,8 @@ inline uint PackedArray::operator [] (uint ii) {
    uint S=b%8;
 
    uint a1 = *((uint*) (charArray+B));
-   a1 = ((a1>>S)<<wordCompLength)>>wordCompLength;
+   a1 >>= S;
+   a1 &= bitRecMask;
    return a1;
 };
 

source/ReadAlign.cpp

@@ -7,11 +7,7 @@ ReadAlign::ReadAlign (Parameters& Pin, Genome &genomeIn, Transcriptome *TrIn, in
                     : mapGen(genomeIn), P(Pin), chunkTr(TrIn)
 {
     readNmates=P.readNmates;
-    winBin = new uintWinBin* [2];
-    winBin[0] = new uintWinBin [P.winBinN];
-    winBin[1] = new uintWinBin [P.winBinN];
-    memset(winBin[0],255,sizeof(winBin[0][0])*P.winBinN);
-    memset(winBin[1],255,sizeof(winBin[0][0])*P.winBinN);
+    winBin.resize(2, FastResetVector<uintWinBin>(P.winBinN, uintWinBinMax));
     //RNGs
     rngMultOrder.seed(P.runRNGseed*(iChunk+1));
     rngUniformReal0to1=std::uniform_real_distribution<double> (0.0, 1.0);

source/ReadAlign.h

@@ -13,6 +13,7 @@
 #include "ChimericDetection.h"
 #include "SoloRead.h"
 #include "ReadAnnotations.h"
+#include "FastResetVector.h"
 
 #include <time.h>
 #include <random>
@@ -111,7 +112,7 @@ class ReadAlign {
 //         uint fragLength[MAX_N_FRAG], fragStart[MAX_N_FRAG]; //fragment Lengths and Starts in read space
 
         //binned alignments
-        uintWinBin **winBin; //binned genome: window ID (number) per bin
+        vector< FastResetVector<uintWinBin> > winBin; //binned genome: window ID (number) per bin
 
         //alignments
         uiPC *PC; //pieces coordinates

source/ReadAlign_createExtendWindowsWithAlign.cpp

@@ -9,7 +9,7 @@ int ReadAlign::createExtendWindowsWithAlign(uint a1, uint aStr) {
 
     uint aBin = (a1 >> P.winBinNbits); //align's bin
     uint iBinLeft=aBin, iBinRight=aBin;
-    uintWinBin* wB=winBin[aStr];
+    FastResetVector<uintWinBin>& wB=winBin[aStr];
     uint iBin=-1, iWin=-1, iWinRight=-1;
 
 

source/ReadAlign_stitchPieces.cpp

@@ -13,9 +13,8 @@
 void ReadAlign::stitchPieces(char **R, uint Lread) {
 
     //zero-out winBin
-    memset(winBin[0],255,sizeof(winBin[0][0])*P.winBinN);
-    memset(winBin[1],255,sizeof(winBin[0][0])*P.winBinN);
-
+    winBin[0].reset();
+    winBin[1].reset();
 
 //     for (uint iWin=0;iWin<nWall;iWin++) {//zero out winBin
 //         if (WC[iWin][WC_gStart]<=WC[iWin][WC_gEnd]) {//otherwise the window is dead