#include "PairedReadsMapping.h"
// Given two mapped paired read set list and the index table,
// this function maps reads in parallel
int parallelMappingPairedReads(vector<string>& readSetsList1, vector<string>& readSetsList2,\
                               CGenome_Index_TableQ& indexTable, MappingOpts P)
{
    unsigned int readStartIndex = P.truncatedReadPrefix;// ignore the first few base pair
    P.bDiscardReadWithN = false;
    ASSERT_TRUE(((int)(readSetsList1.size()) == (int)(readSetsList2.size())), "Read sets are not in paired");
    P.clearOutputFileName(readSetsList1.size() > 1);
    int i;

#ifdef _OPENMP
    int numberOfCPUs = omp_get_num_procs();
    LOG_INFO("\nInfo %d: %d CPUs detected. %s.\n",\
             INFO_LOG, numberOfCPUs, BLANK_LINE);
    #pragma omp parallel for
#endif
//    __OPENMP_FOR_PARALLEL__(#, pragma)
    for (i = 0; i < min((int)readSetsList1.size(), (int)readSetsList2.size()); i++) {
        CPairedReadsMapping mapping(P);
        const char* readSetName1 = (readSetsList1.at(i)).c_str();
        const char* readSetName2 = (readSetsList2.at(i)).c_str();
        if (checkFileExist(readSetName1) && checkFileExist(readSetName2)) {
            CReadInBitsSet readSet1\
            (readSetName1, P.readsFileFormat, readStartIndex, indexTable.uiRead_Length, P.allowedNumOfNinRead);
            CReadInBitsSet readSet2\
            (readSetName2, P.readsFileFormat, readStartIndex, indexTable.uiRead_Length, P.allowedNumOfNinRead);
            if (P.bIgnoreQS) {
                readSet1.ignoreQScores();
                readSet2.ignoreQScores();
            }
            TIME_INFO(mapping.mapPairedReadsInPairedFiles(readSet1, readSet2, indexTable), "Mapping takes");
        }
    }
    return(0);
}

// Given the list of paired read sets in the single file
// format and the index table, this function maps reads int parallel.
int parallelMappingPairedReads(vector<string>& readSetsList,\
                               CGenome_Index_TableQ& indexTable, MappingOpts P)
{
    P.bDiscardReadWithN = false;
    P.clearOutputFileName(readSetsList.size() > 1);
    int i;
#ifdef _OPENMP
    int numberOfCPUs = omp_get_num_procs();
    LOG_INFO("\nInfo %d: %d CPUs detected. %s.\n",\
             INFO_LOG, numberOfCPUs, BLANK_LINE);
    #pragma omp parallel for
#endif
//__OPENMP_FOR_PARALLEL__(#pragma)#ifdef _OPENMP
    for (i = 0; i < (int)readSetsList.size(); i++) {
        CPairedReadsMapping mapping(P);
        const char* readSetName = (readSetsList.at(i)).c_str();
        if (checkFileExist(readSetName)) {
            // assume the paired end reads is in the format
            // that concatenate paired read in 5'-3' and 3' to 5'.
            bool in5to3cat3to5Format = true;
            int allowedNumOfNinRead = P.readLength; // Don't throw reads with N in the paired-end mapping
            CPairedReadsSet pairedReadSet\
            (readSetName, P.readsFileFormat, indexTable.uiRead_Length * 2, in5to3cat3to5Format, allowedNumOfNinRead);
            if (P.bIgnoreQS) {
                pairedReadSet.ignoreQScores();
            }
            TIME_INFO(mapping.mapPairedReads(*pairedReadSet.F_Reads, *pairedReadSet.R_Reads,\
                                             indexTable), "Mapping takes");
        }
    }
    return(0);
}

