// -*- 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 "blt_common/adjust_joint_eprob.hh"

#include <cmath>

#include <algorithm>
#include <iostream>
#include <vector>


typedef unsigned icall_t; 
typedef std::vector<icall_t> icalls_t;



struct sort_icall_by_eprob {
    sort_icall_by_eprob(const snp_pos_info& pi) : _pi(pi) {}

    bool 
    operator()(const icall_t& a, const icall_t& b) const {
        return (_pi.calls[a].qscore > _pi.calls[b].qscore);
    }

    const snp_pos_info& _pi;
};



//#define DEBUG_ADJUST

static
void
adjust_icalls_eprob(const blt_options& opt,
                    icalls_t& ic,
                    const snp_pos_info& pi,
                    std::vector<float>& dependent_eprob) {

    const unsigned ic_size(ic.size());

#ifdef DEBUG_ADJUST 
    for(unsigned i(0);i<ic_size;++i){
        base_call& bi(pi.calls[ic[i]]);
        std::cerr << "BEFORE: " << i << " " << bi.is_neighbor_mismatch << " " << dependent_eprob[ic[i]] << "\n";
    }
#endif

    static const bool is_use_vexp_frac(true);

    // produce weighted fraction of reads with a neighboring mismatch:
    double vexp_frac(opt.bsnp_ssd_no_mismatch);
    if(is_use_vexp_frac){
        static const double lnran(std::log(0.75));
        double num(0);
        double den(0);
        for(unsigned i(0);i<ic_size;++i){
            const base_call& bi(pi.calls[ic[i]]);
            const double eprob(dependent_eprob[ic[i]]);
            const double weight(lnran-std::log(eprob));
            den += weight;
            if(bi.is_neighbor_mismatch) { num += weight; }
        }
        double mismatch_frac(0);
        if(ic_size and (den>0.)) mismatch_frac=(num/den);
        vexp_frac=(1.-mismatch_frac)*opt.bsnp_ssd_no_mismatch+mismatch_frac*opt.bsnp_ssd_one_mismatch;
    }

    static const double dep_converge_prob(0.75);
    const bool is_limit_vexp_iterations(opt.max_vexp_iterations>0);
    const int max_vexp_iterations(opt.max_vexp_iterations);

    const bool is_limit_vexp(opt.is_min_vexp);
    const double min_vexp(opt.min_vexp);

    std::sort(ic.begin(),ic.end(),sort_icall_by_eprob(pi));
    double vexp(1.);
    for(unsigned i(0);i<ic_size;++i){
        const base_call& bi(pi.calls[ic[i]]);
        const double eprob(dependent_eprob[ic[i]]);
        const double val(std::pow(eprob,vexp));
        const double frac((1.-val)/(1.-eprob));
        dependent_eprob[ic[i]]=std::max(eprob,frac*val+(1-frac)*dep_converge_prob);

        if(is_limit_vexp_iterations and (static_cast<int>(i)>=max_vexp_iterations)) continue;

        double next_vexp(vexp);
        if(is_use_vexp_frac) {
            next_vexp *= (1.-vexp_frac);
        } else {
            if(bi.is_neighbor_mismatch) {
                next_vexp *= (1.-opt.bsnp_ssd_one_mismatch);
            } else {
                next_vexp *= (1.-opt.bsnp_ssd_no_mismatch);
            }
        }
        if(is_limit_vexp) {
            vexp = std::max(min_vexp,next_vexp);
        } else {
            vexp = next_vexp;
        }
    }

#ifdef DEBUG_ADJUST
    for(unsigned i(0);i<ic_size;++i){
        base_call& bi(pi.calls[ic[i]]);
        std::cerr << "AFTER: " << i << " " << bi.is_neighbor_mismatch << " " << dependent_eprob[ic[i]] << "\n";
    }
#endif
}



// downweight qvals at a site based on dependency effects:
//
void
adjust_joint_eprob(const blt_options& opt,
                   const snp_pos_info& pi,
                   std::vector<float>& dependent_eprob,
                   const bool is_dependent) {

    const unsigned n_calls(pi.calls.size());
    for(unsigned i(0);i<n_calls;++i){
        dependent_eprob.push_back(pi.calls[i].error_prob());
    }
    
    if(not is_dependent) return;

    // split icalls into fwd and reverse strand and allele types:
    //
    static const unsigned group_size(8); // (is_fwd*base_id)
    icalls_t icalls[group_size];
    for(unsigned i(0);i<n_calls;++i){
        // exclude q2's and filtered bases:
        const base_call& b(pi.calls[i]);
#ifdef NOREFFILTER
        if(b.is_call_filter and (b.base_id != pi.ref_base_id)) continue;
#else
        if(b.is_call_filter) continue;
#endif

        //if(b.is_call_filter or (b.error_prob>=0.5)) continue;
        if(b.qscore<3) continue;

        const unsigned group_index((b.is_fwd_strand)+(2*b.base_id));
        icalls[group_index].push_back(i);
    }

    // process each isize array:
    for(unsigned i(0);i<group_size;++i){
        adjust_icalls_eprob(opt,icalls[i],pi,dependent_eprob);
    } 
}