5
\$\begingroup\$

This code is meant to compute the reverse complement for a given sequence of bases. The code is tested versus some test cases and seems to work so far.

The primary goal is that it should be as fast as possible (e.g. recovery from "file not found" or format errors in the input file is not required).

The code is pretty fast already. But I would like to ask if someone sees any obvious opportunities for a worthwhile optimisation or obvious mistakes of course.

#include <iostream>
#include <chrono>

#include <cstdio>
#include <sys/stat.h>

#include <thread>
#include <queue>
#include <tuple>
#include <algorithm>

#define INPUTFILENAME "~/70000000.txt"
#define OUTPUTFILENAME "~/output.txt"

namespace helpers {
    size_t getFileSize(void)
    {
        struct stat stat_buf;
        int rc = stat(INPUTFILENAME, &stat_buf);
        return rc == 0 ? stat_buf.st_size : -1;
    }
}
char* static_front = 0;
char* static_back = 0;

std::mutex mutex;
std::condition_variable cv;
std::queue<std::tuple<char*, char*>> queue;

static inline void staticTransformation(char* front, char *back) {

    constexpr char complement[] = {
        ['A']='T', ['a']='T',
        ['C']='G', ['c']='G',
        ['G']='C', ['g']='C',
        ['T']='A', ['t']='A',
        ['U']='A', ['u']='A',
        ['M']='K', ['m']='K',
        ['R']='Y', ['r']='Y',
        ['W']='W', ['w']='W',
        ['S']='S', ['s']='S',
        ['Y']='R', ['y']='R',
        ['K']='M', ['k']='M',
        ['V']='B', ['v']='B',
        ['H']='D', ['h']='D',
        ['D']='H', ['d']='H',
        ['B']='V', ['b']='V',
        ['N']='N', ['n']='N'};

    char help;
    while(front < back) {
        if(*front=='\n') {
            front++;
        }
        if(*back=='\n') {
           back--;
        }
        // swap after complement
        help = complement[(*front)];
        *front = complement[(*back)];
        *back = help;
        front++;
        back--;
    }
}

void staticTransformationWorker(int name) {
    char* _front;
    char* _back;
    auto tuple = std::tie(_front, _back);
    while(1) {
        {
            std::unique_lock<std::mutex> lock(mutex);
            cv.wait(lock, []{ return !queue.empty(); });
            tuple = queue.front();
            queue.pop();
        }
        staticTransformation(_front, _back);
    }
}

int main(int argc, const char * argv[]) {
    size_t buffer_size = 0;

    auto start = std::chrono::steady_clock::now();

    unsigned concurentThreadsSupported = std::thread::hardware_concurrency();
    unsigned numberOfThreads = std::min(concurentThreadsSupported,8u);

    std::thread threads[8];

    for(unsigned i=0; i<numberOfThreads; i++)
        threads[i] = std::thread( staticTransformationWorker,i);

    buffer_size = helpers::getFileSize();

    char* buffer = static_cast<char*>(malloc(buffer_size));

    // Read the file and close
    FILE *file = fopen(INPUTFILENAME,"r");
    fread(buffer, sizeof(char), buffer_size, file);
    fclose(file);

    char* sequence_start = 0;  // start of the plain base codes, excl. header and first nl
    char* sequence_end = 0;    // end of the plain base codes, eccl. the last nl

    auto tuple = std::tie(sequence_start, sequence_end);

    for(char* cursor = buffer; cursor < buffer+buffer_size; cursor++) {
        if(*cursor == '>') {   // found header of a sequence
            if(sequence_end > sequence_start) { // Sequenz fertig gecodet
                {
                    std::lock_guard<std::mutex> lock(mutex);
                    queue.push(tuple);
                }
                cv.notify_all();
            }
            while(*cursor != '\n') cursor++;  //  skip the header
            cursor++;  // skip the newline at the end of the header
            sequence_start = cursor;   // store address of first letter code
        }
        if(*cursor == '\n') { // potential end of a sequence
            sequence_end = cursor-1;
            if(*(cursor+1)!='>') // when not a sequence header
                while(*cursor != '\n') cursor++;  //  skip the line
        }
    }

    // Reverse complement transformation in the while-loop is only triggered when a new sequence starts
    // the last transformation won´t be triggered. We can run ist synchronous now, since this task has nothing to do.

    staticTransformation(sequence_start,sequence_end);

    file = fopen(OUTPUTFILENAME,"w");
    fwrite(buffer, sizeof(char), buffer_size, file);
    fclose(file);

    for(unsigned i=0; i<numberOfThreads; i++)
        threads[i].detach();

    free(buffer);

    auto duration = std::chrono::duration_cast< std::chrono::milliseconds>
    (std::chrono::steady_clock::now() - start);
    std::cout << "\n\nZeit: " << duration.count() << " ms" << std::endl;

    return 0;
}