int parallelMappingPairedLongReads(vector<string>& readSetsList1, vector<string>& readSetsList2,\
                                   CGenome_Index_TableQ& indexTable, MappingOpts P)
{
    unsigned int readStartIndex = P.truncatedReadPrefix;// ignore the first few base pair
    P.bDiscardReadWithN = false;
    ASSERT_TRUE(((int)(readSetsList1.size()) == (int)(readSetsList2.size())), "Read sets are not in paired");
    P.clearOutputFileName(readSetsList1.size() > 1);
    int i;

#ifdef _OPENMP
    int numberOfCPUs = omp_get_num_procs();
    LOG_INFO("\nInfo %d: %d CPUs detected. %s.\n",\
             INFO_LOG, numberOfCPUs, BLANK_LINE);
    #pragma omp parallel for
#endif
//    __OPENMP_FOR_PARALLEL__(#, pragma)
    for (i = 0; i < min((int)readSetsList1.size(), (int)readSetsList2.size()); i++) {
        CPairedReadsMapping mapping(P);
        const char* readSetName1 = (readSetsList1.at(i)).c_str();
        const char* readSetName2 = (readSetsList2.at(i)).c_str();
        if (checkFileExist(readSetName1) && checkFileExist(readSetName2)) {
            int allowedNumOfNinRead = P.readLength; // Don't throw reads with N in the paired-end mapping
            CLongReadsSet readSet1(readSetName1, P.readsFileFormat, P.readLength, allowedNumOfNinRead, readStartIndex);
            CLongReadsSet readSet2(readSetName2, P.readsFileFormat, P.readLength, allowedNumOfNinRead, readStartIndex);

            if (P.bIgnoreQS) {
                readSet1.ignoreQScores();
                readSet2.ignoreQScores();
            }
            TIME_INFO(mapping.mapPairedLongReadsInBases(readSet1, readSet2, indexTable), "Mapping takes");
        }
    }
    return(0);
}

CBestPairedMapping::CBestPairedMapping(void)
{
    this->validMappingNo = 0;
    this->bestMappingNo = 0;
    this->minDiff = MAX_READ_LENGTH;
}

CBestPairedMapping::~CBestPairedMapping(void)
{
}

inline void CBestPairedMapping::update\
(CMappingResult &m1, CMappingResult &m2, bool excludeAmbigousRead)
{
    unsigned int diff = m1.uiDiff + m2.uiDiff;
    if (diff < this->minDiff) {
        this->bm1 = m1;
        this->bm2 = m2;
        minDiff = m1.uiDiff + m2.uiDiff;
        this->bestMappingNo = 1;
    } else if (diff == minDiff) {
        this->bestMappingNo++;
    }
    this->validMappingNo++;
}


CPairedReadsMapping::CPairedReadsMapping(void)
{
    this->initialization();
}

CPairedReadsMapping::~CPairedReadsMapping(void)
{
}

CPairedReadsMapping::CPairedReadsMapping(const MappingOpts P): CReadsMapping( P )
{
    this->initialization();
}

void CPairedReadsMapping::initialization(void)
{
    this->noOfPairsInRange = 0;
    this->noOfSingle1stEndMapped = 0;
    this->noOfSingle2ndEndMapped = 0;
    this->noOfAmbiguousPairs = 0;
    this->noOfPairsSepLess = 0;
    this->noOfPairsSepMore = 0;
    this->noOfPairsSepMoreAndLess = 0;
    this->noOfExpMappedPairedStrand = 0;
}
/*
int CPairedReadsMapping::mapPairedReadsInASingleFile\
(CPairedReadsSet& readSet, CGenome_Index_TableQ& table)
{
    CReadInBitsSet& readSet1 = *readSet.F_Reads;
    CReadInBitsSet& readSet2 = *readSet.R_Reads;
    cout << "Start mapping " << readSet1.InputFile << " and " << readSet2.InputFile << endl;
    getReadsFileFormat(readSet.InputFile, opt.readsFileFormat);
    if (this->setUpIO4Aligment(readSet.InputFile, table) != 0) {
        LOG_INFO("\nInfo %d: Fail to setup I/O files.", ERROR_LOG);
        return(1);
    }
    while (readSet.get_next_capacity_reads_pairs_from_single_file()) {
        mapPairedReads(readSet1, readSet2, table);
    }
    this->printMappedPairStats(cout, readSet1, table.uiSubDiffThreshold * 2);
    this->tearDownIO4Aligment();
    return(0);
}
*/
int CPairedReadsMapping::mapPairedReadsInPairedFiles\
(CReadInBitsSet& readSet1, CReadInBitsSet& readSet2, CGenome_Index_TableQ& table)
{
    cout << "Start mapping " << readSet1.InputFile << " and " << readSet2.InputFile << endl;
    getReadsFileFormat(readSet1.InputFile, opt.readsFileFormat);
    if (this->setUpIO4Aligment(readSet1.InputFile, table) != 0) {
        LOG_INFO("\nInfo %d: Fail to setup I/O files.", ERROR_LOG);
    }
    while (readSet1.get_next_capacity_reads(BUFFERED_READS_SIZE, opt.readtag_delimiter) &&
            readSet2.get_next_capacity_reads(BUFFERED_READS_SIZE, opt.readtag_delimiter)) {
        mapPairedReads(readSet1, readSet2, table);
    }
    this->printMappedPairStats(cout, readSet1, table.uiSubDiffThreshold * 2);
    this->tearDownIO4Aligment();
    return(0);
}

