Skip to content
Permalink
Browse files

Fix clang 3.6 compile error

  • Loading branch information...
milot-mirdita committed Aug 2, 2019
1 parent 4778002 commit a98349156f9664cac4e3fc7e6df213169a2b82bc
Showing with 14 additions and 57 deletions.
  1. +14 −57 src/assembler/correctreads.cpp
@@ -1,22 +1,18 @@
#include <stdlib.h>
#include <mmseqs/lib/omptl/omptl_algorithm>
#include "Timer.h"
#include "NucleotideMatrix.h"
#include "ReducedMatrix.h"
#include "Indexer.h"
#include "Sequence.h"
#include "SubstitutionMatrix.h"
#include "simd.h"
#include "Debug.h"
#include "AminoAcidLookupTables.h"
#include "DBReader.h"
#include "DBWriter.h"
#include "LocalParameters.h"
#include "simd.h"

#ifdef OPENMP
#include <omp.h>
#endif

#include <omptl/omptl_algorithm>
#include <cstdlib>

#define HI_NIBBLE(b) (((b) >> 4) & 0x0F)
#define LO_NIBBLE(b) ((b) & 0x0F)

@@ -56,46 +52,6 @@ std::pair<size_t, char> getSubstituion(const char lastLetter, const size_t currK

void printKmer(size_t idx, int size);

// Compute reverse complement of k-mer in 2-bit-per-nucleotide encoding (A: 00, C: 01, G: 10, T: 11)
uint64_t revComplement(const uint64_t kmer, const int k) {
// broadcast 64bit to 128 bit
__m128i x = _mm_cvtsi64_si128(kmer);

// create lookup (set 16 bytes in 128 bit)
// a lookup entry at the index of two nucleotids (4 bit) describes the reverse
// complement of these two nucleotid in the higher 4 bits (lookup1) or in the
// lower 4 bits (lookup2)
__m128i lookup1 = _mm_set_epi8(0x50,0x10,0xD0,0x90,0x40,0x00,0xC0,0x80,0x70,
0x30,0xF0,0xB0,0x60,0x20,0xE0,0xA0);
__m128i lookup2 = _mm_set_epi8(0x05,0x01,0x0D,0x09,0x04,0x00,0x0C,0x08,0x07,
0x03,0x0F,0x0B,0x06,0x02,0x0E,0x0A);

// _mm_set1_epi8: create 128 bit with all bytes set to given value
// here: 0x0F (00001111) and 0xF0 (11110000)
// _mm_and_si128: bitwise AND
__m128i kmer1 = _mm_and_si128(x, _mm_set1_epi8(0x0F)); // get lower 4 bits
__m128i kmer2 = _mm_and_si128(x, _mm_set1_epi8(0xF0)); // get higher 4 bits

// shift right by 2 nucleotids
kmer2 >>= 4;

// use _mm_shuffle_epi8 to look up reverse complement
kmer1 = _mm_shuffle_epi8(lookup1, kmer1);
kmer2 = _mm_shuffle_epi8(lookup2, kmer2);

// _mm_or_si128: bitwise OR
x = _mm_or_si128(kmer1, kmer2);

// set upper 8 bytes to 0 and revert order of lower 8 bytes
x = _mm_shuffle_epi8(x,
_mm_set_epi8(0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0,1,2,3,4,5,6,7));

// shift out the unused nucleotide positions (1 <= k <=32 )
// broadcast 128 bit to 64 bit
return (((uint64_t)_mm_cvtsi128_si64(x)) >> (uint64_t)(64-2*k));

}

int correctreads(int argc, const char **argv, const Command& command) {
LocalParameters& par = LocalParameters::getLocalInstance();
par.kmerSize = 5;
@@ -118,19 +74,23 @@ int correctreads(int argc, const char **argv, const Command& command) {
Kmer * allKmers = new Kmer[seqDb.getAminoAcidDBSize()*2];
// Create a 1D Tensor on length 20 for input data.
Debug(Debug::INFO) << "Extract kmers\n";

Debug::Progress progress(seqDb.getSize());
#pragma omp parallel
{
Kmer * threadKmerBuffer = new Kmer[BUFFER_SIZE];
size_t bufferPos = 0;
unsigned int thread_idx = 0;
#ifdef OPENMP
thread_idx = static_cast<unsigned int>(omp_get_thread_num());
#endif

Kmer* threadKmerBuffer = new Kmer[BUFFER_SIZE];
size_t bufferPos = 0;
Sequence seq(par.maxSeqLen, seqDb.getDbtype(), &subMat, par.kmerSize, false, false);

#pragma omp for schedule(static)
for (size_t id = 0; id < seqDb.getSize(); id++) {
progress.updateProgress();

char *seqData = seqDb.getData(id, thread_idx);
unsigned int seqLen = seqDb.getSeqLens(id)-2;

@@ -152,7 +112,7 @@ int correctreads(int argc, const char **argv, const Command& command) {
threadKmerBuffer[bufferPos].isReverse=false;
bufferPos++;

size_t revKmerIdx = revComplement(kmerIdx, par.kmerSize);
size_t revKmerIdx = Util::revComplement(kmerIdx, par.kmerSize);
// printKmer(kmerIdx, par.kmerSize);
// printKmer(revKmerIdx, par.kmerSize);
threadKmerBuffer[bufferPos].kmer = revKmerIdx;
@@ -174,16 +134,13 @@ int correctreads(int argc, const char **argv, const Command& command) {
memcpy(allKmers + writeOffset, threadKmerBuffer, sizeof(Kmer) * bufferPos);
}

delete [] threadKmerBuffer;
delete[] threadKmerBuffer;
}
Debug(Debug::INFO) << "Done\n";
Debug(Debug::INFO) << "Time for extracting kmers: " << timer.lap() << "\n";
timer.reset();

Debug(Debug::INFO) << "Sort kmer ... ";
timer.reset();
omptl::sort(allKmers, allKmers + offset, Kmer::compareRepSequenceAndIdAndPos);

Debug(Debug::INFO) << "Done\n";
Debug(Debug::INFO) << "Time for sort: " << timer.lap() << "\n";

if(allKmers[0].kmer != allKmers[1].kmer){
@@ -300,4 +257,4 @@ std::pair<size_t, char> getSubstituion(const char lastLetter, const size_t currK
}

#undef HI_NIBBLE
#undef LO_NIBBLE
#undef LO_NIBBLE

0 comments on commit a983491

Please sign in to comment.
You can’t perform that action at this time.