a typical input file would be something like: (much larger of course)

>ONE Homo sapiens alu
GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGGGAGGCCGAGGCGGGCGGA
CGGGCGTGGTGGCGCGCGCCTGTAATCCCAGCTACTCGGGAGGCTGAGGCAGGAGAATCG
CTTGAACCCGGGAGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCCAGCCTG
GGCGACAGAGCGAGACTCCG

One annoying thing is, for example, that running it several times directly in a row gets slower and slower results.


fwrite is about 40.5% and fread is 16% of overall execution time.
Here is the testfile I´m using: http://privatfrickler.de/70000000.txt

\$\endgroup\$
  • \$\begingroup\$ Is using SSE allowed to speed it up? \$\endgroup\$ – harold Aug 19 '18 at 20:16
  • \$\begingroup\$ hmm, dunno .... there is no guarantee for a specific processor, but that would be required wouldn´t it? \$\endgroup\$ – DrSvanHay Aug 19 '18 at 20:22
  • \$\begingroup\$ Yes but only a little, it won't run on a mobile phone or router but x86-without-SSE is only found in museums \$\endgroup\$ – harold Aug 19 '18 at 20:30
  • 1
    \$\begingroup\$ One thing I noticed that I would recommend changing is where you start the timer used to profile the runs. Having it started right at the beginning of main is a bad place for it as you are doing a whole bunch of setup with the timer running that may not run in a consistent time, notably reading the file will likely not be the same every time. Same thing for the timer after completion. I would move the profiling closer to the loop that is actually doing the work in order to get more useful/accurate profiling results. \$\endgroup\$ – DJHenjin Aug 20 '18 at 9:22
  • 2
    \$\begingroup\$ What have you profiled? At the moment, it's hard to tell whether your bottleneck is the file I/O or the processing. It makes a difference, because algorithms that are I/O friendly aren't necessarily as parallel or as cache-friendly, for example. \$\endgroup\$ – Toby Speight Aug 20 '18 at 9:34
5
\$\begingroup\$

The main operation is substitution via a small table, which is also what _mm_shuffle_epi8 does. The low 4 bits of the indexes clash though, and I could not find an easy way to remap them. The clash can also be worked around by remapping the A..O (4x range) and P..Z (5x range) ranges separately with a shuffle for each, then conditional-select based on which range the input character was in. The lower case letters can be handled easily by upper casing them first, which without special characters or validation can be done with just a bitwise AND.

I could not find a way to gracefully handle newlines with SSE, but they can be detected and then handled with scalar code. That slows the code down significantly around a newline though, since a single newline can cause a dozen scalar iterations, which are now also slower than usual since every iteration checks whether it might be the start of a newline-free 16 byte block yet.

My compiler couldn't handle that fancy array initialization so I changed it, of course you can change that back, but this is the code as I tested it:

