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

#include "blt_util/math_util.hh"

#include <cassert>
#include <cmath>
#include <cstdlib>

#include <iostream>
#include <iomanip>


#ifdef DEBUG_INDEL_CALL
#include "blt_util/log.hh"
#endif



/// find more accurate complement of posterior_probability:
static
double
prob_comp(const double* const dindel_prob,
          const unsigned cgt){

    double val(0.);
    for(unsigned gt(0);gt<STAR_DIINDEL::SIZE;++gt){
        if(gt == cgt) continue;
        val += dindel_prob[gt];
    }
    return val;
}



static
double
integrate_out_sites(const starling_deriv_options& client_dopt,
                    const double nonsite_path_lnp,
                    double p_on_site){

    return log_sum((p_on_site + client_dopt.site_lnprior),(nonsite_path_lnp + client_dopt.nonsite_lnprior));
}



///
///
void
starling_indel_call_pprob_digt(const starling_options& client_opt,
                               const starling_deriv_options& client_dopt,
                               const double indel_error_prob,
                               const double ref_error_prob,
                               const indel_key& ik,
                               const indel_data& id,
                               starling_diploid_indel& dindel){

    static const double loghalf(-std::log(2.));

    static const bool is_corrected_het_ratio(true);

    const bool is_breakpoint(ik.is_breakpoint());

    // get likelihood of each genotype:
    double lhood[STAR_DIINDEL::SIZE];
    for(unsigned gt(0);gt<STAR_DIINDEL::SIZE;++gt) lhood[gt] = 0.;

    const double indel_error_lnp(std::log(indel_error_prob));
    const double indel_real_lnp(std::log(1.-indel_error_prob));
    const double ref_error_lnp(std::log(ref_error_prob));
    const double ref_real_lnp(std::log(1.-ref_error_prob));

    typedef read_path_scores::alt_indel_t::const_iterator aiter;

    typedef indel_data::score_t::const_iterator siter;
    siter i(id.read_path_lnp.begin()), i_end(id.read_path_lnp.end());
    for(;i!=i_end;++i){
        const read_path_scores& path_lnp(i->second);

        // get alt path lnp:
        double alt_path_lnp(path_lnp.ref);
        if(is_use_alt_indel and (not path_lnp.alt_indel.empty()) ) {
            aiter j(path_lnp.alt_indel.begin()), j_end(path_lnp.alt_indel.end());
            for(;j!=j_end;++j){
                if(j->second>alt_path_lnp) alt_path_lnp=j->second;
            }
        }

        const double noindel_lnp(log_sum(alt_path_lnp+ref_real_lnp,path_lnp.indel+indel_error_lnp));
        const double hom_lnp(log_sum(alt_path_lnp+ref_error_lnp,path_lnp.indel+indel_real_lnp));

        double log_ref_prob(loghalf);
        double log_indel_prob(loghalf);
        if(is_corrected_het_ratio and (not is_breakpoint)) {
            // get expected het ratio
            const unsigned L(path_lnp.read_length);
            const unsigned B(client_opt.min_read_bp_flank);
            const unsigned spath_expect((2*B)>(L+1) ? 0 : (L+1)-(2*B) );
            const unsigned total_expect((2*spath_expect)+std::min(ik.length,spath_expect));
            if(total_expect>0){
                const double spath_prob(static_cast<double>(spath_expect)/static_cast<double>(total_expect));
                const double lpath_prob(1.-spath_prob);
                if(ik.type==INDEL::INSERT) {
                    log_ref_prob=std::log(spath_prob);
                    log_indel_prob=std::log(lpath_prob);
                } else if(ik.type==INDEL::DELETE){
                    log_ref_prob=std::log(lpath_prob);
                    log_indel_prob=std::log(spath_prob);
                } else {
                    assert(0); //unexpected indel type
                }
            }
        }
        const double het_lnp(log_sum(noindel_lnp+log_ref_prob,hom_lnp+log_indel_prob));

        lhood[STAR_DIINDEL::NOINDEL] += integrate_out_sites(client_dopt,path_lnp.nonsite,noindel_lnp);
        lhood[STAR_DIINDEL::HOM]     += integrate_out_sites(client_dopt,path_lnp.nonsite,hom_lnp);
        lhood[STAR_DIINDEL::HET]     += integrate_out_sites(client_dopt,path_lnp.nonsite,het_lnp);

#ifdef DEBUG_INDEL_CALL
        log_os << std::setprecision(8);
        log_os << "INDEL_CALL i,ref_lnp,indel_lnp,lhood(noindel),lhood(hom),lhood(het): " << i << " " << path_lnp.ref << " " << path_lnp.indel << " " << lhood[STAR_DIINDEL::NOINDEL] << " " << lhood[STAR_DIINDEL::HOM] << " " << lhood[STAR_DIINDEL::HET] << "\n";
#endif
    }

    // mult by prior distro to get unnormalized pprob:
    const double* lnprior(client_dopt.diploid_indel_state_lnprior);
    for(unsigned gt(0);gt<STAR_DIINDEL::SIZE;++gt){
        dindel.pprob[gt] = lhood[gt] + lnprior[gt];
    }

    // scale and exp pprob values:
    dindel.max_gt=0;
    dindel.max2_gt=1;
    double max(dindel.pprob[dindel.max_gt]);
    double max2(dindel.pprob[dindel.max2_gt]);
    for(int gt(1);gt<STAR_DIINDEL::SIZE;++gt){
        if(dindel.pprob[gt] > max){
            max2 = max;
            max = dindel.pprob[gt];
            dindel.max2_gt = dindel.max_gt;
            dindel.max_gt = gt;
        } else if(dindel.pprob[gt] > max2) {
            max2 = dindel.pprob[gt];
            dindel.max2_gt = gt;
        }
    }

    double sum(0.);
    for(unsigned gt(0);gt<STAR_DIINDEL::SIZE;++gt){
        dindel.pprob[gt] = std::exp(dindel.pprob[gt]-max);
        sum += dindel.pprob[gt];
    }

    // normalize:
    sum = 1./sum;
    for(unsigned gt(0);gt<STAR_DIINDEL::SIZE;++gt){
        dindel.pprob[gt] *= sum;
    }
#ifdef DEBUG_INDEL_CALL
        log_os << "INDEL_CALL pprob(noindel),pprob(hom),pprob(het): " << dindel.pprob[STAR_DIINDEL::NOINDEL] << " " << dindel.pprob[STAR_DIINDEL::HOM] << " " << dindel.pprob[STAR_DIINDEL::HET] << "\n";
#endif

    dindel.indel_qphred=error_prob_to_qphred(dindel.pprob[STAR_DIINDEL::NOINDEL]);
    dindel.max_gt_qphred=error_prob_to_qphred(prob_comp(dindel.pprob,dindel.max_gt));

    // old report criteria:
    //dindel.is_indel=(dindel.max_gt != STAR_DIINDEL::NOINDEL);
    dindel.is_indel=(dindel.indel_qphred != 0);
}



