agv_navigation_ws/fast_gicp/thirdparty/nvbio/nvbio-aln-diff/pe_analyzer.cpp

/*
 * nvbio
 * Copyright (c) 2011-2014, NVIDIA CORPORATION. All rights reserved.
 * 
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *    * Redistributions of source code must retain the above copyright
 *      notice, this list of conditions and the following disclaimer.
 *    * Redistributions in binary form must reproduce the above copyright
 *      notice, this list of conditions and the following disclaimer in the
 *      documentation and/or other materials provided with the distribution.
 *    * Neither the name of the NVIDIA CORPORATION nor the
 *      names of its contributors may be used to endorse or promote products
 *      derived from this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL NVIDIA CORPORATION BE LIABLE FOR ANY
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

#include <nvbio-aln-diff/pe_analyzer.h>
#include <nvbio-aln-diff/html.h>
#include <nvbio-aln-diff/utils.h>

namespace nvbio {
namespace alndiff {

PEAnalyzer::PEAnalyzer(Filter& filter, const bool id_check) :
    m_filter( filter ),
    m_id_check( id_check )
{
    n = 0;
    n_mismatched = 0;
}

void PEAnalyzer::push(
    const AlignmentPair& alnL,
    const AlignmentPair& alnR)
{
    if ((m_id_check && alnL[0].read_id  != alnR[0].read_id)  ||
        (m_id_check && alnL[1].read_id  != alnR[1].read_id)  ||
        alnL[0].read_len != alnR[0].read_len ||
        alnL[1].read_len != alnR[1].read_len)
    {
        n_mismatched++;
        return;
    }

    mapped.push( alnL.is_mapped(), alnR.is_mapped() );
    paired.push( alnL.is_mapped_paired(), alnR.is_mapped_paired() );
    unique.push( alnL.is_unique_paired(), alnR.is_unique_paired() );
    ambiguous.push( alnL.is_ambiguous_paired(), alnR.is_ambiguous_paired() );

    if ((alnL.is_mapped_paired()    == true) && (alnR.is_mapped_paired()    == false)) paired_L_not_R_by_mapQ.push( log_bin( alnL.mapQ() ) );
    if ((alnR.is_mapped_paired()    == true) && (alnL.is_mapped_paired()    == false)) paired_R_not_L_by_mapQ.push( log_bin( alnR.mapQ() ) );
    if ((alnL.is_unique_paired()    == true) && (alnR.is_unique_paired()    == false)) unique_L_not_R_by_mapQ.push( log_bin( alnL.mapQ() ) );
    if ((alnR.is_unique_paired()    == true) && (alnL.is_unique_paired()    == false)) unique_R_not_L_by_mapQ.push( log_bin( alnR.mapQ() ) );
    if ((alnL.is_ambiguous_paired() == true) && (alnR.is_ambiguous_paired() == false)) ambiguous_L_not_R_by_mapQ.push( log_bin( alnL.mapQ() ) );
    if ((alnR.is_ambiguous_paired() == true) && (alnL.is_ambiguous_paired() == false)) ambiguous_R_not_L_by_mapQ.push( log_bin( alnR.mapQ() ) );

    // edit-distance correlation stats
    //if (alnL.is_mapped_paired())    sec_ed_by_ed_l.push( alnL.ed()+1, alnL.has_second() ? alnL.sec_ed()+1 : 0u );
    //else                            sec_ed_by_ed_l.push( 0u, 0u );

    // edit-distance correlation stats
    //if (alnR.is_mapped_paired())    sec_ed_by_ed_r.push( alnR.ed()+1, alnR.has_second() ? alnR.sec_ed()+1 : 0u );
    //else                            sec_ed_by_ed_r.push( 0u, 0u );

    if (alnL.is_mapped_paired() && alnR.is_mapped_paired())
    {
        const uint32 mapQ_bin = log_bin( alnR.mapQ() );

        uint32 read_flags = 0u;

        if ((alnL[0].ref_id != alnR[0].ref_id) &&
            (alnL[1].ref_id != alnR[1].ref_id))
        {
            n_different_ref12.push( mapQ_bin );
            n_different_ref.push( mapQ_bin );
            n_distant12.push( mapQ_bin );
            n_distant.push( mapQ_bin );

            if ((alnL.is_unique_paired() == true) &&
                (alnR.is_unique_paired() == true))
            {
                n_different_ref_unique.push( mapQ_bin );
                n_distant_unique.push( mapQ_bin );
            }
            if ((alnL.is_ambiguous_paired() == false) &&
                (alnR.is_ambiguous_paired() == false))
            {
                n_different_ref_not_ambiguous.push( mapQ_bin );
                n_distant_not_ambiguous.push( mapQ_bin );
            }

            read_flags |= Filter::DISTANT;
            read_flags |= Filter::DIFFERENT_REF;
        }
        else if (alnL[0].ref_id != alnR[0].ref_id)
        {
            n_different_ref1.push( mapQ_bin );
            n_different_ref.push( mapQ_bin );
            n_distant1.push( mapQ_bin );
            n_distant.push( mapQ_bin );

            if ((alnL.is_unique_paired() == true) &&
                (alnR.is_unique_paired() == true))
            {
                n_different_ref_unique.push( mapQ_bin );
                n_distant_unique.push( mapQ_bin );
            }
            if ((alnL.is_ambiguous_paired() == false) &&
                (alnR.is_ambiguous_paired() == false))
            {
                n_different_ref_not_ambiguous.push( mapQ_bin );
                n_distant_not_ambiguous.push( mapQ_bin );
            }

            read_flags |= Filter::DISTANT;
            read_flags |= Filter::DIFFERENT_REF;
        }
        else if (alnL[1].ref_id != alnR[1].ref_id)
        {
            n_different_ref2.push( mapQ_bin );
            n_different_ref.push( mapQ_bin );
            n_distant2.push( mapQ_bin );
            n_distant.push( mapQ_bin );

            if ((alnL.is_unique_paired() == true) &&
                (alnR.is_unique_paired() == true))
            {
                n_different_ref_unique.push( mapQ_bin );
                n_distant_unique.push( mapQ_bin );
            }
            if ((alnL.is_ambiguous_paired() == false) &&
                (alnR.is_ambiguous_paired() == false))
            {
                n_different_ref_not_ambiguous.push( mapQ_bin );
                n_distant_not_ambiguous.push( mapQ_bin );
            }

            read_flags |= Filter::DISTANT;
            read_flags |= Filter::DIFFERENT_REF;
        }
        else if (alndiff::distant( alnL[0], alnR[0] ) && alndiff::distant( alnL[1], alnR[1] ))
        {
            n_distant12.push( mapQ_bin );
            n_distant.push( mapQ_bin );

            if ((alnL.is_unique_paired() == true) &&
                (alnR.is_unique_paired() == true))
                n_distant_unique.push( mapQ_bin );
            if ((alnL.is_ambiguous_paired() == false) &&
                (alnR.is_ambiguous_paired() == false))
                n_distant_not_ambiguous.push( mapQ_bin );

            read_flags |= Filter::DISTANT;
        }
        else if (alndiff::distant( alnL[0], alnR[0] ))
        {
            n_distant1.push( mapQ_bin );
            n_distant.push( mapQ_bin );

            if ((alnL.is_unique_paired() == true) &&
                (alnR.is_unique_paired() == true))
                n_distant_unique.push( mapQ_bin );
            if ((alnL.is_ambiguous_paired() == false) &&
                (alnR.is_ambiguous_paired() == false))
                n_distant_not_ambiguous.push( mapQ_bin );

            read_flags |= Filter::DISTANT;
        }
        else if (alndiff::distant( alnL[1], alnR[1] ))
        {
            n_distant2.push( mapQ_bin );
            n_distant.push( mapQ_bin );

            if ((alnL.is_unique_paired() == true) &&
                (alnR.is_unique_paired() == true))
                n_distant_unique.push( mapQ_bin );
            if ((alnL.is_ambiguous_paired() == false) &&
                (alnR.is_ambiguous_paired() == false))
                n_distant_not_ambiguous.push( mapQ_bin );

            read_flags |= Filter::DISTANT;
        }

        if ((alnL[0].is_rc() != alnR[0].is_rc()) &&
            (alnL[1].is_rc() != alnR[1].is_rc()))
        {
            n_discordant12.push( mapQ_bin );
            n_discordant.push( mapQ_bin );

            if ((alnL.is_unique_paired() == true) &&
                (alnR.is_unique_paired() == true))
                n_discordant_unique.push( mapQ_bin );
            if ((alnL.is_ambiguous_paired() == false) &&
                (alnR.is_ambiguous_paired() == false))
                n_discordant_not_ambiguous.push( mapQ_bin );

            read_flags |= Filter::DISCORDANT;
        }
        if (alnL[0].is_rc() != alnR[0].is_rc())
        {
            n_discordant1.push( mapQ_bin );
            n_discordant.push( mapQ_bin );

            if ((alnL.is_unique_paired() == true) &&
                (alnR.is_unique_paired() == true))
                n_discordant_unique.push( mapQ_bin );
            if ((alnL.is_ambiguous_paired() == false) &&
                (alnR.is_ambiguous_paired() == false))
                n_discordant_not_ambiguous.push( mapQ_bin );

            read_flags |= Filter::DISCORDANT;
        }
        else if (alnL[1].is_rc() != alnR[1].is_rc())
        {
            n_discordant2.push( mapQ_bin );
            n_discordant.push( mapQ_bin );

            if ((alnL.is_unique_paired() == true) &&
                (alnR.is_unique_paired() == true))
                n_discordant_unique.push( mapQ_bin );
            if ((alnL.is_ambiguous_paired() == false) &&
                (alnR.is_ambiguous_paired() == false))
                n_discordant_not_ambiguous.push( mapQ_bin );

            read_flags |= Filter::DISCORDANT;
        }

        const uint32 length_bin = read_length_bin( alnL.read_len() );

        // generic stats
        {
            // score
            m_filter( al_stats.higher_score.push( alnL.score(), alnR.score(), length_bin, alnL.mapQ(), alnR.mapQ() ), read_flags, Filter::SCORE, alnL.read_id() );

            // ed
            m_filter( al_stats.lower_ed.push( alnL.ed(), alnR.ed(), length_bin, alnL.mapQ(), alnR.mapQ() ), read_flags, Filter::ED, alnL.read_id() );

            // mapQ
            m_filter( al_stats.higher_mapQ.push( alnL.mapQ(), alnR.mapQ(), length_bin, alnL.mapQ(), alnR.mapQ() ), read_flags, Filter::MAPQ, alnL.read_id() );

            // longer mapping
            al_stats.longer_mapping.push( alnL.mapped_read_bases(), alnR.mapped_read_bases(), length_bin, alnL.mapQ(), alnR.mapQ() );

            al_stats.lower_subs.push( alnL.subs(), alnR.subs(), length_bin, alnL.mapQ(), alnR.mapQ() );
            m_filter( al_stats.lower_mms.push(  alnL.n_mm(), alnR.n_mm(), length_bin, alnL.mapQ(), alnR.mapQ() ), read_flags, Filter::MMS,  alnL.read_id() );
            m_filter( al_stats.lower_ins.push(  alnL.ins(),  alnR.ins(),  length_bin, alnL.mapQ(), alnR.mapQ() ), read_flags, Filter::INS,  alnL.read_id() );
            m_filter( al_stats.lower_dels.push( alnL.dels(), alnR.dels(), length_bin, alnL.mapQ(), alnR.mapQ() ), read_flags, Filter::DELS, alnL.read_id() );

            al_stats.higher_pos.push( alnL[0].pos, alnR[0].pos, length_bin, alnL.mapQ(), alnR.mapQ() );
              // TODO: this is a 2d quantity

            //if (alnL.has_second())
                sec_score_by_score_l.push( log_bin( alnL.score() ), log_bin( alnL.score() - alnL.sec_score() ) );

            //if (alnR.has_second())
                sec_score_by_score_r.push( log_bin( alnR.score() ), log_bin( alnR.score() - alnR.sec_score() ) );
        }
        if (read_flags & Filter::DISTANT)
        {
            // ed
            distant_stats.lower_ed.push( alnL.ed(), alnR.ed(), length_bin, alnL.mapQ(), alnR.mapQ() );

            // mapQ
            distant_stats.higher_mapQ.push( alnL.mapQ(), alnR.mapQ(), length_bin, alnL.mapQ(), alnR.mapQ() );

            // longer mapping
            distant_stats.longer_mapping.push( alnL.mapped_read_bases(), alnR.mapped_read_bases(), length_bin, alnL.mapQ(), alnR.mapQ() );

            distant_stats.lower_subs.push( alnL.subs(), alnR.subs(), length_bin, alnL.mapQ(), alnR.mapQ() );
            distant_stats.lower_mms.push(  alnL.n_mm(), alnR.n_mm(), length_bin, alnL.mapQ(), alnR.mapQ() );
            distant_stats.lower_ins.push(  alnL.ins(),  alnR.ins(),  length_bin, alnL.mapQ(), alnR.mapQ() );
            distant_stats.lower_dels.push( alnL.dels(), alnR.dels(), length_bin, alnL.mapQ(), alnR.mapQ() );

            distant_stats.higher_pos.push( alnL[0].pos,  alnR[0].pos,  length_bin, alnL.mapQ(), alnR.mapQ() );
              // TODO: this is a 2d quantity
        }
        if (read_flags & Filter::DISCORDANT)
        {
            // ed
            discordant_stats.lower_ed.push( alnL.ed(), alnR.ed(), length_bin, alnL.mapQ(), alnR.mapQ() );

            // mapQ
            discordant_stats.higher_mapQ.push( alnL.mapQ(), alnR.mapQ(), length_bin, alnL.mapQ(), alnR.mapQ() );

            // longer mapping
            discordant_stats.longer_mapping.push( alnL.mapped_read_bases(), alnR.mapped_read_bases(), length_bin, alnL.mapQ(), alnR.mapQ() );

            discordant_stats.lower_subs.push( alnL.subs(), alnR.subs(), length_bin, alnL.mapQ(), alnR.mapQ() );
            discordant_stats.lower_mms.push(  alnL.n_mm(), alnR.n_mm(), length_bin, alnL.mapQ(), alnR.mapQ() );
            discordant_stats.lower_ins.push(  alnL.ins(),  alnR.ins(),  length_bin, alnL.mapQ(), alnR.mapQ() );
            discordant_stats.lower_dels.push( alnL.dels(), alnR.dels(), length_bin, alnL.mapQ(), alnR.mapQ() );

            discordant_stats.higher_pos.push( alnL[0].pos,  alnR[0].pos,  length_bin, alnL.mapQ(), alnR.mapQ() );
              // TODO: this is a 2d quantity
        }
    }

    ++n;
}

namespace {

void generate_summary_header(FILE* html_output)
{
    html::tr_object tr( html_output, NULL );
    html::th_object( html_output, html::FORMATTED, NULL, "" );
    html::th_object( html_output, html::FORMATTED, NULL, "better" );
    html::th_object( html_output, html::FORMATTED, "class", "red", NULL, "worse" );
}
template <typename StatsType>
void generate_summary_cell(FILE* html_output, const std::string file_name, const char* type, const StatsType& stats, const uint32 n, const char* cls)
{
    char link_name[1024];
    sprintf( link_name, "<a href=\"%s\">%s</a>", local_file( file_name ), type );
    html::tr_object tr( html_output, "class", cls, NULL );
    html::th_object( html_output, html::FORMATTED, NULL, link_name );
    html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%", 100.0f * float(stats.l.diff_hist.all_but(0))/float(n) );
    html::td_object( html_output, html::FORMATTED, "class", "pink", NULL, "%5.2f %%", 100.0f * float(stats.r.diff_hist.all_but(0))/float(n) );
}

} // anonymous namespace

void PEAnalyzer::generate_report(const char* aln_file_nameL, const char* aln_file_nameR, const char* report)
{
    if (report == NULL)
        return;

    const std::string mapped_bps_name     = generate_file_name( report, "mapped-bps" );
    const std::string score_name          = generate_file_name( report, "score" );
    const std::string ed_name             = generate_file_name( report, "ed" );
    const std::string mapQ_name           = generate_file_name( report, "mapQ" );
    const std::string mms_name            = generate_file_name( report, "mms" );
    const std::string ins_name            = generate_file_name( report, "ins" );
    const std::string dels_name           = generate_file_name( report, "dels" );
    const std::string pos_name            = generate_file_name( report, "pos" );

    const std::string distant_mapped_bps_name       = generate_file_name( report, "distant_stats.mapped-bps" );
    const std::string distant_ed_name               = generate_file_name( report, "distant_stats.ed" );
    const std::string distant_mapQ_name             = generate_file_name( report, "distant_stats.mapQ" );
    const std::string distant_mms_name              = generate_file_name( report, "distant_stats.mms" );
    const std::string distant_ins_name              = generate_file_name( report, "distant_stats.ins" );
    const std::string distant_dels_name             = generate_file_name( report, "distant_stats.dels" );
    const std::string distant_pos_name              = generate_file_name( report, "distant_stats.pos" );

    const std::string discordant_mapped_bps_name    = generate_file_name( report, "discordant_stats.mapped-bps" );
    const std::string discordant_ed_name            = generate_file_name( report, "discordant_stats.ed" );
    const std::string discordant_mapQ_name          = generate_file_name( report, "discordant_stats.mapQ" );
    const std::string discordant_mms_name           = generate_file_name( report, "discordant_stats.mms" );
    const std::string discordant_ins_name           = generate_file_name( report, "discordant_stats.ins" );
    const std::string discordant_dels_name          = generate_file_name( report, "discordant_stats.dels" );
    const std::string discordant_pos_name           = generate_file_name( report, "discordant_stats.pos" );

    generate_table( mapped_bps_name.c_str(),    "paired L & R", "mapped bps",    "mapped bps diff",     al_stats.longer_mapping.bin_type(), al_stats.longer_mapping.l, al_stats.longer_mapping.r,   n, paired.L_and_R );
    generate_table( score_name.c_str(),         "paired L & R", "score",         "score diff",          al_stats.higher_score.bin_type(),   al_stats.higher_score.l,   al_stats.higher_score.r,     n, paired.L_and_R );
    generate_table( ed_name.c_str(),            "paired L & R", "edit distance", "edit distance diff",  al_stats.lower_ed.bin_type(),       al_stats.lower_ed.l,       al_stats.lower_ed.r,         n, paired.L_and_R );
    generate_table( mapQ_name.c_str(),          "paired L & R", "mapQ",          "mapQ diff",           al_stats.higher_mapQ.bin_type(),    al_stats.higher_mapQ.l,    al_stats.higher_mapQ.r,      n, paired.L_and_R );
    generate_table( mms_name.c_str(),           "paired L & R", "mms",           "mms diff",            al_stats.lower_mms.bin_type(),      al_stats.lower_mms.l,      al_stats.lower_mms.r,        n, paired.L_and_R );
    generate_table( ins_name.c_str(),           "paired L & R", "ins",           "ins diff",            al_stats.lower_ins.bin_type(),      al_stats.lower_ins.l,      al_stats.lower_ins.r,        n, paired.L_and_R );
    generate_table( dels_name.c_str(),          "paired L & R", "dels",          "dels diff",           al_stats.lower_dels.bin_type(),     al_stats.lower_dels.l,     al_stats.lower_dels.r,       n, paired.L_and_R );
    generate_table( pos_name.c_str(),           "paired L & R", "position",      "distance",            al_stats.higher_pos.bin_type(),     al_stats.higher_pos.l,     al_stats.higher_pos.r,       n, paired.L_and_R, false );

    generate_table( distant_mapped_bps_name.c_str(),    "distant", "mapped bps",    "mapped bps diff",      distant_stats.longer_mapping.bin_type(), distant_stats.longer_mapping.l, distant_stats.longer_mapping.r,   n, n_distant.count );
    generate_table( distant_ed_name.c_str(),            "distant", "edit distance", "edit distance diff",   distant_stats.lower_ed.bin_type(),       distant_stats.lower_ed.l,       distant_stats.lower_ed.r,         n, n_distant.count );
    generate_table( distant_mapQ_name.c_str(),          "distant", "mapQ",          "mapQ diff",            distant_stats.higher_mapQ.bin_type(),    distant_stats.higher_mapQ.l,    distant_stats.higher_mapQ.r,      n, n_distant.count );
    generate_table( distant_mms_name.c_str(),           "distant", "mms",           "mms diff",             distant_stats.lower_mms.bin_type(),      distant_stats.lower_mms.l,      distant_stats.lower_mms.r,        n, n_distant.count );
    generate_table( distant_ins_name.c_str(),           "distant", "ins",           "ins diff",             distant_stats.lower_ins.bin_type(),      distant_stats.lower_ins.l,      distant_stats.lower_ins.r,        n, n_distant.count );
    generate_table( distant_dels_name.c_str(),          "distant", "dels",          "dels diff",            distant_stats.lower_dels.bin_type(),     distant_stats.lower_dels.l,     distant_stats.lower_dels.r,       n, n_distant.count );
    generate_table( distant_pos_name.c_str(),           "distant", "position",      "distance",             distant_stats.higher_pos.bin_type(),     distant_stats.higher_pos.l,     distant_stats.higher_pos.r,       n, n_distant.count, false );

    generate_table( discordant_mapped_bps_name.c_str(),    "discordant", "mapped bps",   "mapped bps diff",     discordant_stats.longer_mapping.bin_type(), discordant_stats.longer_mapping.l, discordant_stats.longer_mapping.r,   n, n_discordant.count );
    generate_table( discordant_ed_name.c_str(),            "discordant", "edit distance","edit distance diff",  discordant_stats.lower_ed.bin_type(),       discordant_stats.lower_ed.l,       discordant_stats.lower_ed.r,         n, n_discordant.count );
    generate_table( discordant_mapQ_name.c_str(),          "discordant", "mapQ",         "mapQ diff",           discordant_stats.higher_mapQ.bin_type(),    discordant_stats.higher_mapQ.l,    discordant_stats.higher_mapQ.r,      n, n_discordant.count );
    generate_table( discordant_mms_name.c_str(),           "discordant", "mms",          "mms diff",            discordant_stats.lower_mms.bin_type(),      discordant_stats.lower_mms.l,      discordant_stats.lower_mms.r,        n, n_discordant.count );
    generate_table( discordant_ins_name.c_str(),           "discordant", "ins",          "ins diff",            discordant_stats.lower_ins.bin_type(),      discordant_stats.lower_ins.l,      discordant_stats.lower_ins.r,        n, n_discordant.count );
    generate_table( discordant_dels_name.c_str(),          "discordant", "dels",         "dels diff",           discordant_stats.lower_dels.bin_type(),     discordant_stats.lower_dels.l,     discordant_stats.lower_dels.r,       n, n_discordant.count );
    generate_table( discordant_pos_name.c_str(),           "discordant", "position",     "distance",            discordant_stats.higher_pos.bin_type(),     discordant_stats.higher_pos.l,     discordant_stats.higher_pos.r,       n, n_discordant.count, false );

    FILE* html_output = fopen( report, "w" );
    if (html_output == NULL)
    {
        log_warning( stderr, "unable to write HTML report \"%s\"\n", report );
        return;
    }

    std::vector<std::string> log_bins;
    for (uint32 y = 0; y < 11; ++y)
    {
        if (y == 0)
            log_bins.push_back( "0" );
        else if (y == 1)
            log_bins.push_back( "1" );
        else if (y == 2)
            log_bins.push_back( "2" );
        else
        {
            char buffer[16];
            sprintf(buffer, "2^%u", y-1);
            log_bins.push_back( buffer );
        }
    }
    log_bins.push_back( "inf" );

    std::vector<std::string> ed_bins;
    ed_bins.push_back("-");
    for (uint32 y = 0; y < 16; ++y)
    {
        char buffer[16];
        sprintf(buffer, "%u", y);
        ed_bins.push_back( buffer );
    }

    const Histogram<8> cum_different_ref   = reverse_cumulative( n_different_ref );
    const Histogram<8> cum_different_ref12 = reverse_cumulative( n_different_ref12 );
    const Histogram<8> cum_different_ref1  = reverse_cumulative( n_different_ref1 );
    const Histogram<8> cum_different_ref2  = reverse_cumulative( n_different_ref2 );
    const Histogram<8> cum_different_ref_unique         = reverse_cumulative( n_different_ref_unique );
    const Histogram<8> cum_different_ref_not_ambiguous  = reverse_cumulative( n_different_ref_not_ambiguous );

    const Histogram<8> cum_distant   = reverse_cumulative( n_distant );
    const Histogram<8> cum_distant12 = reverse_cumulative( n_distant12 );
    const Histogram<8> cum_distant1  = reverse_cumulative( n_distant1 );
    const Histogram<8> cum_distant2  = reverse_cumulative( n_distant2 );
    const Histogram<8> cum_distant_unique          = reverse_cumulative( n_distant_unique );
    const Histogram<8> cum_distant_not_ambiguous   = reverse_cumulative( n_distant_not_ambiguous );

    const Histogram<8> cum_discordant   = reverse_cumulative( n_discordant );
    const Histogram<8> cum_discordant12 = reverse_cumulative( n_discordant12 );
    const Histogram<8> cum_discordant1  = reverse_cumulative( n_discordant1 );
    const Histogram<8> cum_discordant2  = reverse_cumulative( n_discordant2 );
    const Histogram<8> cum_discordant_unique        = reverse_cumulative( n_discordant_unique );
    const Histogram<8> cum_discordant_not_ambiguous = reverse_cumulative( n_discordant_not_ambiguous );

    const uint32 HI_MAPQ_BIN = 6; // >= 32
    {
        html::html_object html( html_output );
        {
            const char* meta_list = "<meta http-equiv=\"refresh\" content=\"5\" />";

            html::header_object hd( html_output, "nv-aln-diff report", html::style(), meta_list );
            {
                html::body_object body( html_output );

                //
                // alignment stats
                //
                {
                    html::table_object table( html_output, "alignment-stats", "stats", "alignment stats" );
                    {
                        html::tr_object tr( html_output, NULL );
                        html::th_object( html_output, html::FORMATTED, NULL, "" );
                        html::th_object( html_output, html::FORMATTED, NULL, "L = %s", aln_file_nameL );
                        html::th_object( html_output, html::FORMATTED, NULL, "R = %s", aln_file_nameR );
                        html::th_object( html_output, html::FORMATTED, NULL, "L & R" );
                        html::th_object( html_output, html::FORMATTED, NULL, "L \\ R" );
                        html::th_object( html_output, html::FORMATTED, NULL, "R \\ L" );
                    }
                    {
                        html::tr_object tr( html_output, "class", "alt", NULL );
                        html::th_object( html_output, html::FORMATTED, NULL, "mapped");
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%", 100.0f * float(mapped.L) / float(n) );
                        html::td_object( html_output, html::FORMATTED, "class", "pink",  NULL, "%5.2f %%", 100.0f * float(mapped.R) / float(n) );
                        html::td_object( html_output, html::FORMATTED, "class", "green", NULL, "%5.2f %%", 100.0f * float(mapped.L_and_R) / float(n) );
                        html::td_object( html_output, html::FORMATTED,                   NULL, "%5.2f %%", 100.0f * float(mapped.L_not_R) / float(n) );
                        html::td_object( html_output, html::FORMATTED, "class", "pink",  NULL, "%5.2f %%", 100.0f * float(mapped.R_not_L) / float(n) );
                    }
                    {
                        html::tr_object tr( html_output, NULL );
                        html::th_object( html_output, html::FORMATTED, NULL, "paired");
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%", 100.0f * float(paired.L) / float(n) );
                        html::td_object( html_output, html::FORMATTED, "class", "pink",  NULL, "%5.2f %%", 100.0f * float(paired.R) / float(n) );
                        html::td_object( html_output, html::FORMATTED, "class", "green", NULL, "%5.2f %%", 100.0f * float(paired.L_and_R) / float(n) );
                        html::td_object( html_output, html::FORMATTED,                   NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(paired.L_not_R) / float(n), 100.0f * float(paired_L_not_R_by_mapQ[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, "class", "pink",  NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(paired.R_not_L) / float(n), 100.0f * float(paired_R_not_L_by_mapQ[HI_MAPQ_BIN]) / float(n) );
                    }
                    {
                        html::tr_object tr( html_output, "class", "alt", NULL );
                        html::th_object( html_output, html::FORMATTED, NULL, "unique");
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%", 100.0f * float(unique.L) / float(n) );
                        html::td_object( html_output, html::FORMATTED, "class", "pink",  NULL, "%5.2f %%", 100.0f * float(unique.R) / float(n) );
                        html::td_object( html_output, html::FORMATTED, "class", "green", NULL, "%5.2f %%", 100.0f * float(unique.L_and_R) / float(n) );
                        html::td_object( html_output, html::FORMATTED,                   NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(unique.L_not_R) / float(n), 100.0f * float(unique_L_not_R_by_mapQ[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, "class", "pink",  NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(unique.R_not_L) / float(n), 100.0f * float(unique_R_not_L_by_mapQ[HI_MAPQ_BIN]) / float(n) );
                    }
                    {
                        html::tr_object tr( html_output, NULL );
                        html::th_object( html_output, html::FORMATTED, NULL, "ambiguous");
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%", 100.0f * float(ambiguous.L) / float(n) );
                        html::td_object( html_output, html::FORMATTED, "class", "pink",  NULL, "%5.2f %%", 100.0f * float(ambiguous.R) / float(n) );
                        html::td_object( html_output, html::FORMATTED, "class", "green", NULL, "%5.2f %%", 100.0f * float(ambiguous.L_and_R) / float(n) );
                        html::td_object( html_output, html::FORMATTED,                   NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(ambiguous.L_not_R) / float(n), 100.0f * float(ambiguous_L_not_R_by_mapQ[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, "class", "pink",  NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(ambiguous.R_not_L) / float(n), 100.0f * float(ambiguous_R_not_L_by_mapQ[HI_MAPQ_BIN]) / float(n) );
                    }
                }
                //
                // discordance stats
                //
                {
                    html::table_object table( html_output, "discordance-stats", "stats", "discordance stats" );
                    {
                        html::tr_object tr( html_output, NULL );
                        html::th_object( html_output, html::FORMATTED, NULL, "" );
                        html::th_object( html_output, html::FORMATTED, NULL, "items" );
                        html::th_object( html_output, html::FORMATTED, NULL, "%% of total" );
                        html::th_object( html_output, html::FORMATTED, NULL, "mate 1" );
                        html::th_object( html_output, html::FORMATTED, NULL, "mate 2" );
                        html::th_object( html_output, html::FORMATTED, NULL, "mate 1 & 2" );
                        html::th_object( html_output, html::FORMATTED, NULL, "unique" );
                        html::th_object( html_output, html::FORMATTED, NULL, "not-ambiguous" );
                    }
                    {
                        html::tr_object tr( html_output, NULL );
                        html::th_object( html_output, html::FORMATTED, NULL, "different reference" );
                        html::td_object( html_output, html::FORMATTED, NULL, "%.2f M    (%.2f M)",             float(cum_different_ref[0]) * 1.0e-6f,            float(cum_different_ref[HI_MAPQ_BIN]) * 1.0e-6f );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_different_ref[0]) / float(n),   100.0f * float(cum_different_ref[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_different_ref1[0]) / float(n),  100.0f * float(cum_different_ref1[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_different_ref2[0]) / float(n),  100.0f * float(cum_different_ref2[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_different_ref12[0]) / float(n), 100.0f * float(cum_different_ref12[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_different_ref_unique[0]) / float(n),          100.0f * float(cum_different_ref_unique[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_different_ref_not_ambiguous[0]) / float(n),   100.0f * float(cum_different_ref_not_ambiguous[HI_MAPQ_BIN]) / float(n) );
                    }
                    {
                        html::tr_object tr( html_output, "class", "alt", NULL );
                        html::th_object( html_output, html::FORMATTED, NULL, "distant" );
                        html::td_object( html_output, html::FORMATTED, NULL, "%.2f M   (%.2f M)",              float(cum_distant[0]) * 1.0e-6f, float(cum_distant[HI_MAPQ_BIN]) * 1.0e-6f );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_distant[0]) / float(n),   100.0f * float(cum_distant[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_distant1[0]) / float(n),  100.0f * float(cum_distant1[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_distant2[0]) / float(n),  100.0f * float(cum_distant2[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_distant12[0]) / float(n), 100.0f * float(cum_distant12[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_distant_unique[0]) / float(n),          100.0f * float(cum_distant_unique[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_distant_not_ambiguous[0]) / float(n),   100.0f * float(cum_distant_not_ambiguous[HI_MAPQ_BIN]) / float(n) );
                    }
                    {
                        html::tr_object tr( html_output, NULL );
                        html::th_object( html_output, html::FORMATTED, NULL, "discordant" );
                        html::td_object( html_output, html::FORMATTED, NULL, "%.2f M   (%.2f M)",              float(cum_discordant[0]) * 1.0e-6f,            float(cum_discordant[HI_MAPQ_BIN]) * 1.0e-6f );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_discordant[0]) / float(n),   100.0f * float(cum_discordant[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_discordant1[0]) / float(n),  100.0f * float(cum_discordant1[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_discordant2[0]) / float(n),  100.0f * float(cum_discordant2[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_discordant12[0]) / float(n), 100.0f * float(cum_discordant12[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_discordant_unique[0]) / float(n),          100.0f * float(cum_discordant_unique[HI_MAPQ_BIN]) / float(n) );
                        html::td_object( html_output, html::FORMATTED, NULL, "%5.2f %%   (%.3f %%)", 100.0f * float(cum_discordant_not_ambiguous[0]) / float(n),   100.0f * float(cum_discordant_not_ambiguous[HI_MAPQ_BIN]) / float(n) );
                    }
                }
                //
                // summary stats
                //
                {
                    html::table_object table( html_output, "summary-stats", "stats", "summary stats" );
                    // ------------------------------------------- generic -------------------------------------------------- //
                    generate_summary_header( html_output );
                    generate_summary_cell( html_output, mapped_bps_name, "mapped bases", al_stats.longer_mapping,   n, "none" );
                    generate_summary_cell( html_output, score_name,      "score",        al_stats.higher_score,     n, "alt" );
                    generate_summary_cell( html_output, ed_name,         "edit distance",al_stats.lower_ed,         n, "none" );
                    generate_summary_cell( html_output, mapQ_name,       "mapQ",         al_stats.higher_mapQ,      n, "alt" );
                    generate_summary_cell( html_output, mms_name,        "mismatches",   al_stats.lower_mms,        n, "none" );
                    generate_summary_cell( html_output, ins_name,        "insertions",   al_stats.lower_ins,        n, "alt" );
                    generate_summary_cell( html_output, dels_name,       "deletions",    al_stats.lower_dels,       n, "none" );
                    generate_summary_cell( html_output, pos_name,        "distance",     al_stats.higher_pos,       n, "alt" );
                    // ------------------------------------------- distant -------------------------------------------------- //
                    generate_summary_header( html_output );
                    generate_summary_cell( html_output, distant_mapped_bps_name, "mapped bases [distant]", distant_stats.longer_mapping,    n, "none" );
                    generate_summary_cell( html_output, distant_ed_name,         "edit distance [distant]",distant_stats.lower_ed,          n, "alt" );
                    generate_summary_cell( html_output, distant_mapQ_name,       "mapQ [distant]",         distant_stats.higher_mapQ,       n, "none" );
                    generate_summary_cell( html_output, distant_mms_name,        "mismatches [distant]",   distant_stats.lower_mms,         n, "alt" );
                    generate_summary_cell( html_output, distant_ins_name,        "insertions [distant]",   distant_stats.lower_ins,         n, "none" );
                    generate_summary_cell( html_output, distant_dels_name,       "deletions [distant]",    distant_stats.lower_dels,        n, "alt" );
                    generate_summary_cell( html_output, distant_pos_name,        "distance [distant]",     distant_stats.higher_pos,        n, "none" );
                    // ------------------------------------------- discordant ---------------------------------------------- //
                    generate_summary_header( html_output );
                    generate_summary_cell( html_output, discordant_mapped_bps_name, "mapped bases [discordant]", discordant_stats.longer_mapping,   n, "alt" );
                    generate_summary_cell( html_output, discordant_ed_name,         "edit distance [discordant]",discordant_stats.lower_ed,         n, "none" );
                    generate_summary_cell( html_output, discordant_mapQ_name,       "mapQ [discordant]",         discordant_stats.higher_mapQ,      n, "alt" );
                    generate_summary_cell( html_output, discordant_mms_name,        "mismatches [discordant]",   discordant_stats.lower_mms,        n, "none" );
                    generate_summary_cell( html_output, discordant_ins_name,        "insertions [discordant]",   discordant_stats.lower_ins,        n, "alt" );
                    generate_summary_cell( html_output, discordant_dels_name,       "deletions [discordant]",    discordant_stats.lower_dels,       n, "none" );
                    generate_summary_cell( html_output, discordant_pos_name,        "distance [discordant]",     discordant_stats.higher_pos,       n, "alt" );
                }

                // paired L not R
                generate_table(
                    html_output,
                    "by mapQ",
                    "paired (L \\ R) vs (R \\ L)",
                    LOG,
                    paired_L_not_R_by_mapQ,
                    paired_R_not_L_by_mapQ,
                    n,
                    n );

                // unique L not R
                generate_table(
                    html_output,
                    "by mapQ",
                    "unique (L \\ R) vs (R \\ L)",
                    LOG,
                    unique_L_not_R_by_mapQ,
                    unique_R_not_L_by_mapQ,
                    n,
                    n );

                // unique L not R
                generate_table(
                    html_output,
                    "by mapQ",
                    "ambiguous (L \\ R) vs (R \\ L)",
                    LOG,
                    ambiguous_L_not_R_by_mapQ,
                    ambiguous_R_not_L_by_mapQ,
                    n,
                    n );

                // different reference by mapQ
                generate_table(
                    html_output,
                    "paired L & R",
                    "different reference by mapQ",
                    LOG,
                    n_different_ref,
                    paired.L_and_R );

                // different reference by mapQ
                generate_table(
                    html_output,
                    "paired L & R",
                    "distant by mapQ",
                    LOG,
                    n_distant,
                    paired.L_and_R );

                // discordant reference by mapQ
                generate_table(
                    html_output,
                    "paired L & R",
                    "discordant by mapQ",
                    LOG,
                    n_discordant,
                    paired.L_and_R );

                // second score by score
                generate_table2d(
                    html_output,
                    "paired L & R",
                    "second score",
                    sec_score_by_score_l,
                    sec_score_by_score_r,
                    "score",
                    log_bins,
                    log_bins,
                    n,
                    paired.L_and_R,
                    false );

                // second score by score
                /*generate_table2d(
                    html_output,
                    "paired L & R",
                    "second ed",
                    sec_ed_by_ed_l,
                    sec_ed_by_ed_r,
                    "ed",
                    ed_bins,
                    ed_bins,
                    n,
                    n,
                    false );*/
            }
        }
    }
    fclose( html_output );
}

} // namespace alndiff
} // namespace nvbio