void staticTransformationSSSE3(char* front, char *back) {
    char c[128];
    c['A'] = 'T'; c['a'] = 'T';
    c['C'] = 'G'; c['c'] = 'G';
    c['G'] = 'C'; c['g'] = 'C';
    c['T'] = 'A'; c['t'] = 'A';
    c['U'] = 'A'; c['u'] = 'A';
    c['M'] = 'K'; c['m'] = 'K';
    c['R'] = 'Y'; c['r'] = 'Y';
    c['W'] = 'W'; c['w'] = 'W';
    c['S'] = 'S'; c['s'] = 'S';
    c['Y'] = 'R'; c['y'] = 'R';
    c['K'] = 'M'; c['k'] = 'M';
    c['V'] = 'B'; c['v'] = 'B';
    c['H'] = 'D'; c['h'] = 'D';
    c['D'] = 'H'; c['d'] = 'H';
    c['B'] = 'V'; c['b'] = 'V';
    c['N'] = 'N'; c['n'] = 'N';

    char help;
    __m128i newline = _mm_set1_epi8('\n');
    __m128i lowercase = _mm_set1_epi8(~0x20);
    //                               N    M       K          H    G          D    C    B    A
    __m128i sub4x = _mm_set_epi8(0, 'N', 'K', 0, 'M', 0, 0, 'D', 'C', 0, 0, 'H', 'G', 'V', 'T', 0);
    //                                              Y       W    V    U    T    S    R
    __m128i sub5x = _mm_set_epi8(0, 0, 0, 0, 0, 0, 'R', 0, 'W', 'B', 'A', 'A', 'S', 'Y', 0, 0);
    __m128i reverse = _mm_set_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15);
    while (front + 15 < back - 15) {
        __m128i fdata = _mm_loadu_si128((__m128i*)front);
        __m128i bdata = _mm_loadu_si128((__m128i*)(back - 15));
        // if there are no newlines, process 2x16 chars
        // otherwise do a scalar iteration
        if (_mm_movemask_epi8(_mm_or_si128(_mm_cmpeq_epi8(fdata, newline), 
                                           _mm_cmpeq_epi8(bdata, newline))) == 0) {
            // toUpper
            fdata = _mm_and_si128(fdata, lowercase);
            bdata = _mm_and_si128(bdata, lowercase);
            // substitute 4x and 5x with separate tables
            __m128i subf4x = _mm_shuffle_epi8(sub4x, fdata);
            __m128i subf5x = _mm_shuffle_epi8(sub5x, fdata);
            __m128i maskf = _mm_cmpgt_epi8(fdata, _mm_set1_epi8(0x4F));
            // note: could use _mm_blendv_epi8 with SSE4.1
            __m128i subf = _mm_or_si128(_mm_and_si128(subf5x, maskf), _mm_andnot_si128(maskf, subf4x));
            __m128i fpart = _mm_shuffle_epi8(subf, reverse);

            // substitute 4x and 5x with separate tables
            __m128i subb4x = _mm_shuffle_epi8(sub4x, bdata);
            __m128i subb5x = _mm_shuffle_epi8(sub5x, bdata);
            __m128i maskb = _mm_cmpgt_epi8(bdata, _mm_set1_epi8(0x4F));
            // note: could use _mm_blendv_epi8 with SSE4.1
            __m128i subb = _mm_or_si128(_mm_and_si128(subb5x, maskb), _mm_andnot_si128(maskb, subb4x));
            __m128i bpart = _mm_shuffle_epi8(subb, reverse);

            // write results swapped
            _mm_storeu_si128((__m128i*)front, bpart);
            _mm_storeu_si128((__m128i*)(back - 15), fpart);

            front += 16;
            back -= 16;
        }
        else {
            if (*front == '\n')
                front++;
            if (*back == '\n')
                back--;
            help = c[(*front)];
            *front = c[(*back)];
            *back = help;
            front++;
            back--;
        }
    }
    while (front < back) {
        if (*front == '\n')
            front++;
        if (*back == '\n')
            back--;
        help = c[(*front)];
        *front = c[(*back)];
        *back = help;
        front++;
        back--;
    }
}

By the way you need something like #include <tmmintrin.h> in order to access the SSSE3 intrinsics, and a compiler-dependent compile flag may need to be used when compiler, such as an appropriate -march= or -mssse3 (not needed on MSVC).

The performance is strongly dependent on the number and distribution of newlines, and will also depend on the microarchitecture. For example for 16KB of data without newlines, on my 4770K Haswell the SSSE3 version is about 10 times as fast, but if there is 1 newline per 256 bytes it already goes down to only 6 times as fast.


Newline handling was a bit more important than I estimated, it can be improved by always doing the SIMD character substitution and reverse but then only writing back up to the newline. Then the next newline test will step over that newline immediately, instead of some scalar iterations having to happen first.

Finding the first newline from both sides can be done by bit-scanning the comparison mask from both sides, for which you could use these functions:

// same semantics as __builtin_ctz,
// but also works on MSVC
int countTrailingZeros(uint32_t x) {
#ifdef _MSC_VER
    unsigned long res;
    _BitScanForward(&res, x);
    return res;
#else
    return __builtin_ctz(x);
#endif
}

