<commit_before>/*++

Module Name:

    SmarterPairedEndAligner.cpp

Abstract:

    A more sophisticated paired-end aligner.

Authors:

    <NAME>, February, 2012

Environment:

    User mode service.

Revision History:

--*/

#include "stdafx.h"
#include <math.h>
#include "SmarterPairedEndAligner.h"
#include "LandauVishkin.h"

using namespace std;


#ifdef TRACE_PAIRED_ALIGNER
#define TRACE printf
#else
#define TRACE(...) {}
#endif


SmarterPairedEndAligner::SmarterPairedEndAligner(
         GenomeIndex  *index_,
         unsigned      maxReadSize_,
         unsigned      confDiff_,
         unsigned      maxHits_,
         unsigned      maxK_,
         unsigned      maxSeeds_,
         unsigned      minSpacing_,
         unsigned      maxSpacing_,
         unsigned      adaptiveConfDiffThreshold_)
 :   index(index_), maxReadSize(maxReadSize_), confDiff(confDiff_), maxHits(maxHits_),
     maxK(maxK_), maxSeeds(maxSeeds_), minSpacing(minSpacing_), maxSpacing(maxSpacing_),
     adaptiveConfDiffThreshold(adaptiveConfDiffThreshold_), seedLen(index_->getSeedLength()),
     lv(2 * MAX_BUCKETS)
{
    // Initialize the bucket data structures.
    buckets = new Bucket[MAX_BUCKETS];
    for (int r = 0; r < 2; r++) {
        for (int rc = 0; rc < 2; rc++) {
            bucketLocations[r][rc].reserve(MAX_BUCKETS);
            bucketTable[r][rc].reserve(2 * MAX_BUCKETS);
        }
    }
    candidates.reserve(MAX_BUCKETS);
    
    // Initialize the complements array.
    memset(complement, 0, sizeof(complement));
    complement['A'] = 'T';
    complement['C'] = 'G';
    complement['G'] = 'C';
    complement['T'] = 'A';
    complement['N'] = 'N';
    complement['n'] = 'n';
    
    // Initialize the seed offsets used when we wrap around in testing seeds.
    vector<bool> used(seedLen + 1, false);
    used[0] = true;
    used[seedLen] = true;      // Sentinel to make the calculation work nicely
    wrapOffset[0] = 0;
    for (int i = 1; i < seedLen; i++) {
        // For the next wrap offset, pick the offset that is most distant from used offsets.
        int bestDist = 0;
        int bestP = 0;
        for (int p = 1; p < seedLen; p++) {
            int dist = 0;
            while (!used[p - dist] && !used[p + dist]) {
                dist++;
            }
            if (dist > bestDist) {
                bestDist = dist;
                bestP = p;
            }
        }
        wrapOffset[i] = bestP;
        used[bestP] = 1;
    }
#ifdef TRACE_PAIRED_ALIGNER
    printf("Wrap offsets:");
    for (int i = 0; i < seedLen; i++) {
        printf(" %d", wrapOffset[i]);
    }
    printf("\n");
#endif

    // Create single-end aligners.
    singleAligner = new BaseAligner(index, confDiff + 1, maxHits, maxK, maxReadSize,
                                    maxSeeds, 1000000 /* lvLimit */, adaptiveConfDiffThreshold);
    mateAligner = new BaseAligner(index, confDiff, maxHits, maxK, maxReadSize,
                                  maxSeeds, 1000000 /* lvLimit */, adaptiveConfDiffThreshold, &lv);
}


SmarterPairedEndAligner::~SmarterPairedEndAligner()
{
    delete[] buckets;
}