int CPairedReadsMapping::mapPairedReads(CReadInBitsSet& readSet1, CReadInBitsSet& readSet2,
                                        CGenome_Index_TableQ& table)
{
    vector<CReadInBits>::iterator it1 = readSet1.pReadsSet->begin();
    vector<CReadInBits>::iterator it2 = readSet2.pReadsSet->begin();
    for (int i = 0; it1 != readSet1.pReadsSet->end()
            && it2 != readSet2.pReadsSet->end(); it1++, it2++, i++) {
        alignmentsQ[0].read = *it1;
        alignmentsQ[1].read = *it2;
        readSet1.get_read_id(i, alignmentsQ[0].tag);
        readSet2.get_read_id(i, alignmentsQ[1].tag);
        alignmentsQ[0].qualityScores = readSet1.getQScoresPtr(i);
        alignmentsQ[1].qualityScores = readSet2.getQScoresPtr(i);

        if (table.bMapReadInColors) {
            table.queryReadColors(*it1, alignmentsQ[0], true, true);
            table.queryReadColors(*it1, alignmentsQ[0], false, false);
            table.queryReadColors(*it2, alignmentsQ[1], true, true);
            table.queryReadColors(*it2, alignmentsQ[1], false, false);
        } else {
            table.queryReadBases(*it1, alignmentsQ[0], true, true);
            table.queryReadBases(*it1, alignmentsQ[0], false, false);
            table.queryReadBases(*it2, alignmentsQ[1], true, true);
            table.queryReadBases(*it2, alignmentsQ[1], false, false);
        }
        this->dealMappedPairedReads(table);
    }
    return(0);
}

inline bool isExpPairedMappedStrand(bool firstEndFirst, char end1Strand, char end2Strand, bool bSOLiD)
{
    if (bSOLiD) {
        if (firstEndFirst && (end1Strand == '+') && (end2Strand == '+')) {
            return(true);
        } else if (!firstEndFirst && (end1Strand == '-') && (end2Strand == '-')) {
            return(true);
        }
    } else {
        if (firstEndFirst && (end1Strand == '+') && (end2Strand == '-')) {
            return(true);
        } else if (!firstEndFirst && (end1Strand == '-') && (end2Strand == '+')) {
            return(true);
        }
    }
    return(false);
}

inline bool isValidPaired(CMappingResult& m1, CMappingResult& m2, MappingOpts& opts)
{
    bool allowSameStrand = !(opts.frOnly);
    bool allowDiffStrand = !(opts.ffOnly);
    if (strcmp(m1.RNAME, m2.RNAME) == 0) {
        if(m1.strand == m2.strand) {
            return(allowSameStrand);
        } else {
            return(allowDiffStrand);
        }
    }
    return(false);
}

// TODO becareful about the big separation that overflow to negative value
inline int getSep(unsigned int uiPos1, unsigned int uiPos2, bool expM2gtM1)
{
    int range = 0;
    if (uiPos1 < uiPos2) {
        range = (int)(uiPos2 - uiPos1);
        if(!expM2gtM1) {
            range *= -1;
        }
    } else {
        range = (int)(uiPos1 - uiPos2);
        if(expM2gtM1) {
            range *= -1;
        }
    }
    return(range);
}


inline int set_ISIZE(CMappingResult &m1, CMappingResult &m2, int readLength)
{
    bool expM2gtM1 = (m1.strand = '+');
    int sep = getSep(m1.uiPOS, m2.uiPOS, expM2gtM1);
    m2.ISIZE = m1.ISIZE = sep + readLength;
    return(sep);
}