// same semantics as __builtin_clz,
// but also works on MSVC
int countLeadingZeros(uint32_t x) {
#ifdef _MSC_VER
    unsigned long res;
    _BitScanReverse(&res, x);
    return res ^ 31;
#else
    return __builtin_clz(x);
#endif
}

The conditional writeback can be done with the same type of conditional-select also used to merge the results of the 4x and 5x range character substitutions, with the mask depending on the number of characters up to the newline. There are various ways to calculate that mask, the easiest is looking it up in a small table. The reverse mask, for conditional writeback at the back, can be found in the same table by using 16 - index and then using the mask with inverted interpretation (conditional-select with reversed operands).

void staticTransformationSSSE3(char* front, char *back) {
    char c[128];
    c['A'] = 'T'; c['a'] = 'T';
    c['C'] = 'G'; c['c'] = 'G';
    c['G'] = 'C'; c['g'] = 'C';
    c['T'] = 'A'; c['t'] = 'A';
    c['U'] = 'A'; c['u'] = 'A';
    c['M'] = 'K'; c['m'] = 'K';
    c['R'] = 'Y'; c['r'] = 'Y';
    c['W'] = 'W'; c['w'] = 'W';
    c['S'] = 'S'; c['s'] = 'S';
    c['Y'] = 'R'; c['y'] = 'R';
    c['K'] = 'M'; c['k'] = 'M';
    c['V'] = 'B'; c['v'] = 'B';
    c['H'] = 'D'; c['h'] = 'D';
    c['D'] = 'H'; c['d'] = 'H';
    c['B'] = 'V'; c['b'] = 'V';
    c['N'] = 'N'; c['n'] = 'N';

    __m128i lengthMaskLUT[16];
    {
        __m128i m = _mm_setzero_si128();
        __m128i ones = _mm_set1_epi8(0xFF);
        for (size_t i = 0; i < 16; i++) {
            lengthMaskLUT[i] = m;
            m = _mm_alignr_epi8(m, ones, 15);
        }
    }

    char help;
    __m128i newline = _mm_set1_epi8('\n');
    __m128i lowercase = _mm_set1_epi8(~0x20);
    //                               N    M       K             H    G          D    C    B    A
    __m128i sub4x = _mm_set_epi8(0, 'N', 'K', 0, 'M', '\n', 0, 'D', 'C', 0, 0, 'H', 'G', 'V', 'T', 0);
    //                                              Y       W    V    U    T    S    R
    __m128i sub5x = _mm_set_epi8(0, 0, 0, 0, 0, 0, 'R', 0, 'W', 'B', 'A', 'A', 'S', 'Y', 0, 0);
    __m128i reverse = _mm_set_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15);
    while (front + 15 < back - 15) {
        // ignore newline
        if (*front == '\n') {
            front++;
            continue;
        }
        if (*back == '\n') {
            back--;
            continue;
        }
        __m128i fdata = _mm_loadu_si128((__m128i*)front);
        __m128i bdata = _mm_loadu_si128((__m128i*)(back - 15));
        uint32_t newlineMaskF = _mm_movemask_epi8(_mm_cmpeq_epi8(fdata, newline));
        uint32_t newlineMaskB = _mm_movemask_epi8(_mm_cmpeq_epi8(bdata, newline));
        // toUpper
        fdata = _mm_and_si128(fdata, lowercase);
        bdata = _mm_and_si128(bdata, lowercase);
        // substitute 4x and 5x with separate tables
        __m128i subf4x = _mm_shuffle_epi8(sub4x, fdata);
        __m128i subf5x = _mm_shuffle_epi8(sub5x, fdata);
        __m128i maskf = _mm_cmpgt_epi8(fdata, _mm_set1_epi8(0x4F));
        // note: could use _mm_blendv_epi8 with SSE4.1
        __m128i subf = _mm_or_si128(_mm_and_si128(subf5x, maskf), _mm_andnot_si128(maskf, subf4x));
        __m128i fpart = _mm_shuffle_epi8(subf, reverse);

        // substitute 4x and 5x with separate tables
        __m128i subb4x = _mm_shuffle_epi8(sub4x, bdata);
        __m128i subb5x = _mm_shuffle_epi8(sub5x, bdata);
        __m128i maskb = _mm_cmpgt_epi8(bdata, _mm_set1_epi8(0x4F));
        // note: could use _mm_blendv_epi8 with SSE4.1
        __m128i subb = _mm_or_si128(_mm_and_si128(subb5x, maskb), _mm_andnot_si128(maskb, subb4x));
        __m128i bpart = _mm_shuffle_epi8(subb, reverse);

        // if there is no newline in the block, write the whole thing
        // otherwise write up to the newline
        if ((newlineMaskF | newlineMaskB) == 0) {
            _mm_storeu_si128((__m128i*)front, bpart);
            _mm_storeu_si128((__m128i*)(back - 15), fpart);
            front += 16;
            back -= 16;
        }
        else {
            // find newline-free prefix
            // must have non-zero length otherwise the "ignore newline"
            // step would have triggered
            int prefixLength = std::min(
                countTrailingZeros(newlineMaskF | 0x10000),
                countLeadingZeros((newlineMaskB << 16) | 0x8000));
            // merge with mask, don't use _mm_maskmoveu_si128 since it's an NT store
            __m128i mergedAtF = _mm_or_si128(
                _mm_and_si128(bpart, lengthMaskLUT[prefixLength]),
                _mm_andnot_si128(lengthMaskLUT[prefixLength], fdata));
            _mm_storeu_si128((__m128i*)front, mergedAtF);
            // reverse mask is the complement of the mask at the negated index
            __m128i mergedAtB = _mm_or_si128(
                _mm_and_si128(bdata, lengthMaskLUT[16 - prefixLength]),
                _mm_andnot_si128(lengthMaskLUT[16 - prefixLength], fpart));
            _mm_storeu_si128((__m128i*)(back - 15), mergedAtB);
            front += prefixLength;
            back -= prefixLength;
        }
    }
    while (front < back) {
        if (*front == '\n')
            front++;
        if (*back == '\n')
            back--;
        help = c[(*front)];
        *front = c[(*back)];
        *back = help;
        front++;
        back--;
    }
}
\$\endgroup\$
  • 1
    \$\begingroup\$ +1 Thank you so much that you did so much work! I will need days to understand that code but I verified it with my testfiles and it works! The speed increase is about 80ms per 100 MB. Thats not as much as could have been expected but there is a newline every 63 bytes. But anyway the 80ms break the sound barrier and I would have been happy to just save 5 or 10 ms :-) Thank you! \$\endgroup\$ – DrSvanHay Aug 20 '18 at 7:42
  • \$\begingroup\$ @DrSvanHay OK I improved that a bit for more newlines, I have some trickier ideas but I haven't worked those out yet \$\endgroup\$ – harold Aug 20 '18 at 10:39
  • \$\begingroup\$ Yeah, thats another 20ms :-) and another 5 days to understand that :-) I really appreciate your efford! Thank you! \$\endgroup\$ – DrSvanHay Aug 20 '18 at 13:14
