#pragma once
/******
 * Purpose: This program is design to do mapping SOLiD or Illumina Reads,
 * The input is a set of short reads file and a Genome_Indes_TableQ
 * The mapping result is output to a file
 *  Author: Yangho Chen
 */
#include "LongReadsSet.h"
#include "PairedReadsSet.h"
#include "ReadsMappingStats.h"
#include "MappingResult.h"
#include "Genome_Index_TableQ.h"
#include "ParameterList.h"
#include "AlignmentsQ.h"
#include "ReadInBitsSet.h"
#include "ReadInBits.h"
#include "ShortReadUtil.h"
#include "ColorSpaceRead.h"
#include "GenomeNTdata.h"
#include "chromosomeNTdata.h"
#include "stdafx.h"
#include <time.h>
#include <vector>
#ifdef WIN32
#include <direct.h>
#endif

// Use OpenMP is gcc version is later than 4.2
#ifdef __GNUC__
#ifdef __GNUC_PATCHLEVEL__
#define __GNUC_VERSION__ (__GNUC__ * 10000 \
                            + __GNUC_MINOR__ * 100 \
                            + __GNUC_PATCHLEVEL__)
#else
#define __GNUC_VERSION__ (__GNUC__ * 10000 \
                            + __GNUC_MINOR__ * 100)
# endif

#if __GNUC_VERSION__ >= 40200
#include <omp.h>
#endif
#else
#ifdef _MSC_VER
#if _MSC_VER >= 2000
#include <omp.h>
#endif
#endif
#endif

// Macro for Parallelization with OpenMP
#ifndef _OPENMP
#define __OPENMP_FOR_PARALLEL__(SHARP_SIGN, openmp_flag) \
{LOG_INFO("Info %d: Use single CPU because OpenMP is not available.\n", CONFIG_LOG);}
#else
#define __OPENMP_FOR_PARALLEL__(SHARP_SIGN, openmp_flag) {\
    int numberOfCPUs = omp_get_num_procs();\
    std::cout << numberOfCPUs << " CPUs." << BLANK_LINE << "\n";\
    SHARP_SIGN##openmp_flag omp parallel for\
}
#endif

// For both counter for chromosome and tolerated substitution error
const int ALIGNMENT_RESULT_FILE_BUFFER_SIZE = 1000000;

int parallelMappingLongReads(vector<string>& readSetsList,\
                             CGenome_Index_TableQ& indexTable, MappingOpts P);

int parallelMapping(vector<string>& readSetsList,\
                    CGenome_Index_TableQ& indexTable, MappingOpts P);

/*
 * This class maps Illumina or SOliD short reads to reference genome.
 * It switches different ways to do read mapping, according to different options.
 */
class CReadsMapping: public CReadsMappingStats
{
public:
    CReadsMapping(void);
    CReadsMapping(MappingOpts P);
    ~CReadsMapping(void);
    MappingOpts opt;
    char cOutputFormat;
    int mapReadsSets(const char* ReadsSetsList, CGenome_Index_TableQ& table,  bool bDiscardReadsWN = true);
    int mapReads(CReadInBitsSet& readsSet, const CGenome_Index_TableQ& table);
    int mapLongReads(CLongReadsSet& pairedReadSet, const CGenome_Index_TableQ& table);
    int queryALongReadInColors(CReadInBits& r1stHalf, CReadInBits& r2ndHalf, const CGenome_Index_TableQ& table, CAlignmentsQ& aQue) const;
    int queryALongReadInBase(CReadInBits& r1stHalf, CReadInBits& r2ndHalf, const CGenome_Index_TableQ& table, CAlignmentsQ& aQue) const;
protected:
    int printMapInfo(CReadInBitsSet& readsSet, int seedOpt);
    int setUpIO4Aligment(const char* Solexafile, const CGenome_Index_TableQ& table);
    int tearDownIO4Aligment(void);
    FileOutputBuffer*  AlignResult;
    FileOutputBuffer*  AmbiguousReads;
    FileOutputBuffer*  BadReads;
    FileOutputBuffer*  MissReads;
    time_t start, end;

