// -*- mode: c++; indent-tabs-mode: nil; -*-
//
// Copyright 2009 Illumina, Inc.
//
// This software is covered by the "Illumina Genome Analyzer Software
// License Agreement" and the "Illumina Source Code License Agreement",
// and certain third party copyright/licenses, and any user of this
// source file is bound by the terms therein (see accompanying files
// Illumina_Genome_Analyzer_Software_License_Agreement.pdf and
// Illumina_Source_Code_License_Agreement.pdf and third party
// copyright/license notices).
//
//

/// \file
///
/// \author Chris Saunders
///


#include "alignment/GlobalUtilities.hh"
#include "blt_util/blt_exception.hh"
#include "blt_util/log.hh"
#include "blt_util/parse_util.hh"
#include "blt_util/seq_util.hh"
#include "starling/align_path.hh"

#include "boost/lexical_cast.hpp"

#include <cassert>
#include <cctype>

#include <algorithm>
#include <iostream>
#include <sstream>



enum {
    INDEL_BEGIN='^',
    INDEL_END='$'
};



static
void
unknown_md_error(const char* const md,
                 const char* const mdptr){
    
    std::ostringstream oss;
    oss << "ERROR: can't parse match descriptor string: " << md << "\n"
        << "\tunexpected character: '" << *mdptr << "' at position: " << (mdptr-md+1) << "\n";
    throw blt_exception(oss.str().c_str());
}



static
void
unknown_cigar_error(const char* const cigar,
                    const char* const cptr){
    
    std::ostringstream oss;
    oss << "ERROR: can't parse cigar string: " << cigar << "\n"
        << "\tunexpected character: '" << *cptr << "' at position: " << (cptr-cigar+1) << "\n";
    throw blt_exception(oss.str().c_str());
}