4
\$\begingroup\$

You're using the C++ versions of the C headers (e.g. <cstdio>). That's a good thing, as it puts the identifiers into the std namespace where they can be unambiguous.

It seems that your implementation also defines these names in the global namespace, but you cannot depend on that. So you need to qualify the names you use: std::size_t, std::malloc (why not new[]?), std::FILE, std::fopen, std::fread, std::fclose, std:free (and any that I've missed).

Oh, and sizeof (char) is always 1, since sizeof reports in units of char.

Have you considered using OpenMP rather than creating and destroying threads directly? It may also help with vectorizing your code, too.

\$\endgroup\$
  • \$\begingroup\$ +1 Stupid mistake leaving out the std::. Very good point! Thank you. The sizeof is there just to remind me, that the parameter is a size :-) \$\endgroup\$ – DrSvanHay Aug 20 '18 at 13:17
3
\$\begingroup\$

I see a number of things that may help you improve your program.

Include relevant #include files

The code uses std::mutex which is defined in <mutex> and std::condition_variable which is in <condition_variable> but the code doesn't include those headers. It should.

Avoid compiler extensions

The current code includes this code:

constexpr char complement[] = {
    ['A']='T', ['a']='T',
    // etc.
    ['N']='N', ['n']='N'};

The problem with that is that designated initialization is, at this time a C++20 proposal and even if adopted, would not support this syntax.