    const static unsigned int CHECK_POINTS = 1000000;
    int dealMappedRead(const CGenome_Index_TableQ& table, CAlignmentsQ& aQue);
    int dealMappedLongRead(const CGenome_Index_TableQ& table, CAlignmentsQ& aQue, CMappingResult& m);
    int dealMissedRead(CMappingResult& m);
    int dealMissedRead(bool bMapReadInColors, const char* readName, CReadInBits r, const char* qs = NULL);
    int dealAmbiguousRead(CMappingResult& m);
    int dealAmbiguousRead(bool bMapReadInColors, const char* readName, CReadInBits r, const char* qs = NULL);
    int printSingleEndReads(CMappingResult& m);
    int printLogFile(const char* inputFile);
    inline void printCheckPointInfo(int i);
protected:
    void initialization(void);
    inline void printRead(FileOutputBuffer* FileBuf, CMappingResult& m);
    inline void printRead(FileOutputBuffer* FileBuf, bool bMapReadInColors, const char* readName, CReadInBits r, const char* qs);
    inline void getLongBaseReadInfo(CReadInBitsSet& readsSet1stHalf, CReadInBitsSet& readsSet2ndHalf, int readId,\
                                    CReadInBits& r1stHalf, CReadInBits& r2ndHalf, CMappingResult& m);
    inline void getLongColorReadInfo(CReadInBitsSet& readSet1stHalf, CReadInBitsSet& readSet2ndHalf, int readId, \
                                     CReadInBits& r1stHalf, CReadInBits& r2ndHalf, CMappingResult& m);
    string getMappingFileN(const char* caReadsSetName, const CGenome_Index_TableQ& table);
    unsigned int checkPairedReadSetSize(CReadInBitsSet& firstHalfSet, CReadInBitsSet& SecondHalfSet);
};

inline void CReadsMapping::getLongBaseReadInfo\
(CReadInBitsSet& readSet1stHalf, CReadInBitsSet& readSet2ndHalf, int readId, \
 CReadInBits& r1stHalf, CReadInBits& r2ndHalf, \
 CMappingResult& m)
{
    m.uiReadLength = this->opt.readLength; // Check the read length
    // Get read tag
    if ((int)readSet1stHalf.pReadsID->size() > readId) {
        strcpy(m.QNAME, readSet1stHalf.pReadsID->at(readId).id);
    } else {
        sprintf(m.QNAME, "Read_%d", this->iReadCounter + readId);
    }
    // Get read sequence
    //unsigned int halfReadLength = this->opt.readLength / 2;
    unsigned int halfReadLength = this->opt.anchorLength;
    unsigned int secondHalfStart = this->opt.readLength - this->opt.anchorLength;
    r1stHalf.decode(m.caRead);
    r2ndHalf.decode(&(m.caRead[secondHalfStart]));
    // Get the quality score
    if (readSet1stHalf.pQualScores != NULL) {
        const char* qS1stHalf = readSet1stHalf.pQualScores->qScores((unsigned int)readId);
        trQScores(halfReadLength, 0, qS1stHalf, m.rawScores);
        const char* qS2ndHalf = readSet2ndHalf.pQualScores->qScores((unsigned int)readId);
        trQScores(halfReadLength, 0, qS2ndHalf, &m.rawScores[secondHalfStart]);
        trQScores(m.uiReadLength, SolexaScoreEncodingShift, m.rawScores, m.QScores);
    }
    if (this->cOutputFormat == 's') {
        m.getReverseReadandQual();
    }
}

inline void CReadsMapping::getLongColorReadInfo\
(CReadInBitsSet& readSet1stHalf, CReadInBitsSet& readSet2ndHalf, int readId, \
 CReadInBits& r1stHalf, CReadInBits& r2ndHalf, \
 CMappingResult& m)
{
    m.uiReadLength = this->opt.readLength; // Check the read length
    // Get read tag
    if ((int)readSet1stHalf.pReadsID->size() > readId) {
        strcpy(m.QNAME, readSet1stHalf.pReadsID->at(readId).id);
    } else {
        sprintf(m.QNAME, "Read_%d", this->iReadCounter + readId);
    }
    // Get read sequence
    unsigned int halfReadLength = this->opt.anchorLength;
    unsigned int secondHalfStart = this->opt.readLength - this->opt.anchorLength;
    decodeColors(m.caRead, r1stHalf);
    decodePureColors(&(m.caRead[secondHalfStart]), r2ndHalf);
    // Get the quality score for SOLiD
    /* TODO fix the SOLiD score (base score and color space score*/
    if (readSet1stHalf.pQualScores != NULL) {
        const char* qS1stHalf = readSet1stHalf.pQualScores->qScores((unsigned int)readId);
        trQScores(halfReadLength, 0, qS1stHalf, m.rawScores);
        const char* qS2ndHalf = readSet2ndHalf.pQualScores->qScores((unsigned int)readId);
        trQScores(halfReadLength, 0, qS2ndHalf, &m.rawScores[secondHalfStart]);
        trQScores(m.uiReadLength, Phred_SCALE_QUAL_SHIFT, m.rawScores, m.QScores);
    }
    // TODO for sam Format, one need reversed quality and Seq
    if (this->cOutputFormat == 's') {
        m.getReverseReadandQual();
    }
}