namespace ALIGNPATH {


static
void
apath_push(path_t& apath,
           path_segment& ps,
           const align_t t){
    
    if( (0==ps.length) or (ps.type==t) ) return;
    apath.push_back(ps);
    ps.clear();
}



static
void
export_md_to_apath_impl(const char* md,
                        path_t& apath) {

    using casava::blt_util::parse_unsigned;

    const char* mdptr(md);
    path_segment ps;

    while(*mdptr) {
        if       (isdigit(*mdptr)) {
            apath_push(apath,ps,MATCH);
            const unsigned mlen(parse_unsigned(mdptr));
            ps.length += mlen;
            ps.type = MATCH;

        } else if(is_valid_base(*mdptr)){
            apath_push(apath,ps,MATCH);
            mdptr++;
            ps.length++;
            ps.type = MATCH;

        } else if(*mdptr == INDEL_BEGIN){
            mdptr++; // eat INDEL_BEGIN

            while(*mdptr != INDEL_END){
                if       (isdigit(*mdptr)){
                    apath_push(apath,ps,INSERT);
                    const unsigned mlen(parse_unsigned(mdptr));
                    ps.length=mlen;
                    ps.type=INSERT;
                    
                } else if(is_valid_base(*mdptr)){
                    apath_push(apath,ps,DELETE);
                    mdptr++;
                    ps.length++;
                    ps.type=DELETE;

                } else {
                    unknown_md_error(md,mdptr);
                }
            }

            mdptr++; // eat INDEL_END

        } else {
            unknown_md_error(md,mdptr);
        }
    }

    apath_push(apath,ps,NONE);
}



void
export_md_to_apath(const char* md,
                   const bool is_fwd_strand,
                   path_t& apath,
                   const bool is_edge_deletion_error){

    // to make best use of previous code, we parse the MD in the
    // alignment direction and then orient apath to the forward strand
    // as a second step if required
    //    
    assert(NULL != md);

    apath.clear();
    export_md_to_apath_impl(md,apath);

    unsigned as(apath.size());

    if( ((as>0) and (apath.front().type == DELETE)) or
        ((as>1) and (apath.back().type == DELETE)) ) {
        std::ostringstream oss;
        if(is_edge_deletion_error) {
            oss << "ERROR: ";
        } else {
            oss << "WARNING: ";
        }
        std::string cigar;
        apath_to_cigar(apath,cigar);
        oss << "alignment path: " << cigar << " contains meaningless edge deletion.\n";
        if(is_edge_deletion_error) {
            throw blt_exception(oss.str().c_str());
        } else {
            log_os << oss.str();
            path_t apath2;
            for(unsigned i(0);i<as;++i){
                if(((i==0) or ((i+1)==as)) and
                   apath[i].type == DELETE) continue;
                apath2.push_back(apath[i]);
            }
            apath=apath2;
            as=apath.size();
        }
    }

    if( (not is_fwd_strand) and (as>1) ) {
        std::reverse(apath.begin(),apath.end());
    }
}



void
apath_to_cigar(const path_t& apath,
               std::string& cigar) {

    cigar.clear();

    const unsigned as(apath.size());
    for(unsigned i(0);i<as;++i){
        const path_segment& ps(apath[i]);
        cigar += boost::lexical_cast<std::string>(ps.length);
        cigar.push_back(segment_type_to_cigar_code(ps.type));
    }
}

void
apath_to_export_md(path_t& apath,
                   const char* ref_seq,
                   const char* ref_end,
                   const int32_t ref_pos,
                   const std::string& read_bases,
                   const bool is_fwd_strand,
                   std::string& md) {

    md.clear();

    if(is_fwd_strand) {

        const char* pRead      = read_bases.c_str();
        const char* pReference = ref_seq + ref_pos - 1;
        fwd_apath_to_export_md(apath, ref_seq, pReference, ref_end, pRead, md);

    } else {

        uint32_t numRefBases = 0;
        path_t::const_iterator pCIter;
        for(pCIter = apath.begin(); pCIter != apath.end(); ++pCIter) {
            if((pCIter->type == MATCH) || (pCIter->type == DELETE)) {
                numRefBases += pCIter->length;
            }
        }

        const char* pRead      = read_bases.c_str();
        const char* pReference = ref_seq + ref_pos + numRefBases - 2;
        rev_apath_to_export_md(apath, ref_seq, pReference, ref_end, pRead, md);
    }
}

void
fwd_apath_to_export_md(path_t& apath,
                       const char* ref_begin,
                       const char* ref_bases,
                       const char* ref_end,
                       const char* read_bases,
                       std::string& md) {

    // process the align path
    bool foundUnsupportedCigar = false;
    path_t::const_iterator pCIter;
    for(pCIter = apath.begin(); pCIter != apath.end(); ++pCIter) {

        if(pCIter->type == DELETE) {

            // handle deletion
            md.push_back('^');
            for(uint32_t i = 0; i < pCIter->length; ++i, ++ref_bases) {
                md.push_back(*ref_bases);
            }
            md.push_back('$');

        } else if(pCIter->type == INSERT) {

            // handle insertion
            md.push_back('^');
            md += boost::lexical_cast<std::string>(pCIter->length);
            read_bases += pCIter->length;
            md.push_back('$');

        } else if(pCIter->type == MATCH) {

            // handle match/mismatch
            uint32_t numMatchingBases = 0;
            for(uint32_t i = 0; i < pCIter->length; ++i, ++ref_bases, ++read_bases) {

                // handle circular genome
                if((ref_bases < ref_begin) || (ref_bases > ref_end)) {
                    md.push_back('N');
                    continue;
                }

                if(*ref_bases != *read_bases) {

                    // write the number of preceding matching bases
                    if(numMatchingBases != 0) {
                        md += boost::lexical_cast<std::string>(numMatchingBases);
                        numMatchingBases = 0;
                    }

                    // output the mismatched base
                    md.push_back(*ref_bases);

                } else ++numMatchingBases;
            }

            // write the number of trailing matching bases
            if(numMatchingBases != 0) {
                md += boost::lexical_cast<std::string>(numMatchingBases);
            }

        } else {

            // handle unsupported CIGAR operation
            foundUnsupportedCigar = true;
            break;
        }
    }

    if(foundUnsupportedCigar) md = "UNSUPPORTED";
}

void
rev_apath_to_export_md(path_t& apath,
                       const char* ref_begin,
                       const char* ref_bases,
                       const char* ref_end,
                       const char* read_bases,
                       std::string& md) {

    // process the align path
    bool foundUnsupportedCigar = false;
    path_t::const_reverse_iterator pCRIter;
    for(pCRIter = apath.rbegin(); pCRIter != apath.rend(); ++pCRIter) {

        if(pCRIter->type == DELETE) {

            // handle deletion
            md.push_back('^');
            for(uint32_t i = 0; i < pCRIter->length; ++i, --ref_bases) {
                md.push_back(reverseCharASCII[(uint)*ref_bases]);
            }
            md.push_back('$');

        } else if(pCRIter->type == INSERT) {

            // handle insertion
            md.push_back('^');
            md += boost::lexical_cast<std::string>(pCRIter->length);
            read_bases += pCRIter->length;
            md.push_back('$');

        } else if(pCRIter->type == MATCH) {

            // recreate the the match descriptor for this non-INDEL region
            uint32_t numMatchingBases = 0;
            for(uint32_t i = 0; i < pCRIter->length; ++i, --ref_bases, ++read_bases) {

                // handle circular genome
                if((ref_bases < ref_begin) || (ref_bases > ref_end)) {
                    md.push_back('N');
                    continue;
                }

                const char rcRefBase = reverseCharASCII[(uint)*ref_bases];

                if(rcRefBase != *read_bases) {

                    // write the number of preceding matching bases
                    if(numMatchingBases != 0) {
                        md += boost::lexical_cast<std::string>(numMatchingBases);
                        numMatchingBases = 0;
                    }

                    // output the mismatched base
                    md.push_back(rcRefBase);

                } else ++numMatchingBases;
            }

            // write the number of trailing matching bases
            if(numMatchingBases != 0) {
                md += boost::lexical_cast<std::string>(numMatchingBases);
            }
        
        } else {

            // handle unsupported CIGAR operation
            foundUnsupportedCigar = true;
            break;
        }
    }

    if(foundUnsupportedCigar) md = "UNSUPPORTED";
}

void
cigar_to_apath(const char* cigar,
               path_t& apath){

    using casava::blt_util::parse_unsigned;

    assert(NULL != cigar);

    apath.clear();

    path_segment lps;
    const char* cptr(cigar);
    while(*cptr) {
        path_segment ps;
        // expect sequences of digits and cigar codes:
        if((not isdigit(*cptr)) or (*cptr=='0')) unknown_cigar_error(cigar,cptr);
        ps.length = parse_unsigned(cptr);
        ps.type = cigar_code_to_segment_type(*cptr);
        if(ps.type == NONE) unknown_cigar_error(cigar,cptr);
        cptr++;
        if((ps.type == PAD) or (ps.length == 0)) continue;

        if(ps.type != lps.type){
            if(lps.type != NONE) apath.push_back(lps);
            lps = ps;
        } else {
            lps.length += ps.length;
        }
    }
    
    if(lps.type != NONE) apath.push_back(lps);
}

    

unsigned
apath_read_length(const path_t& apath) {

    unsigned val(0);
    const unsigned as(apath.size());
    for(unsigned i(0);i<as;++i){
        const path_segment& ps(apath[i]);
        if(not is_segment_type_read_length(ps.type)) continue;
        val += ps.length;
    }
    return val;
}



unsigned
apath_ref_length(const path_t& apath) {

    unsigned val(0);
    const unsigned as(apath.size());
    for(unsigned i(0);i<as;++i){
        const path_segment& ps(apath[i]);
        if(not is_segment_type_ref_length(ps.type)) continue;
        val += ps.length;
    }
    return val;
}



unsigned
apath_read_lead_size(const path_t& apath){
    unsigned val(0);
    const unsigned as(apath.size());
    for(unsigned i(0);i<as;++i){
        const path_segment& ps(apath[i]);
        if(not is_segment_type_unaligned_read_edge(ps.type)) return val;
        val += ps.length;
    }
    return val;
}



unsigned
apath_read_trail_size(const path_t& apath){
    unsigned val(0);
    const unsigned as(apath.size());
    for(unsigned i(0);i<as;++i){
        const path_segment& ps(apath[as-i-1]);
        if(not is_segment_type_unaligned_read_edge(ps.type)) return val;
        val += ps.length;
    }
    return val;
}



unsigned
apath_soft_clip_lead_size(const path_t& apath){
    unsigned val(0);
    const unsigned as(apath.size());
    for(unsigned i(0);i<as;++i){
        const path_segment& ps(apath[i]);
        if(SOFT_CLIP != ps.type) return val;
        val += ps.length;
    }
    return val;
}



unsigned
apath_soft_clip_trail_size(const path_t& apath){
    unsigned val(0);
    const unsigned as(apath.size());
    for(unsigned i(0);i<as;++i){
        const path_segment& ps(apath[as-i-1]);
        if(SOFT_CLIP != ps.type) return val;
        val += ps.length;
    }
    return val;
}



unsigned
apath_insert_lead_size(const path_t& apath){
    unsigned val(0);
    const unsigned as(apath.size());
    for(unsigned i(0);i<as;++i){
        const path_segment& ps(apath[i]);
        if(INSERT != ps.type) return val;
        val += ps.length;
    }
    return val;
}



unsigned
apath_insert_trail_size(const path_t& apath){
    unsigned val(0);
    const unsigned as(apath.size());
    for(unsigned i(0);i<as;++i){
        const path_segment& ps(apath[as-i-1]);
        if(INSERT != ps.type) return val;
        val += ps.length;
    }
    return val;
}



void
apath_soft_clip_clipper(path_t& apath,
                        unsigned& lead,
                        unsigned& trail){

    lead=0;
    trail=0;
    bool is_lead(true);
    path_t apath2;
    const unsigned as(apath.size());
    for(unsigned i(0);i<as;++i){
        const path_segment& ps(apath[i]);
        if(SOFT_CLIP == ps.type) {
            if(is_lead) { lead += ps.length; }
            else        { trail += ps.length; }
        } else {
            is_lead=false;
            trail=0;
            apath2.push_back(ps);
        }
    }
    apath=apath2;
}



void
apath_soft_clip_adder(path_t& apath,
                      const unsigned lead,
                      const unsigned trail){

    path_segment ps;
    ps.type = SOFT_CLIP;
    path_t apath2;
    if(lead>0) {
        ps.length = lead;
        apath2.push_back(ps);
    }
    apath2.insert(apath2.end(),apath.begin(),apath.end());
    if(trail>0) {
        ps.length = trail;
        apath2.push_back(ps);
    }
    apath=apath2;
}



std::pair<unsigned,unsigned>
get_nonclip_end_segments(const path_t& apath) {
    const unsigned as(apath.size());
    std::pair<unsigned,unsigned> res(as,as);
    bool is_nonclip(false);
    for(unsigned i(0);i<as;++i){
        const path_segment& ps(apath[i]);
        if(not (ps.type == SOFT_CLIP or
                ps.type == HARD_CLIP)) {
            if(not is_nonclip) res.first=i;
            is_nonclip=true;
            res.second=i;
        }
    }
    return res;
}



unsigned
apath_exon_count(const path_t& apath){
    unsigned val(1);
    const unsigned as(apath.size());
    for(unsigned i(0);i<as;++i){
        if(apath[i].type==SKIP) val++;
    }
    return val;
}



bool
is_soft_clipped(const path_t& apath){
    const unsigned as(apath.size());
    return ((as>0) and ((apath[0].type == SOFT_CLIP) or (apath[as-1].type==SOFT_CLIP)));
}



bool
is_edge_clipped(const path_t& apath){
    const unsigned as(apath.size());
    return ((as>0) and
            ((is_segment_type_unaligned_read_edge(apath[0].type)) or
             (is_segment_type_unaligned_read_edge(apath[as-1].type))));
}



bool
is_seq_swap(const path_t& apath) {
    bool is_last_indel(false);
    const unsigned as(apath.size());
    for(unsigned i(0);i<as;++i){
        const path_segment& ps(apath[as-i-1]);
        const bool is_indel((ps.type==INSERT) or (ps.type==DELETE));
        if(is_indel and is_last_indel) return true;
        is_last_indel=is_indel;
    }
    return false;
}



// take a first shot at the really simple stuff:
//
// 1) clipping only occurs on the edge
// 2) delete and skip cannot occur on edge
// 3) no unknown segments
// 4) no repeated segments
// 5) must contain at least one match segment
//
bool
is_apath_invalid(const path_t& apath,
                 const unsigned seq_length) {

    bool is_match(false);
    align_t last_type(NONE);
    const unsigned as(apath.size());
    for(unsigned i(0);i<as;++i){
        const path_segment& ps(apath[i]);
        
        if(ps.type==NONE) return true;
        if((i!=0) and ps.type==last_type) return true;

        if((i==0) or ((i+1)==as)){
            if((ps.type==DELETE) or (ps.type==SKIP)) return true;
        } else {
            if((ps.type==SOFT_CLIP) or (ps.type==HARD_CLIP)) return true;
        }

        if((not is_match) and (ps.type==MATCH)) is_match=true;
    }

    if(not is_match) return true;

    if(seq_length != apath_read_length(apath)) return true;

    return false;
}



bool
is_apath_starling_invalid(const path_t& apath) {

    const unsigned as(apath.size());
    for(unsigned i(0);i<as;++i){
        const path_segment& ps(apath[i]);
        
        if(ps.type==HARD_CLIP) return true;
        if(ps.type==PAD) return true;
    }
    return false;
}



}  // namespace ALIGNPATH