Add error handling

At the moment, the getFileSize function returns -1 if the call to stat fails. This results in malloc being called with -1 which is translated, on my 64-bit machine, into 18446744073709551615 and so, of course, malloc fails. Unfortunately, the return value of malloc is also not checked for error and this causes a crash. Similarly, fread and fopen can fail; their return values should be checked.

Fix the bugs

The main code includes this check:

if(*(cursor+1)!='>') // when not a sequence header
    while(*cursor != '\n') cursor++;  //  skip the line

However, because cursor goes all the way to the end of the buffer, this will inevitably access one byte beyond the end of the buffer which is a bug. Further, there is no bound on the inner loop which is also suspect.

Don't detach threads

This code deploys multiple threads and then detaches them all. This is not a good idea generally, and is definitely a problem here. Here's how it's used within the code:

for(unsigned i=0; i<numberOfThreads; i++)
    threads[i].detach();

free(buffer);

auto duration = std::chrono::duration_cast< std::chrono::milliseconds>
(std::chrono::steady_clock::now() - start);

This is exactly why detach should only be used in rare circumstances! Each of the threads was launched with a pair of pointers into buffer. After all threads are detached, the buffer is destroyed, so those pointers that the threads are possibly still using are no longer valid. This is a serious error! Better is to restructure so that the threads end when they are finished with their data and then use join to wait for all threads to complete. It also means that with the code as it is, you're only timing the amount of time to launch the threads and have no way to tell if they've completed or not.

Encapsulate in an object

Right now the code to handle file structure details is partly in some free-standing functions and partly in main. It would be much neater and easier to understand if something like a DNASequence object were used. That way all of the ugly details about how the file is processed would be neatly encapsulated within a single object implementation, and even details about how many threads to launch could reasonably be delegated.

Consider using standard library functions

There are two essential components to the task of converting a DNA sequence to its reverse complement: reversing and complementing. These could be done in either order. Using objects, here's how I'd write write a main loop:

int main(int argc, char *argv[]) {
    if (argc != 3) {
        std::cerr << "Usage: revdna inputfile outputfile\n";
        return 1;
    }
    std::ifstream in{argv[1]};
    std::ofstream out{argv[2]};
    DNASequence dna;
    while (in >> dna) {
        out << dna.reverse_complement();
    }
}

Note that there is no mention of tasks or threads here, since the reading of the file is sequential anyway. Given the simplicity of the transformation required, I would expect that the program would largely be I/O bound and not computationally bound, so speeding the transformation is unlikely to benefit the overall performance as much as simply having faster hardware (e.g. solid-state drives).

I'd probably write the reverse_complement function using std::reverse and std::toupper with a custom locale.

An example

I reworked the code using most of the suggestions above except using std::reverse and a custom locale. Here is the result:

#include <iostream>
#include <fstream>
#include <string>
#include <limits>
#include <algorithm>

class DNASequence {
public:
    friend std::istream& operator>>(std::istream& in, DNASequence& dna);
    friend std::ostream& operator<<(std::ostream& out, const DNASequence& dna);
    DNASequence& reverse_complement();
    static constexpr char delim = '>';
private:
    static char complement(char ch);
    std::string header;
    std::string dna; 
};

std::istream& operator>>(std::istream& in, DNASequence& dna) {
    // look for a header
    in.ignore(std::numeric_limits<std::streamsize>::max(), dna.delim); 
    std::getline(in, dna.header);
    // fetch the DNA sequence
    std::getline(in, dna.dna, dna.delim);
    in.putback(dna.delim);
    return in;
}

std::ostream& operator<<(std::ostream& out, const DNASequence& dna) {
    return out << dna.delim << dna.header << '\n' << dna.dna;
}

