// -*- 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_util.hh"
#include "starling_read_util.hh"

#include "starling/starling_shared.hh"

#include <cassert>

#include <sstream>



void
get_read_key_from_export_line(const export_line_parser& exl,
                              std::string& key){

    // note that this has been changed to the slower numerical version to 
    // guarantee that matching reads will be identified as such even after GROUPER
    // removes leading zeros
    //
    std::ostringstream oss;
    oss << exl.machine() << '_' << exl.run_number() << ':'
        << exl.lane() << ':' << exl.tile() << ':'
        << exl.x_coordinate() << ':' << exl.y_coordinate();
    key = oss.str();
}



pos_t
get_alignment_buffer_pos(const alignment& al){
    
    const pos_t lead(apath_read_lead_size(al.path));
    return al.pos-lead;
}



struct ddata {

    ddata(const unsigned read_size,
          const unsigned flank_size)
        : is_totaled(false), 
          fs(flank_size), fs2(fs*2),
          delta_size(std::max(1+fs2,read_size)-fs2) {
        for(unsigned i(0);i<delta_size;++i) { delta[i] = 0; }
    }

    void
    inc(const unsigned start_pos,
        const unsigned length) {
        assert(not is_totaled);
        delta[std::max(fs2,start_pos)-fs2] += 1;
        if((start_pos+length)<delta_size){ delta[start_pos+length] -= 1; }
    }

    int
    get(const unsigned pos) {
        if(not is_totaled) { total(); is_totaled=true; }
        return delta[std::min(delta_size-1,std::max(fs,pos)-fs)];
    }

private:
    void
    total() {
        for(unsigned i(1);i<delta_size;++i){ delta[i] += delta[i-1]; }
    }

    bool is_totaled;
    const unsigned fs;
    const unsigned fs2;
    const unsigned delta_size;
    int delta[MAX_READ_SIZE];
};




/// \brief generate mismatch-filter map in a way that scales
/// efficiently for large mismatch density windows:
///
//
// basic algo (described with 1-based indices):
//
// S=length(seq) F=length(flank)
//
// 1. delta list length is DL=max(1,S-2F)
// 2. for each position:
//	1. if position is mismatch, add 1 to delta list at position max(1,p-2F), and subtract 1 from delta list at position p+1 if p+1 <= DL
// 3. sum delta list to match-count list
// 4. match(pos) = delta[min(DL,max(1,pos-F))]
//
//..for window size 1:
//
// match:       M
// pos  : 1  2  3  4  5
// delta: x  1  0  0  x
//
void
create_mismatch_filter_map(const blt_options& client_opt,
                           const alignment& al,
                           const bam_seq_base& ref_seq,
                           const bam_seq_base& read_seq,
                           const unsigned read_begin,
                           const unsigned read_end,
                           bool* mismatch_filter_map,
                           int* mismatch_count_ns) {

    assert(read_begin<=read_end);
    const unsigned read_size(read_seq.size());
    assert(read_end<=read_size);
    const unsigned fs(client_opt.max_win_mismatch_flank_size);
    ddata dd(read_size,fs);

    bool is_mismatch[MAX_READ_SIZE];
    for(unsigned i(0);i<read_size;++i) is_mismatch[i]=false;

    pos_t ref_head_pos(al.pos);
    unsigned read_head_pos(0);

    using namespace ALIGNPATH;

    const unsigned as(al.path.size());
    for(unsigned i(0);i<as;++i){
        const path_segment& ps(al.path[i]);

        // only fair treatment for soft-clipped sequence is
        // to treat it as an insert (ie. as a mismatch):
        //
        if       ((ps.type == INSERT) or (ps.type == SOFT_CLIP)) {
            if((i!=0) and ((i+1)!=as)) dd.inc(read_head_pos,ps.length);
            read_head_pos += ps.length;

        } else if(ps.type == DELETE) {
            dd.inc(read_head_pos,0);
            ref_head_pos += ps.length;

        } else if(ps.type == MATCH) {
            for(unsigned j(0);j<ps.length;++j){
                const unsigned read_pos(read_head_pos+j);
                if((read_pos < read_begin) or (read_pos >= read_end)) continue; // allow for read end trimming
                const pos_t ref_pos(ref_head_pos+static_cast<pos_t>(j));
                
                if(read_seq.get_char(read_pos) != ref_seq.get_char(ref_pos)){
                    is_mismatch[read_pos]=true;
                    dd.inc(read_pos,1);
                } 
            }   
            read_head_pos += ps.length;
            ref_head_pos += ps.length;

        } else {
            assert(0);
        }
    }

    // set mismatch filter value:
    const int max_pass(static_cast<int>(client_opt.max_win_mismatch));
    for(unsigned i(0);i<read_size;++i){
        const int del(dd.get(i));
        mismatch_count_ns[i] = del - is_mismatch[i];
        mismatch_filter_map[i] = (max_pass < del);
    }
}