From 59001c205c70d241ab33cab2893fb78f53ad637c Mon Sep 17 00:00:00 2001
From: John Marshall <John.W.Marshall@glasgow.ac.uk>
Date: Thu, 21 Mar 2019 14:40:32 +0000
Subject: [PATCH] Maintain ordering of reads in ReadBuffer::flush()

Having `m_filteredReads.push_back(m_rawreads[i])` in a loop that is
parallel for `i` meant that entries could be pushed onto `m_filteredReads`
in arbitrary order when multiple threads are used.

Instead store pointers in a fixed-size array in the parallel loop,
and squash them sequentially into `m_filteredReads` afterwards.
Fixes a bug whereby incorrect output was produced when multiple
threads were used (differences observed beyond six threads).
---
 src/read_buffer.cpp | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/src/read_buffer.cpp b/src/read_buffer.cpp
index bb2c6be5..5878ad67 100644
--- a/src/read_buffer.cpp
+++ b/src/read_buffer.cpp
@@ -36,6 +36,8 @@ ReadBuffer::~ReadBuffer()
 void ReadBuffer::flush()
 {
    // std::cout << "in flush " << std::endl;
+   std::vector<SPLIT_READ*> filteredReads(m_currentsize);
+
    #pragma omp parallel for
    for (int i=0; i<m_currentsize ; i++ ) {
       // std::cout << "before GetCloseEnd " << std::endl;
@@ -55,11 +57,11 @@ void ReadBuffer::flush()
          if (m_rawreads[i].hasCloseEnd()) {
             updateReadAfterCloseEndMapping(m_rawreads[i]);
 
-            #pragma omp critical
+            //#pragma omp critical
             {
                //			g_ReadSeq2Index.insert(std::pair<std::string, unsigned> (m_rawreads[i].UnmatchedSeq, m_filteredReads.size()));
                m_rawreads[i].SampleName2Number.insert(std::pair <std::string, unsigned> (m_rawreads[i].Tag, 1));
-               m_filteredReads.push_back(m_rawreads[i]);
+               filteredReads[i] = &m_rawreads[i];
                //std::cout << "with closed end: " << m_rawreads[i].Name << std::endl;
             }
 
@@ -94,6 +96,11 @@ void ReadBuffer::flush()
 
    }
    //std::cout << "end of flush " << std::endl;
+
+   for (size_t i = 0; i < filteredReads.size(); i++)
+      if (filteredReads[i])
+         m_filteredReads.push_back(*filteredReads[i]);
+
    m_rawreads.clear();
    m_currentsize = 0;
    std::cout << "The number of one end mapped read: " << m_filteredReads.size() << std::endl;