void
write_starling_diploid_indel(const starling_diploid_indel& dindel,
                             const starling_indel_report_info& iri,
                             std::ostream& os) {

    os << std::setprecision(10) << std::fixed;

    os << "type: " << iri.desc
       << " ref_upstream: " << iri.ref_upstream
       << " ref/indel: " << iri.ref_seq << "/" << iri.indel_seq
       << " ref_downstream: " << iri.ref_downstream
       << " Q(indel): " << dindel.indel_qphred
       << " max_gtype: " << STAR_DIINDEL::label(dindel.max_gt)
       << " Q(max_gtype): " << dindel.max_gt_qphred
       << " bp1_reads: " << iri.depth
       << " alt_reads: " << iri.n_q30_ref_reads+iri.n_q30_alt_reads
       << " indel_reads: " << iri.n_q30_indel_reads
       << " other_reads: " << iri.n_other_reads
       << " repeat_unit: " << iri.repeat_unit
       << " ref_repeat_count: " << iri.ref_repeat_count
       << " indel_repeat_count: " << iri.indel_repeat_count;

    os.unsetf(std::ios::fixed);
}



void
write_starling_diploid_indel_file(const starling_diploid_indel& dindel,
                                  const starling_indel_report_info& iri,
                                  std::ostream& os) {

    os << std::setprecision(10) << std::fixed;

    os << iri.desc
       << '\t' << iri.ref_upstream
       << '\t' << iri.ref_seq << "/" << iri.indel_seq
       << '\t' << iri.ref_downstream
       << '\t' << dindel.indel_qphred
       << '\t' << STAR_DIINDEL::label(dindel.max_gt)
       << '\t' << dindel.max_gt_qphred
       << '\t' << iri.depth
       << '\t' << iri.n_q30_ref_reads+iri.n_q30_alt_reads
       << '\t' << iri.n_q30_indel_reads
       << '\t' << iri.n_other_reads
       << '\t' << iri.repeat_unit
       << '\t' << iri.ref_repeat_count
       << '\t' << iri.indel_repeat_count;

    os.unsetf(std::ios::fixed);
}