int CPairedReadsMapping::printValidMappedPair(const CGenome_Index_TableQ& table, CMappingResult& m1, CMappingResult& m2, int validMappedPairNo)
{
    for (unsigned int i = 0; i < this->alignmentsQ[0].load; i++) {
        for (unsigned int j = 0; j < this->alignmentsQ[1].load; j++) {
            getSingleMappingIndex(*table.pgenomeNT, this->alignmentsQ[0], i, m1);
            getSingleMappingIndex(*table.pgenomeNT, this->alignmentsQ[1], j, m2);
            // Require reads mapped to the same ref sequence
            if (isValidPaired(m1, m2, opt)) {
                int sep =  set_ISIZE(m1, m2, (int)opt.readLength);
                if (opt.disLB <= sep && sep <= opt.disUB) {
                    this->printAMappedPair(table, m1, m2, validMappedPairNo);
                }
            }
        }
    }
    return(validMappedPairNo);
}

int CPairedReadsMapping::printBestMappedPair(const CGenome_Index_TableQ& table, CMappingResult& m1, CMappingResult& m2, int minMismatchNo, int bestMappedPairNo)
{
    for (unsigned int i = 0; i < this->alignmentsQ[0].load; i++) {
        for (unsigned int j = 0; j < this->alignmentsQ[1].load; j++) {
            getSingleMappingIndex(*table.pgenomeNT, this->alignmentsQ[0], i, m1);
            getSingleMappingIndex(*table.pgenomeNT, this->alignmentsQ[1], j, m2);
            // Require reads mapped to the same ref sequence
            if (isValidPaired(m1, m2, opt)) {
                int sep =  set_ISIZE(m1, m2, (int)opt.readLength);
                if (opt.disLB <= sep && sep <= opt.disUB) {
                    if ((int)(m1.uiDiff + m2.uiDiff) == minMismatchNo) {
                        this->printAMappedPair(table, m1, m2, bestMappedPairNo);
                    }
                }
            }
        }
    }
    return(bestMappedPairNo);
}

int CPairedReadsMapping::dealMappedPairedReads(CGenome_Index_TableQ& table)
{
    bool mapSOLiDRead = table.bMapReadInColors;
    bool bNoMapping = true, sepMore = false, sepLess = false, pairedOnExpStrand = false;
    bool bPrintAllMapping = opt.bGetAllAlignments && !opt.bExcludeAmbiguousPaired;
    CBestPairedMapping bestMP;
    CMappingResult m1(this->alignmentsQ[0], opt.readLength);
    CMappingResult m2(this->alignmentsQ[1], opt.readLength);
    bool samFormat = (this->cOutputFormat == 's');
    if (!table.bMapReadInColors) {
        // get mapping info
        /*
        getQscores4Solexa(this->alignmentsQ[0], m1, samFormat);
        getQscores4Solexa(this->alignmentsQ[1], m2, samFormat);
        */
        getReadQscores4Solexa(this->alignmentsQ[0], m1, samFormat);
        getReadQscores4Solexa(this->alignmentsQ[1], m2, samFormat);
    }

    for (unsigned int i = 0; i < this->alignmentsQ[0].load; i++) {
        for (unsigned int j = 0; j < this->alignmentsQ[1].load; j++) {
            getSingleMappingIndex(*table.pgenomeNT, this->alignmentsQ[0], i, m1);
            getSingleMappingIndex(*table.pgenomeNT, this->alignmentsQ[1], j, m2);
            // Require reads mapped to the same ref sequence
            if (isValidPaired(m1, m2, opt)) {
                int sep =  set_ISIZE(m1, m2, (int)opt.readLength);
                if (opt.disLB <= sep && sep <= opt.disUB) {
                    bNoMapping = false;
                    bestMP.update(m1, m2, opt.bExcludeAmbiguousReads);
                    /*
                    if (bPrintAllMapping) {
                    	this->printAMappedPair(table, m1, m2, );
                    }*/
                } else {
                    sepMore |= (sep > opt.disUB);
                    sepLess |= (sep < opt.disLB);
                }
                bool firstEndFirst = m1.uiPOS > m2.uiPOS;
                pairedOnExpStrand |= (isExpPairedMappedStrand(firstEndFirst, m1.strand, m2.strand, mapSOLiDRead));
            }
        }
    }
    if (bNoMapping) {
        dealNoMapping(table, m1, m2);
        this->bookNoMappedKeepPairs(sepMore, sepLess, pairedOnExpStrand);
    } else {
        bool isSamFormat = (this->cOutputFormat == 's');
        this->getPairedRInfo(table, m1, m2, isSamFormat);
        if(bPrintAllMapping) {
            this->printValidMappedPair(table, m1, m2, bestMP.validMappingNo);
        } else {
            this->dealBestMapping(table, bestMP, m1, m2);
        }
        this->bookKeepMappedPairs(bestMP);
    }
    return(0);
}