inline void CReadsMapping::printCheckPointInfo(int readNo)
{
    if (readNo % this->CHECK_POINTS == 0) {
        printf("Mapping no %u reads.\r", this->iReadCounter + readNo);
        fflush(stdout);
    }
}

// The following function fill different part of CMappingResult
inline
bool getSingleMappingIndex(CGenomeNTdata& pgenomeNT, CAlignmentsQ& aQue, int mappingIndex, CMappingResult &m)
{
    m.uiDiff = aQue.asdiff[mappingIndex];
    m.MultipleMappedNo = aQue.load;
    m.strand = (mappingIndex >= (int)aQue.ForwardAlignmentLoad) ? '-' : '+' ;

    m.uiGlobalMappedPos = aQue.aiHitIndex[mappingIndex];
    m.uiRefId = pgenomeNT.genomeIndex2chrID(m.uiGlobalMappedPos);
    m.uiPOS = pgenomeNT.genomeLocusID2chrIndex(m.uiGlobalMappedPos);
    ChrIndex2GeneName& geneVec = pgenomeNT.paChromosomes[m.uiRefId]->geneVec;
    if (geneVec.table.size() > 0) {
        CGene g = geneVec.query(m.uiPOS);
        strcpy(m.RNAME, g.name.c_str());
        if(!g.isValid) { // this is for handeling mapping exception for mapping 2 NULL_REGION
            char* chrName = pgenomeNT.paChromosomes[m.uiRefId]->caInputFileName;
            char* strBuf = &m.RNAME[strlen(g.name.c_str())];
            myStrCpy(strBuf, chrName, FILENAME_MAX/2);
            return(false);
        }
        // g.startIndex is the translated index of m.uiPOS. Not the start index of gene
        m.uiPOS = g.startIndex;
    } else {
        sprintf(m.RNAME, "%d", m.uiRefId);
    }
    return(true);
}
void getQscores4Solexa(CAlignmentsQ& aQue, CMappingResult& m, bool samFormat);
void getReadQscores4Solexa(CAlignmentsQ& aQue, CMappingResult& m, bool samFormat);
void getSingleMappingSeqAndQ4SOLiD\
(const CGenome_Index_TableQ& table, CAlignmentsQ& aQue, CMappingResult& m, bool samFormat);
void getSingleMappingSeq4Solexa(const CGenome_Index_TableQ& table, CMappingResult& m, bool samFormat);
void getLongMappingInfo(const CGenome_Index_TableQ& table, CAlignmentsQ& aQue, bool samFormat,\
                        unsigned int mappingId, CMappingResult& m);
void getSingleMappingInfo(const CGenome_Index_TableQ& table, CAlignmentsQ& aQue,\
                          unsigned int mappingId, CMappingResult& m, bool samFormat);

inline void printSamHeader(FileOutputBuffer* AlignResult, vector<CGene>&refs, const char* RG, const char* CL)
{
    sprintf(AlignResult->caBufp, "@HD\tVN:0.1.5c\tSO:queryname\n");
    AlignResult->UpdateSize();
    for (vector<CGene>::iterator it = refs.begin(); it != refs.end(); it++ ) {
        sprintf(AlignResult->caBufp, "@SQ\tSN:%s\tLN:%u\n", it->name.c_str(), it->startIndex);
        AlignResult->UpdateSize(); // startIndex is actually the length of reference.
    }
    sprintf(AlignResult->caBufp, "%s\n@PG\tID:PerM\tVN:0.4.0\tCL:\"%s\"\n", RG, CL);
    AlignResult->UpdateSize();
}

inline const char* getLongRefSeq(const CGenome_Index_TableQ& table, CMappingResult& m, bool bNoRef)
{
    if(bNoRef) {
        m.caRef[0] = '\0';
    } else {
        unsigned int secondHalfStart = m.uiReadLength - table.uiRead_Length;
        CReadInBits ref1stHalf = table.pgenomeNTInBits->getSubstringInBits\
                                 (m.uiGlobalMappedPos, table.uiRead_Length);
        CReadInBits ref2ndHalf = table.pgenomeNTInBits->getSubstringInBits\
                                 (m.uiGlobalMappedPos + secondHalfStart, table.uiRead_Length);
        ref1stHalf.decode(m.caRef);
        ref2ndHalf.decode(&(m.caRef[secondHalfStart]));
        if(m.strand == '-') {
            reverseComplementKmer(m.caRef); // reverse complement reference
        }
    }
    return(m.caRef);
}

bool wrongIndex(const CReadInBitsSet& readsSet, const CGenome_Index_TableQ& table);