char DNASequence::complement(char ch) {
    switch (ch) {
        case 'A':
        case 'a':
            ch = 'T';
            break;
        case 'C':
        case 'c':
            ch ='G';
            break;
        case 'G':
        case 'g':
            ch ='C'; 
            break;
        case 'T':
        case 't':
            ch ='A';
            break;
        case 'U':
        case 'u':
            ch ='A';
            break;
        case 'M':
        case 'm':
            ch ='K';
            break;
        case 'R':
        case 'r':
            ch ='Y';
            break;
        case 'W':
        case 'w':
            ch ='W';
            break;
        case 'S':
        case 's':
            ch ='S';
            break;
        case 'Y':
        case 'y':
            ch ='R';
            break;
        case 'K':
        case 'k':
            ch ='M';
            break;
        case 'V':
        case 'v':
            ch ='B';
            break;
        case 'H':
        case 'h':
            ch ='D';
            break;
        case 'D':
        case 'd':
            ch ='H';
            break;
        case 'B':
        case 'b':
            ch ='V';
            break;
        case 'N':
        case 'n':
            ch ='N';
            break;
        default:
            ch = '?';  // this is an error!
    }
    return ch;
}

DNASequence& DNASequence::reverse_complement() {
    auto front = dna.begin();
    auto back = front + dna.size() - 1;
    while (back > front) {
        if (*front == '\n') {
            ++front;
        } else if (*back == '\n') {
            --back;
        } else {
            *back = complement(*back);
            *front = complement(*front);
            std::swap(*back, *front);
            ++front;
            --back;
        }
    }
    return *this;
}

int main(int argc, char *argv[]) {
    if (argc != 3) {
        std::cerr << "Usage: revdna inputfile outputfile\n";
        return 1;
    }
    std::ifstream in{argv[1]};
    std::ofstream out{argv[2]};
    DNASequence dna;
    while (in >> dna) {
        out << dna.reverse_complement();
    }
}

Results

I downloaded a FASTA format file for Monodelphis domestica chromosome 1 which is a 724M file. Running this program on my machine, (a 64-bit Linux box), I get the following times:

real    0m5.988s
user    0m0.802s
sys     0m1.123s

I think you'll find that when actually measured, iostreams, when implemented and used in a rational way, are not necessarily any slower than the old FILE I/O. See this question and as always, don't assume; measure!

More measurements and notes

In the interests of further investigating timing implications, I've done some more investigation. Specifically, I've concatenated the FASTA file above with that of Drosophila melanogaster chromosome 3L to obtain a file 787253613 byte long (751Mib). I altered the original as little as possible to create a single-threaded version that compiles cleanly with gcc. That version is here:

#include <iostream>
#include <chrono>
#include <cstdio>
#include <sys/stat.h>

static char complement[256];

void comp_init() {
    complement['A']='T'; complement['a']='T';
    complement['C']='G'; complement['c']='G';
    complement['G']='C'; complement['g']='C';
    complement['T']='A'; complement['t']='A';
    complement['U']='A'; complement['u']='A';
    complement['M']='K'; complement['m']='K';
    complement['R']='Y'; complement['r']='Y';
    complement['W']='W'; complement['w']='W';
    complement['S']='S'; complement['s']='S';
    complement['Y']='R'; complement['y']='R';
    complement['K']='M'; complement['k']='M';
    complement['V']='B'; complement['v']='B';
    complement['H']='D'; complement['h']='D';
    complement['D']='H'; complement['d']='H';
    complement['B']='V'; complement['b']='V';
    complement['N']='N'; complement['n']='N';
}

namespace helpers {
    size_t getFileSize(const char *filename)
    {
        struct stat stat_buf;
        int rc = stat(filename, &stat_buf);
        return rc == 0 ? stat_buf.st_size : -1;
    }
}

static inline void staticTransformation(char* front, char *back) {
    char help;
    while(front < back) {
        if(*front=='\n') {
            front++;
        }
        if(*back=='\n') {
           back--;
        }
        // swap after complement
        help = complement[(*front)];
        *front = complement[(*back)];
        *back = help;
        front++;
        back--;
    }
}