void CPairedReadsMapping::dealNoMapping(const CGenome_Index_TableQ& table, CMappingResult& m1, CMappingResult& m2)
{
    setSamFlags4OnlyOneEndMapped(m1, true);
    const int endId = 1;
    this->dealMappedSingleRead(table, this->alignmentsQ[0], m1, endId == 1);
    setSamFlags4OnlyOneEndMapped(m2,false);
    this->dealMappedSingleRead(table, this->alignmentsQ[1], m2, endId == 2);
}

int CPairedReadsMapping::dealBestMapping(const CGenome_Index_TableQ& table, CBestPairedMapping& bestMP, CMappingResult& m1, CMappingResult& m2)
{
    if(bestMP.validMappingNo == 1) { // -A -e & 1 mapping
        this->printAMappedPair(table, bestMP.bm1, bestMP.bm2, 1);
    } else if (opt.bExcludeAmbiguousPaired && !opt.bGetAllAlignments && bestMP.bestMappingNo == 1) { // -e 1 best mapping
        this->printAMappedPair(table, bestMP.bm1, bestMP.bm2, 1);
    } else if (!opt.bExcludeAmbiguousPaired && !opt.bGetAllAlignments) { // -B
        this->printBestMappedPair(table, m1, m2, bestMP.minDiff, bestMP.bestMappingNo);
    }
    return(bestMP.bestMappingNo);
}

inline void CPairedReadsMapping::getPairedRInfo(const CGenome_Index_TableQ& table, CMappingResult &m1, CMappingResult &m2, bool samFormat)
{
    if (this->opt.bPrintAlignments) {
        if (table.bMapReadInColors) {
            getSingleMappingSeqAndQ4SOLiD(table, this->alignmentsQ[0], m1, samFormat);
            getSingleMappingSeqAndQ4SOLiD(table, this->alignmentsQ[1], m2, samFormat);
        } else {
            bool bNoRef = samFormat || !opt.bPrintRef4PairedInMapping;
            if(opt.bPrintPairedRQ) {
                // Get quality score
                getQscores4Solexa(this->alignmentsQ[0], m1, samFormat);
                getQscores4Solexa(this->alignmentsQ[1], m2, samFormat);
            }
            if(opt.bMappedLongRead) {
                getLongRefSeq(table, m1, bNoRef);
                getLongRefSeq(table, m2, bNoRef);
            } else {
                getSingleMappingSeq4Solexa(table, m1, bNoRef);
                getSingleMappingSeq4Solexa(table, m2, bNoRef);
            }
        }
    }
}

// F3read and R3read are index in AlignmentQ
void CPairedReadsMapping::printAMappedPair\
(const CGenome_Index_TableQ& table, CMappingResult &m1, CMappingResult &m2, int noPairedLoc)
{
    if (this->opt.bPrintAlignments) {
        set_ISIZE(m1, m2, (int)opt.readLength);
        bool samFormat = (this->cOutputFormat == 's');
        // getPairedRInfo(table, m1, m2, samFormat);
        if (samFormat) { // sam format
            printAPairedMappingInSam(this->AlignResult, m1, m2);
        } else {
            // string category = getCategory(m1.strand,m2.strand, m1.ISIZE, opt.disLB, opt.disUB);
            printAPairedMappingInPerM(this->AlignResult, m1, m2, noPairedLoc, opt.bPrintNM);
        }
    }
}