void SmarterPairedEndAligner::align(Read *read0, Read *read1, PairedAlignmentResult *result)
{
    Read *reads[2] = {read0, read1};
    int numNotFound = 0;
    int numCertainlyNotFound = 0;
    int numIgnoredMulti = 0;  // MultiHits where all seeds returned too many hits
    int numSingleWithNotFound = 0;

    clearState();
    
    result->status[0] = NotFound;
    result->status[1] = NotFound;

    result->location[0] = 0xFFFFFFFF;
    result->location[1] = 0xFFFFFFFF;

    unsigned bestScore[2] = {INFINITE_SCORE, INFINITE_SCORE};
    unsigned singleLoc[2] = {0xFFFFFFFF, 0xFFFFFFFF};
    bool singleIsRC[2] = {false, false};
    bool certainlyNotFound[2] = {false, false};
    bool bestScoreCertain[2] = {false, false};
    
    TRACE("Aligning read pair %.*s:\n%.*s\n%.*s\n",
            reads[0]->getIdLength(), reads[0]->getId(),
            reads[0]->getDataLength(), reads[0]->getData(),
            reads[1]->getDataLength(), reads[1]->getData());

    if (read0->getDataLength() < 50 && read1->getDataLength() < 50) {
        TRACE("Reads are both too short -- returning");
        return;
    }

#define SINGLE_FIRST
#ifdef SINGLE_FIRST
    // First try aligning each end alone and looking for its mate nearby if it was confident.
    // TODO: Guess which end will be faster to search for first based on it seeds?
    // TODO: We'll want to search for chimeras here when the second read is not found.
    for (int r = 0; r < 2; r++) {
        unsigned loc0, loc1;
        bool rc0, rc1;
        int score0, score1;
        AlignmentResult status0, status1;
        singleAligner->setReadId(r);
        status0 = singleAligner->AlignRead(reads[r], &loc0, &rc0, &score0);
        bestScore[r] = (score0 <= (int) (maxK + confDiff ? score0 : INFINITE_SCORE));
        bestScoreCertain[r] = singleAligner->checkedAllSeeds();
        TRACE("Read %d returned %s (%d) at loc %u-%d\n", r, AlignmentResultToString(status0), score0, loc0, rc0);
        if (isOneLocation(status0)) {
            // Let's just search for the other read nearby.
            if (reads[1-r]->getDataLength() < 50) {
                TRACE("Read %d was found but %d is too short, so returning SingleHit/NotFound\n", r, 1-r);
                result->status[r] = SingleHit;
                result->location[r] = loc0;
                result->isRC[r] = rc0;
                result->status[1-r] = NotFound;
                return;
            }
            mateAligner->setMaxK(maxK - score0 + 1);
            mateAligner->setReadId(1-r);
            status1 = mateAligner->AlignRead(reads[1-r], &loc1, &rc1, &score1, maxSpacing, loc0, !rc0);
            TRACE("Mate %d returned %s at loc %u-%d\n", 1-r, AlignmentResultToString(status1), loc1, rc1);
            if (/*status1 != MultipleHits &&*/ score0 + score1 <= (int)maxK) {
                result->status[r] = SingleHit;
                result->location[r] = loc0;
                result->isRC[r] = rc0;
                result->status[1-r] = isOneLocation(status1) ? SingleHit : status1;
                result->location[1-r] = loc1;
                result->isRC[1-r] = rc1;
                return;
            } else if(status1 == NotFound) {
                // We found read r at one location and didn't find the mate nearby. Let's remember because
                // if the mate is not found anywhere else, we can just return a location for read r.
                singleLoc[r] = loc0;
                singleIsRC[r] = rc0;
                numSingleWithNotFound++;
            }
        } else if (status0 == NotFound) {
            numNotFound++;
            certainlyNotFound[r] = singleAligner->checkedAllSeeds();
            numCertainlyNotFound += certainlyNotFound[r] ? 1 : 0;
        } else if (status0 == MultipleHits && loc0 == 0xFFFFFFFF) {
            numIgnoredMulti++;
            result->status[r] = MultipleHits;
        }
    }
    
    if (numNotFound + numIgnoredMulti == 2) {
        // We couldn't find either read even as MultipleHits, so there's no way we'll do the pair.
        return;
    }

    if (numCertainlyNotFound == 1) {
        // One read was NotFound and we checked all its seeds, so there's no way we will find it as
        // either a bucket element or bucket mate in alignTogether. Let's just return the result for
        // the other one.
        TRACE("One read is certainly not found: %d %d %d, %d %d %d\n",
                certainlyNotFound[0], singleLoc[0], bestScore[0],
                certainlyNotFound[1], singleLoc[1], bestScore[1]);
        for (int r = 0; r < 2; r++) {
            if (certainlyNotFound[r]) {
                result->status[r] = NotFound;
            } else {
                result->status[r] = (singleLoc[r] != 0xFFFFFFFF && bestScore[r] <= 0.6 * maxK) ? SingleHit : NotFound;
                result->location[r] = singleLoc[r];
                result->isRC[r] = singleIsRC[r];
            }
        }
        return;
    }

    if (numSingleWithNotFound == 1 && numNotFound == 1) {
        // One read was Single and the other was NotFound, neither near its mate nor elsewhere.
        TRACE("One read was single and its mate was NotFound either nearby or by itself\n");
        for (int r = 0; r < 2; r++) {
            if (singleLoc[r] == 0xFFFFFFFF) {
                result->status[r] = NotFound;
            } else {
                result->status[r] = (bestScore[r] <= 0.6 * maxK ? SingleHit : NotFound);
                result->location[r] = (bestScore[r] <= 0.6 * maxK ? singleLoc[r] : 0xFFFFFFFF);
                result->isRC[r] = singleIsRC[r];
            }
        }
        return;
    }

    if (read0->getDataLength() < 50 || read1->getDataLength() < 50) {
        TRACE("Can't go further because one read is too short");
        return;
    }

    if (bestScore[0] != INFINITE_SCORE && bestScore[1] != INFINITE_SCORE &&
            bestScore[0] + bestScore[1] > 1.5 * maxK) {
        TRACE("Best scores in each direction add up to more than 1.5 * maxK; returning");
        return;
    }

    int lowerBound[2] = {0, 0};
    if (bestScoreCertain[0]) {
        lowerBound[0] = bestScore[0] < confDiff + 1 ? 0 : bestScore[0];
    }
    if (bestScoreCertain[1]) {
        lowerBound[1] = bestScore[1] < confDiff + 1 ? 0 : bestScore[1];
    }

    if (lowerBound[0] + lowerBound[1] > (int)maxK) {<commit_msg>Bug fix - set result score on early exit from align<commit_after>/*++

Module Name:

    SmarterPairedEndAligner.cpp

Abstract:

    A more sophisticated paired-end aligner.

Authors:

    <NAME>, February, 2012

Environment:

    User mode service.

Revision History:

--*/

#include "stdafx.h"
#include <math.h>
#include "SmarterPairedEndAligner.h"
#include "LandauVishkin.h"

using namespace std;


#ifdef TRACE_PAIRED_ALIGNER
#define TRACE printf
#else
#define TRACE(...) {}
#endif


SmarterPairedEndAligner::SmarterPairedEndAligner(
         GenomeIndex  *index_,
         unsigned      maxReadSize_,
         unsigned      confDiff_,
         unsigned      maxHits_,
         unsigned      maxK_,
         unsigned      maxSeeds_,
         unsigned      minSpacing_,
         unsigned      maxSpacing_,
         unsigned      adaptiveConfDiffThreshold_)
 :   index(index_), maxReadSize(maxReadSize_), confDiff(confDiff_), maxHits(maxHits_),
     maxK(maxK_), maxSeeds(maxSeeds_), minSpacing(minSpacing_), maxSpacing(maxSpacing_),
     adaptiveConfDiffThreshold(adaptiveConfDiffThreshold_), seedLen(index_->getSeedLength()),
     lv(2 * MAX_BUCKETS)
{
    // Initialize the bucket data structures.
    buckets = new Bucket[MAX_BUCKETS];
    for (int r = 0; r < 2; r++) {
        for (int rc = 0; rc < 2; rc++) {
            bucketLocations[r][rc].reserve(MAX_BUCKETS);
            bucketTable[r][rc].reserve(2 * MAX_BUCKETS);
        }
    }
    candidates.reserve(MAX_BUCKETS);
    
    // Initialize the complements array.
    memset(complement, 0, sizeof(complement));
    complement['A'] = 'T';
    complement['C'] = 'G';
    complement['G'] = 'C';
    complement['T'] = 'A';
    complement['N'] = 'N';
    complement['n'] = 'n';
    
    // Initialize the seed offsets used when we wrap around in testing seeds.
    vector<bool> used(seedLen + 1, false);
    used[0] = true;
    used[seedLen] = true;      // Sentinel to make the calculation work nicely
    wrapOffset[0] = 0;
    for (int i = 1; i < seedLen; i++) {
        // For the next wrap offset, pick the offset that is most distant from used offsets.
        int bestDist = 0;
        int bestP = 0;
        for (int p = 1; p < seedLen; p++) {
            int dist = 0;
            while (!used[p - dist] && !used[p + dist]) {
                dist++;
            }
            if (dist > bestDist) {
                bestDist = dist;
                bestP = p;
            }
        }
        wrapOffset[i] = bestP;
        used[bestP] = 1;
    }
#ifdef TRACE_PAIRED_ALIGNER
    printf("Wrap offsets:");
    for (int i = 0; i < seedLen; i++) {
        printf(" %d", wrapOffset[i]);
    }
    printf("\n");
#endif

    // Create single-end aligners.
    singleAligner = new BaseAligner(index, confDiff + 1, maxHits, maxK, maxReadSize,
                                    maxSeeds, 1000000 /* lvLimit */, adaptiveConfDiffThreshold);
    mateAligner = new BaseAligner(index, confDiff, maxHits, maxK, maxReadSize,
                                  maxSeeds, 1000000 /* lvLimit */, adaptiveConfDiffThreshold, &lv);
}


SmarterPairedEndAligner::~SmarterPairedEndAligner()
{
    delete[] buckets;
}


void SmarterPairedEndAligner::align(Read *read0, Read *read1, PairedAlignmentResult *result)
{
    Read *reads[2] = {read0, read1};
    int numNotFound = 0;
    int numCertainlyNotFound = 0;
    int numIgnoredMulti = 0;  // MultiHits where all seeds returned too many hits
    int numSingleWithNotFound = 0;

    clearState();
    
    result->status[0] = NotFound;
    result->status[1] = NotFound;

    result->location[0] = 0xFFFFFFFF;
    result->location[1] = 0xFFFFFFFF;

    unsigned bestScore[2] = {INFINITE_SCORE, INFINITE_SCORE};
    unsigned singleLoc[2] = {0xFFFFFFFF, 0xFFFFFFFF};
    bool singleIsRC[2] = {false, false};
    bool certainlyNotFound[2] = {false, false};
    bool bestScoreCertain[2] = {false, false};
    
    TRACE("Aligning read pair %.*s:\n%.*s\n%.*s\n",
            reads[0]->getIdLength(), reads[0]->getId(),
            reads[0]->getDataLength(), reads[0]->getData(),
            reads[1]->getDataLength(), reads[1]->getData());

    if (read0->getDataLength() < 50 && read1->getDataLength() < 50) {
        TRACE("Reads are both too short -- returning");
        return;
    }

#define SINGLE_FIRST
#ifdef SINGLE_FIRST
    // First try aligning each end alone and looking for its mate nearby if it was confident.
    // TODO: Guess which end will be faster to search for first based on it seeds?
    // TODO: We'll want to search for chimeras here when the second read is not found.
    for (int r = 0; r < 2; r++) {
        unsigned loc0, loc1;
        bool rc0, rc1;
        int score0, score1;
        AlignmentResult status0, status1;
        singleAligner->setReadId(r);
        status0 = singleAligner->AlignRead(reads[r], &loc0, &rc0, &score0);
        bestScore[r] = (score0 <= (int) (maxK + confDiff ? score0 : INFINITE_SCORE));
        bestScoreCertain[r] = singleAligner->checkedAllSeeds();
        TRACE("Read %d returned %s (%d) at loc %u-%d\n", r, AlignmentResultToString(status0), score0, loc0, rc0);
        if (isOneLocation(status0)) {
            // Let's just search for the other read nearby.
            if (reads[1-r]->getDataLength() < 50) {
                TRACE("Read %d was found but %d is too short, so returning SingleHit/NotFound\n", r, 1-r);
                result->status[r] = SingleHit;
                result->location[r] = loc0;
                result->isRC[r] = rc0;
                result->score[r] = score0;
                result->status[1-r] = NotFound;
                return;
            }
            mateAligner->setMaxK(maxK - score0 + 1);
            mateAligner->setReadId(1-r);
            status1 = mateAligner->AlignRead(reads[1-r], &loc1, &rc1, &score1, maxSpacing, loc0, !rc0);
            TRACE("Mate %d returned %s at loc %u-%d\n", 1-r, AlignmentResultToString(status1), loc1, rc1);
            if (/*status1 != MultipleHits &&*/ score0 + score1 <= (int)maxK) {
                result->status[r] = SingleHit;
                result->location[r] = loc0;
                result->isRC[r] = rc0;
                result->score[r] = score0;
                result->status[1-r] = isOneLocation(status1) ? SingleHit : status1;
                result->location[1-r] = loc1;
                result->isRC[1-r] = rc1;
                result->score[1-r] = score1;
                return;
            } else if(status1 == NotFound) {
                // We found read r at one location and didn't find the mate nearby. Let's remember because
                // if the mate is not found anywhere else, we can just return a location for read r.
                singleLoc[r] = loc0;
                singleIsRC[r] = rc0;
                numSingleWithNotFound++;
            }
        } else if (status0 == NotFound) {
            numNotFound++;
            certainlyNotFound[r] = singleAligner->checkedAllSeeds();
            numCertainlyNotFound += certainlyNotFound[r] ? 1 : 0;
        } else if (status0 == MultipleHits && loc0 == 0xFFFFFFFF) {
            numIgnoredMulti++;
            result->status[r] = MultipleHits;
        }
    }
    
    if (numNotFound + numIgnoredMulti == 2) {
        // We couldn't find either read even as MultipleHits, so there's no way we'll do the pair.
        return;
    }

    if (numCertainlyNotFound == 1) {
        // One read was NotFound and we checked all its seeds, so there's no way we will find it as
        // either a bucket element or bucket mate in alignTogether. Let's just return the result for
        // the other one.
        TRACE("One read is certainly not found: %d %d %d, %d %d %d\n",
                certainlyNotFound[0], singleLoc[0], bestScore[0],
                certainlyNotFound[1], singleLoc[1], bestScore[1]);
        for (int r = 0; r < 2; r++) {
            if (certainlyNotFound[r]) {
                result->status[r] = NotFound;
            } else {
                result->status[r] = (singleLoc[r] != 0xFFFFFFFF && bestScore[r] <= 0.6 * maxK) ? SingleHit : NotFound;
                result->location[r] = singleLoc[r];
                result->isRC[r] = singleIsRC[r];
                result->score[r] = bestScore[r];
            }
        }
        return;
    }

    if (numSingleWithNotFound == 1 && numNotFound == 1) {
        // One read was Single and the other was NotFound, neither near its mate nor elsewhere.
        TRACE("One read was single and its mate was NotFound either nearby or by itself\n");
        for (int r = 0; r < 2; r++) {
            if (singleLoc[r] == 0xFFFFFFFF) {
                result->status[r] = NotFound;
            } else {
                result->status[r] = (bestScore[r] <= 0.6 * maxK ? SingleHit : NotFound);
                result->location[r] = (bestScore[r] <= 0.6 * maxK ? singleLoc[r] : 0xFFFFFFFF);
                result->isRC[r] = singleIsRC[r];
                result->score[r] = bestScore[r];
            }
        }
        return;
    }

    if (read0->getDataLength() < 50 || read1->getDataLength() < 50) {
        TRACE("Can't go further because one read is too short");
        return;
    }

    if (bestScore[0] != INFINITE_SCORE && bestScore[1] != INFINITE_SCORE &&
            bestScore[0] + bestScore[1] > 1.5 * maxK) {
        TRACE("Best scores in each direction add up to more than 1.5 * maxK; returning");
        return;
    }

    int lowerBound[2] = {0, 0};
    if (bestScoreCertain[0]) {
        lowerBound[0] = bestScore[0] < confDiff + 1 ? 0 : bestScore[0];
    }
    if (bestScoreCertain[1]) {
        lowerBound[1] = bestScore[1] < confDiff + 1 ? 0 : bestScore[1];
    }

    if (lowerBound[0] + lowerBound[1] > (int)maxK) {
{'fraction_non_alphanumeric': 0.12885709810473556, 'fraction_numerical': 0.019840233905915522, 'mean_word_length': 2.1539601514902027, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 104, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 16, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7230147', 'n_tokens_mistral': 6275, 'n_tokens_neox': 5934, 'n_words': 1607}