/**
 * Project     : CASAVA
 * Module      : $RCSfile: ClusterMergerImpl.cpp,v $
 * @author     : Richard Shaw
 * Copyright   : Copyright (c) Illumina 2010. All rights reserved.
 *
 ** 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).
 *
 */

/*****************************************************************************/

#include <cassert>
#include <fstream>
#include <iterator>

#include "common/Exceptions.hh"
#include "common/SequenceUtils.hh"

#include "variance/ClusterMergerImpl.hh"

// #define DEBUG_CLM 1

/*****************************************************************************/

namespace ca {
namespace variance_detection {

/*****************************************************************************/

using namespace std;
using namespace casava::common;

/*****************************************************************************/

ClusterMergerImpl::ClusterMergerImpl(const std::string& inputFileName,
                                     const std::string& outputFileName,
                                     bool doAdjMerge,
                                     int maxAdjMergeInterBreakPtDist,
                                     unsigned int maxLinksPerCluster)
    : myInputFileName(inputFileName),
      myOutputFileName(outputFileName),
      myDoAdjMerge(doAdjMerge),
      myMaxAdjMergeInterBreakPtDist(maxAdjMergeInterBreakPtDist),
      myMaxLinksPerCluster(maxLinksPerCluster)
{
    if (myInputFileName.empty()) {
        BOOST_THROW_EXCEPTION(CasavaException(EINVAL,
                                              "No value given for "
                                              "input file"));
    }

    if (myOutputFileName.empty()) {
        BOOST_THROW_EXCEPTION(CasavaException(EINVAL,
                                              "No value given for "
                                              "output file"));
    }
}

/*****************************************************************************/

void ClusterMergerImpl::importClusters(PairStats& pairStats)
{
    ifstream inputFile(myInputFileName.c_str());

    if (!inputFile.is_open()) {
        BOOST_THROW_EXCEPTION(CasavaException(EINVAL,
                                              "Failed to open input file: "
                                              + myInputFileName));
    }

    unsigned int numGroups(0);

    while (!inputFile.eof()) {
        AnomReadVec anomReadVec(++numGroups);
        inputFile >> anomReadVec;

        if (!anomReadVec.empty()) {
            anomReadVec.updateBreakPoints(pairStats);
            myAnomReadVecVec.push_back(anomReadVec);
        }
    }

    inputFile.close();
}

/*****************************************************************************/

void ClusterMergerImpl::mergeClusters()
{
    std::string label;

    if (myDoAdjMerge) {
        // Compatible adjacency merging

        // Sort by left breakpoint and look for adjacency merge candidates.
        findAdjacencyMergeCandidates(true, myAdjacencyMergeCandMap);

        // Sort by left breakpoint and look for adjacency merge candidates.
        findAdjacencyMergeCandidates(false, myAdjacencyMergeCandMap);

        // Dump original adjacency merge candidate info.
        label = "Original adjacency merge candidates";
        dumpMergeCandidates(label, myAdjacencyMergeCandMap);
        dumpMergeCandidateStats(label, myAdjacencyMergeCandMap);

        // Consolidate adjacency merges
        consolidateMerges(myAdjacencyMergeCandMap);

        // Dump post-consolidation adjacency merge candidate info.
        label = "Post-consolidation adjacency merge candidates";
#ifdef DEBUG_CLM
        dumpMergeCandidates(label, myAdjacencyMergeCandMap);
#endif // DEBUG_CLM
        dumpMergeCandidateStats(label, myAdjacencyMergeCandMap);

        // Do the merging by adjacency.
        doMerges(myAdjacencyMergeCandMap);
    }


    // Shared read merging

    buildReadClusterIdsMap(myReadClusterIdsMap);
    dumpReadClusterIdsMapStats(myReadClusterIdsMap);

    // Look for shared read merge candidates.
    findSharedReadMergeCandidates(myReadClusterIdsMap,
                                  mySharedReadMergeCandMap);
    killSuperHubClusters(mySharedReadMergeCandMap);

    // Dump original shared read merge candidate info.
    label = "Original shared read merge candidates";
#ifdef DEBUG_CLM
    dumpMergeCandidates(label, mySharedReadMergeCandMap);
#endif // DEBUG_CLM
    dumpMergeCandidateStats(label, mySharedReadMergeCandMap);

    // Consolidate shared read merges
    consolidateMerges(mySharedReadMergeCandMap);

    // Dump post-consolidation shared read merge candidate info.
    label = "Post-consolidation shared read merge candidates";
#ifdef DEBUG_CLM
    dumpMergeCandidates(label, mySharedReadMergeCandMap);
#endif // DEBUG_CLM
    dumpMergeCandidateStats(label, mySharedReadMergeCandMap);

    // Do the merging by shared reads
    doMerges(mySharedReadMergeCandMap);

    // Confirm that there are no duplicated reads remaining.
    ReadClusterIdsMap postMergeReadClusterIdsMap;
    std::cout << "Post-merge read cluster stats :-" << std::endl;
    buildReadClusterIdsMap(postMergeReadClusterIdsMap);
    dumpReadClusterIdsMapStats(postMergeReadClusterIdsMap);
}

/*****************************************************************************/

void ClusterMergerImpl::generateOutputFile()
{
    ofstream outputFile(myOutputFileName.c_str());

    if (!outputFile.is_open()) {
        BOOST_THROW_EXCEPTION(CasavaException(EINVAL,
                                              "Failed to open output file: "
                                              + myOutputFileName));
    }

    for (AnomReadVecVecIter readVecIter(myAnomReadVecVec.begin());
         readVecIter != myAnomReadVecVec.end(); ++readVecIter) {

        if (myKilledClusterIdSet.find(readVecIter->myGroupNum)
            == myKilledClusterIdSet.end()) {
            // Output the cluster as long as it has not been killed.
            outputFile << *readVecIter;
        }
    }

    outputFile.close();
}

/*****************************************************************************/

void ClusterMergerImpl::buildReadClusterIdsMap(ReadClusterIdsMap&
                                               readClusterIdsMap)
{
    for (AnomReadVecVecIter readVecIter(myAnomReadVecVec.begin());
         readVecIter != myAnomReadVecVec.end(); ++readVecIter) {
        AnomReadVec& anomReadVec(*readVecIter);

        for (AnomReadVecCIter anomReadCIter(anomReadVec.begin());
             anomReadCIter != anomReadVec.end(); ++anomReadCIter) {
            const Singleton* anomReadPtr(*anomReadCIter);
            const std::string readId(SequenceUtils::getReadId(anomReadPtr->
                                                              alignment_));
            readClusterIdsMap[readId].push_back(anomReadVec.myGroupNum);
        }
    }
}

/*****************************************************************************/

void ClusterMergerImpl::dumpReadClusterIdsMapStats(ReadClusterIdsMap&
                                                   readClusterIdsMap)
{
    unsigned int numSingleClusterReads(0);
    unsigned int numTwoClusterReads(0);
    unsigned int numGtTwoClusterReads(0);

    for (ReadClusterIdsMapCIter readClustersCIter(readClusterIdsMap.begin());
         readClustersCIter != readClusterIdsMap.end(); ++readClustersCIter) {
        const unsigned int numClusters(readClustersCIter->second.size());

        switch (numClusters) {
        case 1:
            ++numSingleClusterReads;
            break;
        case 2:
            ++numTwoClusterReads;
#ifdef DEBUG_CLM
            {
                std::cout << readClustersCIter->first << " :";
                const ClusterIdVec& clusterIdVec(readClustersCIter->second);
                std::ostream_iterator<ClusterId> id_ostrm(std::cout, " ");
                copy(clusterIdVec.begin(), clusterIdVec.end(), id_ostrm);
                std::cout << std::endl;
            }
#endif // DEBUG_CLM
            break;
        default:
            ++numGtTwoClusterReads;
            break;
        }
    }

    std::cout << readClusterIdsMap.size() << " reads :-" << std::endl
              << "In  1 cluster  : " << numSingleClusterReads << std::endl
              << "In  2 clusters : " << numTwoClusterReads << std::endl
              << "In >2 clusters : " << numGtTwoClusterReads << std::endl;
}

/*****************************************************************************/

void
ClusterMergerImpl::findSharedReadMergeCandidates(ReadClusterIdsMap&
                                                 readClusterIdsMap,
                                                 MergeCandMap& mergeCandMap)
{
    // Link clusters that share reads.
    for (ReadClusterIdsMapCIter readClustersCIter(readClusterIdsMap.begin());
         readClustersCIter != readClusterIdsMap.end(); ++readClustersCIter) {
        ClusterIdVec clusterIdVec(readClustersCIter->second);

        if (clusterIdVec.size() > 1) {
            // Key by Anom Pair cluster ID -> extract it.
            ClusterId anomPairClusterId(0);

            for (ClusterIdVecIter clusterIdIter(clusterIdVec.begin());
                 clusterIdIter != clusterIdVec.end(); ++clusterIdIter) {
                const ClusterId clusterId(*clusterIdIter);

                assert(myAnomReadVecVec[clusterId - 1].myGroupNum
                       == clusterId);

                if (myAnomReadVecVec[clusterId - 1].myAnomReadType
                    != AnomalousRead::ShadowOrSemiAligned) {
                    anomPairClusterId = clusterId;
                    clusterIdVec.erase(clusterIdIter);
                    break;
                }
            }

            // Now add the other cluster ID - set prevents duplicate storage.
            ClusterIdSet& clusterIdSet(mergeCandMap[anomPairClusterId]);
            copy(clusterIdVec.begin(), clusterIdVec.end(),
                 inserter(clusterIdSet, clusterIdSet.begin()));

#ifdef DEBUG_CLM
            AnomReadVec&
                absorberCluster(myAnomReadVecVec[anomPairClusterId - 1]);

            for (ClusterIdVecCIter
                     absorbeeClusterIdCIter(clusterIdVec.begin());
                 absorbeeClusterIdCIter != clusterIdVec.end();
                 ++absorbeeClusterIdCIter) {
                ClusterId absorbeeClusterId(*absorbeeClusterIdCIter);
                const AnomReadVec&
                    absorbeeCluster(myAnomReadVecVec[absorbeeClusterId - 1]);

                const unsigned int
                    numCommonReads(absorberCluster.
                                   numCommonReads(absorbeeCluster));

                cout << "Absorber " << absorberCluster.myAnomReadType
                     << "[" << anomPairClusterId << "] "
                     << " absorbee " << absorbeeCluster.myAnomReadType
                     << "[" << absorbeeClusterId << "] "
                     << " numCommonReads " << numCommonReads << endl;
            }
#endif // DEBUG_CLM
        }
    }
}

/*****************************************************************************/

void ClusterMergerImpl::killSuperHubClusters(MergeCandMap& mergeCandMap)
{
    for (MergeCandMapIter pairIter(mergeCandMap.begin());
         pairIter != mergeCandMap.end(); ) {
        const ClusterId clusterId(pairIter->first);
        const unsigned int numLinks(pairIter->second.size());


        if (numLinks > myMaxLinksPerCluster) {
            cout << "Killing Super Hub Cluster " << clusterId << " of type "
                 << myAnomReadVecVec[clusterId - 1].myAnomReadType
                 << " with " << numLinks << " links to other clusters and "
                 << myAnomReadVecVec[clusterId - 1].size() << " reads"
                 << std::endl;

            // Remove the cluster from the key side of the map & record it.
            // May still be within lists of clusters on the value side.
            mergeCandMap.erase(pairIter++);
            myKilledClusterIdSet.insert(clusterId);
        } else {
            ++pairIter;
        }
    }

    // Second pass to remove the cluster from any sets on the value side.
    for (MergeCandMapIter pairIter(mergeCandMap.begin());
         pairIter != mergeCandMap.end(); ) {

        for (ClusterIdSetCIter
                 killedClusterIdCiter(myKilledClusterIdSet.begin());
             killedClusterIdCiter != myKilledClusterIdSet.end();
             ++killedClusterIdCiter) {
            pairIter->second.erase(*killedClusterIdCiter);
        }

        // If there are no longer any clusters in the set, remove the entry.
        if (pairIter->second.empty()) {
            const ClusterId clusterId(pairIter->first);
            std::cout << "After removing super hub clusters, cluster "
                      << clusterId << " has no links left. Removing entry"
                      << std::cout;
            mergeCandMap.erase(pairIter++);
        } else {
            ++pairIter;
        }
    }
}

/*****************************************************************************/

void
ClusterMergerImpl::findAdjacencyMergeCandidates(bool leftNotRight,
                                                MergeCandMap& mergeCandMap)
{
    if (leftNotRight) {
        sort(myAnomReadVecVec.begin(), myAnomReadVecVec.end(),
             SortAnomReadVecByLeftBrk());
    } else {
        sort(myAnomReadVecVec.begin(), myAnomReadVecVec.end(),
             SortAnomReadVecByRightBrk());
    }

    AnomReadVecVecIter lastReadVecIter(myAnomReadVecVec.begin());
    ClusterId lastClusterId(0);
    ClusterId anomPairClusterId(0);
    ClusterId semiClusterId(0);
    AnomalousRead::Type lastAnomReadType(AnomalousRead::ShadowOrSemiAligned);
    int lastBrkPtPos(-1000000);

    for (AnomReadVecVecIter readVecIter(myAnomReadVecVec.begin());
         readVecIter != myAnomReadVecVec.end(); ++readVecIter) {
        bool isMergeCand(false);

        const ClusterId clusterId(readVecIter->myGroupNum);
        const AnomalousRead::Type anomReadType(readVecIter->myAnomReadType);
        const int brkPtPos(leftNotRight
                           ? readVecIter->myLeftBreakPos
                           : readVecIter->myRightBreakPos);

        if (lastClusterId != 0) { // Two clusters to consider?
            const int separ(brkPtPos - lastBrkPtPos);

            if ((anomReadType == AnomalousRead::ShadowOrSemiAligned)
                && (lastAnomReadType != AnomalousRead::ShadowOrSemiAligned)) {
                isMergeCand = true;
                anomPairClusterId = lastClusterId;
                semiClusterId = clusterId;
            } else if ((anomReadType != AnomalousRead::ShadowOrSemiAligned)
                       && (lastAnomReadType
                           == AnomalousRead::ShadowOrSemiAligned)) {
                isMergeCand = true;
                anomPairClusterId = clusterId;
                semiClusterId = lastClusterId;
            }

            isMergeCand = (isMergeCand
                           && (separ <= myMaxAdjMergeInterBreakPtDist));

            if (isMergeCand) {
                cerr << "Found "
                     << lastAnomReadType << "[" << lastClusterId << "]"
                     << " " << anomReadType << "[" << clusterId << "]"
                     << " merge cand : separ " << separ << endl;

                const unsigned int
                    numCommonReads(readVecIter->
                                   numCommonReads(*lastReadVecIter));

                cerr << "numCommonReads : " << numCommonReads << endl;
            }

            if (isMergeCand) {
                mergeCandMap[anomPairClusterId].insert(semiClusterId);
            }
        }

        lastReadVecIter = readVecIter;
        lastClusterId = clusterId;
        lastAnomReadType = anomReadType;
        lastBrkPtPos = brkPtPos;
    }
}

/*****************************************************************************/

void ClusterMergerImpl::makeReverseMergeMap(MergeCandMap& mergeCandMap,
                                            MergeCandMap& reverseMergeMap)
{
    for (MergeCandMapCIter mergeCandCIter(mergeCandMap.begin());
         mergeCandCIter != mergeCandMap.end(); ++mergeCandCIter) {
        const ClusterId absorberClusterId(mergeCandCIter->first);
        const ClusterIdSet absorbeeClusterIdSet(mergeCandCIter->second);

        for (ClusterIdSetCIter clusterIdCIter(absorbeeClusterIdSet.begin());
             clusterIdCIter != absorbeeClusterIdSet.end(); ++clusterIdCIter) {
            const ClusterId absorbeeClusterId(*clusterIdCIter);
            reverseMergeMap[absorbeeClusterId].insert(absorberClusterId);
        }
    }
}

/*****************************************************************************/

void ClusterMergerImpl::absorbOneMerge(MergeCandMap& mergeCandMap,
                                       MergeCandMap& reverseMergeMap,
                                       AbsorptionMap& absorbedByMap,
                                       ClusterId absorberClusterId,
                                       ClusterId absorbeeClusterId)
{
#ifdef DEBUG_CLM
    cout << "absorbOneMerge : absorber " << absorberClusterId
         << ", absorbee " << absorbeeClusterId << std::endl;
#endif // DEBUG_CLM

    // Sanity check : nothing to do if cluster trying to absorb itself.

    if (absorberClusterId == absorbeeClusterId) {
#ifdef DEBUG_CLM
        cout << "absorbOneMerge : absorber " << absorberClusterId
             << " trying to absorb itself - nothing to do" << std::endl;
#endif // DEBUG_CLM
        return;
    }

    // Check if specified absorber has itself already been absorbed.
    // If so, recurse.

    AbsorptionMapCIter absorbedByCIter(absorbedByMap.find(absorberClusterId));

    if (absorbedByCIter != absorbedByMap.end()) {
        ClusterId newAbsorberClusterId(absorbedByCIter->second);

#ifdef DEBUG_CLM
        cout << "Absorber " << absorberClusterId << " previously absorbed by "
             << newAbsorberClusterId << " : recursing for absorption of "
             << absorbeeClusterId << std::endl;
#endif // DEBUG_CLM

        absorbOneMerge(mergeCandMap, reverseMergeMap, absorbedByMap,
                       newAbsorberClusterId, absorbeeClusterId);
        return;
    }

    ClusterIdSet& absorberClusterIdSet(mergeCandMap[absorberClusterId]);

    // Has the absorber previously absorbed the absorbee?
    if (absorberClusterIdSet.find(absorbeeClusterId)
        != absorberClusterIdSet.end()) {
#ifdef DEBUG_CLM
        cout << "Absorber " << absorberClusterId << " had previously absorbed "
             << absorbeeClusterId << " : nothing to do " << std::endl;
#endif // DEBUG_CLM
        return;
    }

    // Has the absorbee already been absorbed by another absorber?
    absorbedByCIter = absorbedByMap.find(absorbeeClusterId);

    if (absorbedByCIter != absorbedByMap.end()) {
        // Yes, so now the absorber becomes the absorbee of the previous
        // absorber.
        ClusterId newAbsorberClusterId(absorbedByCIter->second);

#ifdef DEBUG_CLM
        cout << "Absorbee " << absorbeeClusterId << " previously absorbed by "
             << newAbsorberClusterId << " so absorber "
             << absorberClusterId << " must now be absorbed by it instead"
             << std::endl;
#endif // DEBUG_CLM

        absorbOneMerge(mergeCandMap, reverseMergeMap, absorbedByMap,
                       newAbsorberClusterId, absorberClusterId);
        return;
    }


    // Find the absorbee and its absorbees.
    MergeCandMapIter absorbeeIter(mergeCandMap.find(absorbeeClusterId));

    if (absorbeeIter == mergeCandMap.end()) {
#ifdef DEBUG_CLM
        cout << "Absorbee " << absorbeeClusterId
             << " : previously absorbed so " << absorberClusterId
             << " cannot absorb it" << std::endl;
#endif // DEBUG_CLM
        return;
    }


    // Absorb the cluster ID itself.
    absorberClusterIdSet.insert(absorbeeClusterId);

    // Absorb all the cluster IDs that absorbeeClusterId was going to absorb.
    ClusterIdSet& absorbeeClusterIdSet(absorbeeIter->second);
    copy(absorbeeClusterIdSet.begin(), absorbeeClusterIdSet.end(),
         inserter(absorberClusterIdSet, absorberClusterIdSet.begin()));

#ifdef DEBUG_CLM
    cout << "Absorber " << absorberClusterId << " absorbed "
         << absorbeeClusterId << " and all its absorbees : ";
    std::ostream_iterator<ClusterId> id_ostrm(std::cout, " ");
    copy(absorbeeClusterIdSet.begin(), absorbeeClusterIdSet.end(),
                 id_ostrm);
    cout << std::endl;
#endif // DEBUG_CLM

    // Erase the absorbee entry
    mergeCandMap.erase(absorbeeIter);

    // Update the absorbed-by map
    absorbedByMap[absorbeeClusterId] = absorberClusterId;
}

/*****************************************************************************/
// Find out which candidate merges are linked by overlapping membership.

void ClusterMergerImpl::consolidateMerges(MergeCandMap& mergeCandMap)
{
    MergeCandMap reverseMergeMap;
    makeReverseMergeMap(mergeCandMap, reverseMergeMap);


    AbsorptionMap absorbedByMap;

    for (MergeCandMapCIter revMergeCandCIter(reverseMergeMap.begin());
         revMergeCandCIter != reverseMergeMap.end(); ++revMergeCandCIter) {
        ClusterIdSet absorberClusterIdSet(revMergeCandCIter->second);

        if (absorberClusterIdSet.size() > 1) {
#ifdef DEBUG_CLM
            const ClusterId absorbeeClusterId(revMergeCandCIter->first);
            std::cout << "Cluster : " << absorbeeClusterId
                      << " is to be merged into the following clusters :-"
                      << std::endl;
#endif // DEBUG_CLM

            // Extract the first one as the super absorber.
            const ClusterId
                superAbsorberClusterId(*absorberClusterIdSet.begin());
            absorberClusterIdSet.erase(absorberClusterIdSet.begin());

#ifdef DEBUG_CLM
            std::cout << "  Super Absorber Cluster ID : "
                      << superAbsorberClusterId << std::endl
                      << "  Other absorbers :";
            std::ostream_iterator<ClusterId> id_ostrm(std::cout, " ");
            copy(absorberClusterIdSet.begin(), absorberClusterIdSet.end(),
                 id_ostrm);
            std::cout << std::endl;
#endif // DEBUG_CLM

            for (ClusterIdSetCIter
                     absorberClusterIdCIter(absorberClusterIdSet.begin());
                 absorberClusterIdCIter != absorberClusterIdSet.end();
                 ++absorberClusterIdCIter) {
                const ClusterId absorberClusterId(*absorberClusterIdCIter);

                // SuperAbsorber gets Absorber (unless SuperAbsorber itself has
                // already been absorbed - handled within absorbOneMerge).
                absorbOneMerge(mergeCandMap, reverseMergeMap,
                               absorbedByMap,
                               superAbsorberClusterId,
                               absorberClusterId);
            }

#ifdef DEBUG_CLM
            std::cout << std::endl;
#endif // DEBUG_CLM
        }
    }

    // DEBUG
    MergeCandMap postReverseMergeMap;
    makeReverseMergeMap(mergeCandMap, postReverseMergeMap);
    unsigned int numFailures(0);

    for (MergeCandMapCIter revMergeCandCIter(postReverseMergeMap.begin());
         revMergeCandCIter != postReverseMergeMap.end(); ++revMergeCandCIter) {
        ClusterIdSet absorberClusterIdSet(revMergeCandCIter->second);

        if (absorberClusterIdSet.size() > 1) {
            std::cout << "*** CONSOLIDATION FAILURE FOR ABSORBEE "
                      << revMergeCandCIter->first << " : Absorbers : ";
            std::ostream_iterator<ClusterId> id_ostrm(std::cout, " ");
            copy(absorberClusterIdSet.begin(), absorberClusterIdSet.end(),
                 id_ostrm);
            std::cout << std::endl;
            ++numFailures;
        }
    }

    if (numFailures == 0) {
        std::cout << "CONSOLIDATION OK" << std::endl;
    }
}

/*****************************************************************************/

void
ClusterMergerImpl::dumpMergeCandidates(const std::string& label,
                                       const MergeCandMap& mergeCandMap)
{
    std::cout << "Linked clusters for " << label << " :-" << std::endl;

    for (MergeCandMapCIter mergeCandCIter(mergeCandMap.begin());
         mergeCandCIter != mergeCandMap.end(); ++mergeCandCIter) {
        const ClusterId absorberClusterId(mergeCandCIter->first);
        const AnomReadVec&
            absorberCluster(myAnomReadVecVec[absorberClusterId - 1]);
        std::cout << absorberClusterId << " "
                  << absorberCluster.myAnomReadType << " ["
                  << absorberCluster.myLeftBreakPos << ","
                  << absorberCluster.myRightBreakPos << "]"
                  << " := ";

        const ClusterIdSet absorbeeClusterIdSet(mergeCandCIter->second);

        for (ClusterIdSetCIter clusterIdCIter(absorbeeClusterIdSet.begin());
             clusterIdCIter != absorbeeClusterIdSet.end(); ++clusterIdCIter) {
            const ClusterId absorbeeClusterId(*clusterIdCIter);
            const AnomReadVec&
                absorbeeCluster(myAnomReadVecVec[absorbeeClusterId - 1]);
            std::cout << " " << absorbeeClusterId << " "
                      << absorbeeCluster.myAnomReadType << " ["
                      << absorbeeCluster.myLeftBreakPos << ","
                      << absorbeeCluster.myRightBreakPos << "]";
        }

        std::cout << std::endl;
    }
}

/*****************************************************************************/

void
ClusterMergerImpl::dumpMergeCandidateStats(const std::string& label,
                                           const MergeCandMap& mergeCandMap)
{
    // Need to put clusters back in GroupNum order so can use GroupNum - 1
    // as index.
    sort(myAnomReadVecVec.begin(), myAnomReadVecVec.end(),
         SortAnomReadVecByGroupNum());

    typedef std::map<unsigned int, unsigned int> Counts;
    typedef Counts::const_iterator CountsCIter;

    typedef std::map<AnomalousRead::Type, Counts> MergeTypeCountsMap;
    typedef MergeTypeCountsMap::const_iterator MergeTypeCountsMapCIter;
    MergeTypeCountsMap mergeTypeCountsMap;

    for (MergeCandMapCIter mergeCandCIter(mergeCandMap.begin());
         mergeCandCIter != mergeCandMap.end(); ++mergeCandCIter) {
        const ClusterId absorberClusterId(mergeCandCIter->first);
        const AnomReadVec&
            absorberCluster(myAnomReadVecVec[absorberClusterId - 1]);
        AnomalousRead::Type absorberType(absorberCluster.myAnomReadType);

        const ClusterIdSet absorbeeClusterIdSet(mergeCandCIter->second);
        const unsigned int numAbsorbees(absorbeeClusterIdSet.size());

        ++mergeTypeCountsMap[absorberType][numAbsorbees];
    }

    std::cout << "Stats for " << label << std::endl;

    for (MergeTypeCountsMapCIter
             mergeTypeCountsCIter(mergeTypeCountsMap.begin());
         mergeTypeCountsCIter != mergeTypeCountsMap.end();
         ++mergeTypeCountsCIter) {
        std::cout << mergeTypeCountsCIter->first << " :";
        const Counts& counts(mergeTypeCountsCIter->second);

        for (CountsCIter countsCIter(counts.begin());
             countsCIter != counts.end(); ++countsCIter) {
            std::cout << " " << countsCIter->first << ":"
                      << countsCIter->second;
        }

        std::cout << std::endl;
    }
}

/*****************************************************************************/

void ClusterMergerImpl::doMerges(MergeCandMap& mergeCandMap)
{
    // Need to put clusters back in GroupNum order so can use GroupNum - 1
    // as index.
    sort(myAnomReadVecVec.begin(), myAnomReadVecVec.end(),
         SortAnomReadVecByGroupNum());

    unsigned int numClusterAbsorptionAttempts(0);
    unsigned int numClustersAbsorbed(0);

    for (MergeCandMapIter mergeCandIter(mergeCandMap.begin());
         mergeCandIter != mergeCandMap.end(); ++mergeCandIter) {
        ClusterId absorberClusterId(mergeCandIter->first);
        AnomReadVec& absorberCluster(myAnomReadVecVec[absorberClusterId - 1]);
        const ClusterIdSet& absorbeeClusterIdSet(mergeCandIter->second);

        for (ClusterIdSetCIter clusterIdCIter(absorbeeClusterIdSet.begin());
             clusterIdCIter != absorbeeClusterIdSet.end(); ++clusterIdCIter) {
            ++numClusterAbsorptionAttempts;
            ClusterId absorbeeClusterId(*clusterIdCIter);
            AnomReadVec&
                absorbeeCluster(myAnomReadVecVec[absorbeeClusterId - 1]);

            if (!absorbeeCluster.empty()) {
                absorberCluster.absorb(absorbeeCluster);
                ++numClustersAbsorbed;
#ifdef DEBUG_CLM
                cerr << "Cluster " << absorberClusterId
                     << " absorbed " << absorbeeClusterId << std::endl;
#endif // DEBUG_CLM
            }
        }
    }

    cerr << "Num cluster absorption attempts : "
         << numClusterAbsorptionAttempts << std::endl
         << "Num clusters absorbed : " << numClustersAbsorbed << std::endl;
}

/*****************************************************************************/

} // end namespace variance_detection
} // end namespace ca

/*****************************************************************************/