int CPairedReadsMapping::dealMappedSingleRead\
(const CGenome_Index_TableQ& table, CAlignmentsQ &Q, CMappingResult &m, bool bFirstEnd)
{
    bool samFormat = (this->cOutputFormat == 's');
    if (Q.load > 0 && samFormat) {
        for (unsigned int i = 0; i < Q.load && i < Q.iMaxCapacity; i++) {
            if (this->opt.bPrintAlignments) {
                getSingleMappingInfo(table, Q, i, m, samFormat);
                this->printSingleEndReads(m);
            }
        }
    } else if (this->opt.bPrintUnMappedReads) {
        // check if the seq is not decoded for SOLiD read?
        // check if code has never reach here?
        this->dealMissedRead(m);
    }
    if (Q.load > 0) {
        if (bFirstEnd) {
            this->noOfSingle1stEndMapped++;
        } else {
            this->noOfSingle2ndEndMapped++;
        }
    }
    return(this->iMissReadCounter++);
}

void CPairedReadsMapping::printMappedPairStats\
(ostream& out, CReadInBitsSet& readSet, unsigned int uiSubThreshold)
{
    string readSetName = getBasename(readSet.InputFile);
    out << '\n';
    out << readSetName.c_str() << ", #Pairs, " << readSet.uiNo_of_Reads << ", ";
    out << "#Mapped Pairs, " << this->noOfPairsInRange << ", ";
    out << "#Multi-mapped Pairs, " << this->noOfAmbiguousPairs << "\n";
    out << readSetName << ", ";
    out << "#Pairs sep more, " << this->noOfPairsSepMore << ", ";
    out << "#Pairs sep less, " << this->noOfPairsSepLess << ", ";
    out << "#Pairs sep more and less, " << this->noOfPairsSepMoreAndLess << endl;
    out << "#Pairs on exp strands " << this->noOfExpMappedPairedStrand << endl;
    out << readSetName << ", ";
    out << "#Mapped single 1st end, " << this->noOfSingle1stEndMapped << endl;
    out << "#Mapped single 2nd end, " << this->noOfSingle2ndEndMapped << endl;
    out << readSetName << ", ";
    unsigned int i;
    for (i = 0; i <= uiSubThreshold; i++) {
        out << "Sub" << i << ", " << iMapDiffCount[i] << ", " ;
    }
    out << endl;
}
unsigned int CPairedReadsMapping::getPairedReadSetSize
(CReadInBitsSet& setA1, CReadInBitsSet& setA2, CReadInBitsSet& setB1, CReadInBitsSet& setB2)
{
    int setAsize = this->checkPairedReadSetSize(setA1, setA2);
    int setBsize = this->checkPairedReadSetSize(setB1, setB2);
    return(min(setAsize, setBsize));
}

int CPairedReadsMapping::dealMappedLongPairedRead
(CAlignmentsQ& q1, CAlignmentsQ& q2, CMappingResult& m1, CMappingResult& m2,
 const CGenome_Index_TableQ& table)
{
    bool bNoMapping = true, sepMore = false, sepLess = false, pairedOnExpStrand = false;
    bool bPrintAllMapping = opt.bGetAllAlignments && !opt.bExcludeAmbiguousPaired;
    CBestPairedMapping bestMP;

    // bool samFormat = (this->cOutputFormat == 's');
    for (unsigned int i = 0; i < q1.load; i++) {
        for (unsigned int j = 0; j < q2.load; j++) {
            // get mapping info
            getSingleMappingIndex(*table.pgenomeNT, q1, i, m1);
            getSingleMappingIndex(*table.pgenomeNT, q2, j, m2);

            if (isValidPaired(m1, m2, opt)) {
                int sep =  set_ISIZE(m1, m2, (int)opt.readLength);
                if (opt.disLB <= sep && sep <= opt.disUB) {
                    bNoMapping = false;
                    bestMP.update(m1, m2, opt.bExcludeAmbiguousPaired);
                    /*
                    if (bPrintAllMapping) {
                        this->printAMappedPair(table, m1, m2, );
                    }*/
                } else {
                    sepMore |= (sep > opt.disUB);
                    sepLess |= (sep < opt.disLB);
                }
                bool firstEndFirst = m1.uiPOS > m2.uiPOS;
                bool mapSOLiDRead = false;
                pairedOnExpStrand |= (isExpPairedMappedStrand(firstEndFirst, m1.strand, m2.strand, mapSOLiDRead));
            }
        }
    }
    if (opt.bExcludeAmbiguousReads && bestMP.bestMappingNo == 1) {
        this->printAMappedPair(table, bestMP.bm1, bestMP.bm2, bestMP.bestMappingNo);
    }
    if (bNoMapping) {
        dealNoMapping(table, m1, m2);
        this->bookNoMappedKeepPairs(sepMore, sepLess, pairedOnExpStrand);
    } else {
        // TODO get quality score, read sequence?
        if(bPrintAllMapping) {
            this->printValidMappedPair(table, m1, m2, bestMP.validMappingNo);
        } else {
            this->dealBestMapping(table, bestMP, m1, m2);
        }
        this->bookKeepMappedPairs(bestMP);
        // Book keep mapping
    }
    return(0);
}