int main(int argc, const char * argv[]) {
    if (argc != 3) {
        std::cerr << "Usage: revdna inputfile outputfile\n";
        return 1;
    }
    comp_init();
    size_t buffer_size = 0;
    auto start = std::chrono::steady_clock::now();
    buffer_size = helpers::getFileSize(argv[1]);
    char* buffer = static_cast<char*>(malloc(buffer_size));
    // Read the file and close
    FILE *file = fopen(argv[1],"r");
    fread(buffer, sizeof(char), buffer_size, file);
    fclose(file);
    char* sequence_start = 0;  // start of the plain base codes, excl. header and first nl
    char* sequence_end = 0;    // end of the plain base codes, eccl. the last nl

    for(char* cursor = buffer; cursor < buffer+buffer_size; cursor++) {
        if(*cursor == '>') {   // found header of a sequence
            if(sequence_end > sequence_start) { // Sequenz fertig gecodet
                staticTransformation(sequence_start, sequence_end);
            }
            while(*cursor != '\n') cursor++;  //  skip the header
            cursor++;  // skip the newline at the end of the header
            sequence_start = cursor;   // store address of first letter code
        }
        if(*cursor == '\n') { // potential end of a sequence
            sequence_end = cursor-1;
            if(*(cursor+1)!='>') // when not a sequence header
                while(*cursor != '\n') cursor++;  //  skip the line
        }
    }
    staticTransformation(sequence_start,sequence_end);
    file = fopen(argv[2],"w");
    fwrite(buffer, sizeof(char), buffer_size, file);
    fclose(file);
    free(buffer);
    auto duration = std::chrono::duration_cast< std::chrono::milliseconds>(
            std::chrono::steady_clock::now() - start);
    std::cout << "\n\nZeit: " << duration.count() << " ms\n";
}

g++ vs clang++ and iostream vs FILE

Compiled with -O3 using g++ 8.1.1 on a 64-bit Linux (Redhat distribution) the iostream version of the code:

real    0m6.846s
user    0m0.862s
sys     0m1.095s

These numbers are fairly consistent across multiple runs.

Using clang++ 6.0.1

real    0m6.776s
user    0m1.156s
sys     0m1.188s

With the FILE-based code above with g++:

real    0m6.362s
user    0m0.780s
sys     0m0.842s

And with clang++:

real    0m6.298s
user    0m0.790s
sys     0m0.832s

All were compiled with the same flags (-O3 and -lstdc++) so were using the same optimization level and library. Also, all versions of the program produced identical outputs which were then verified by running again on the output and comparing the doubly-transformed file to the original (they were, as expected, identical.) As measured by perf, the iostream versions varied by about 3%, and had around 446,000 page faults, while the FILE versions varied by around 1% and had around 116 page faults.

Differences which may account for timing differences are that the iostream version reads in one block at a time (8K by default on this machine) and uses a std::string which is resized multiple times during the read, while the FILE version reads the whole file at once (unless it can't and then it fails with a segfault) and performs the reversing in-place on a single buffer.

Significant speed-up

Curiously, if I insert this line within main (in any version):

std::remove(argv[2]);

The real time drops from around 6.5 to around 1.7 seconds.

\$\endgroup\$
  • \$\begingroup\$ While I appreciate most of your comments - thank you! Using ifstream and ofstream is by far off limits. They are so slow that any optimization attempt would be futile. The primary goal is speed, as I mentioned above. \$\endgroup\$ – DrSvanHay Aug 20 '18 at 19:48
  • 1
    \$\begingroup\$ @DrSvanHay, how did you determine that <iostream> classes are slow enough to affect the speed of your program? Without measurement, even educated guesses are just guesses... \$\endgroup\$ – Toby Speight Aug 20 '18 at 20:12
  • \$\begingroup\$ the version presented here is the last iteration of a whole set of experiments. I tested several other options using c and i tested all possible versions using iostream since i did not believe that it would be so slow. I´ll post a number in a few minutes. \$\endgroup\$ – DrSvanHay Aug 20 '18 at 20:15
  • \$\begingroup\$ Your program takes 1900ms on my machine for a 70MB input file and my code needs 730ms. Thats quite a difference. Especially when you consider that your code ignores the last sequence in a file \$\endgroup\$ – DrSvanHay Aug 20 '18 at 20:36
  • \$\begingroup\$ honestly: Adding e.g. an error-checking if clause to the filesize determinations does cost maybe 100ns, so thats a really good idea. Checking that the tasks are finished before writing the file is a good idea, thats another 100ns, so no problem. Your "Fix the bugs": The code gets even faster by fixing that, great! And object encapsulation? I think that comes for free or it costs some ns, so yes! Great Idea. I really appreciate your suggestions! But I spent days for researching the file io and thats measured in and out, no assumptions made. \$\endgroup\$ – DrSvanHay Aug 20 '18 at 20:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.