int CPairedReadsMapping::mapPairedLongReadsInBases
(CLongReadsSet& longReadSet1, CLongReadsSet& longReadSet2, const CGenome_Index_TableQ& table)
{
    unsigned int uiReadLength = this->opt.readLength;
    CAlignmentsQ& aQue1 = this->alignmentsQ[0];
    CAlignmentsQ& aQue2 = this->alignmentsQ[1];
    CReadInBitsSet& readSetA1stHalf = *(longReadSet1.F_Reads);
    CReadInBitsSet& readSetA2ndHalf = *(longReadSet1.R_Reads);
    CReadInBitsSet& readSetB1stHalf = *(longReadSet2.F_Reads);
    CReadInBitsSet& readSetB2ndHalf = *(longReadSet2.R_Reads);
    const char* readSetName = readSetA1stHalf.InputFile;
    // Flag that set the alignment is ambiguous or not
    getReadsFileFormat(readSetName, opt.readsFileFormat);
    string seedStr = seedSymbol(table.chosenSeedId);
    printf("Mapping %s (%u-bp reads) with %s seed.\n", \
           readSetName, uiReadLength, seedStr.c_str());

    this->initializeStatsCounter();
    if (this->setUpIO4Aligment(readSetName, table) != 0) {
        LOG_INFO("\nInfo %d: Fail to setup I/O files.", ERROR_LOG);
        return(1);
    }
    // alignmentsQ[0].setQueue_All_Best_OneFlag('A');
    while (get_next_capacity_long_paired_reads(longReadSet1, longReadSet2) > 0) {
        int bufferedReadNo =
            this->getPairedReadSetSize(readSetA1stHalf, readSetA2ndHalf, readSetB1stHalf, readSetB2ndHalf);
        vector<CReadInBits>::iterator itA1, itA2, itB1, itB2;
        itA1 = readSetA1stHalf.pReadsSet->begin();
        itA2 = readSetA2ndHalf.pReadsSet->begin();
        itB1 = readSetB1stHalf.pReadsSet->begin();
        itB2 = readSetB2ndHalf.pReadsSet->begin();
        for (int i = 0; i < bufferedReadNo; i++, itA1++, itA2++, itB1++, itB2++) {
            this->printCheckPointInfo(i);
            CMappingResult m1, m2;
            this->getLongBaseReadInfo(readSetA1stHalf, readSetA2ndHalf, i, *itA1, *itA2, m1);
            this->getLongBaseReadInfo(readSetB1stHalf, readSetB2ndHalf, i, *itB1, *itB2, m2);
            this->queryALongReadInBase(*itA1, *itA2, table, aQue1);
            this->queryALongReadInBase(*itB1, *itB2, table, aQue2);

            // statistics and output
            if (aQue1.load > 0 && aQue2.load > 0) {
                this->dealMappedLongPairedRead(aQue1, aQue2, m1, m2, table);
            } else if (this->opt.bPrintUnMappedReads) {
                // TODO to be done
                // dealMissedPairedRead(m1, m2);
            }
        }
        this->iReadCounter += bufferedReadNo;
    }
    this->tearDownIO4Aligment();
    this->iBadReadCounter = longReadSet1.uiNo_of_Bad_Reads + longReadSet2.uiNo_of_Bad_Reads;
    this->printLogFile(readSetName);
    return(0);
}