xref: /dports/biology/lamarc/lamarc-2.1.8/src/tree/rectree.cpp

// $Id: rectree.cpp,v 1.77 2012/02/16 19:10:59 bobgian Exp $

/*
  Copyright 2002  Peter Beerli, Mary Kuhner, Jon Yamato and Joseph Felsenstein

  This software is distributed free of charge for non-commercial use
  and is copyrighted.  Of course, we do not guarantee that the software
  works, and are not responsible for any damage you may cause or have.
*/

#include <cassert>

#include "local_build.h"

#include "errhandling.h"                // Can throw overrun_error.
#include "force.h"                      // For TimeSize object.
#include "range.h"
#include "runreport.h"
#include "tree.h"
#include "mathx.h"
#include "branchbuffer.h"               // Used in SetTargetLinksFrom for BranchBuffer::ExtractConstBranches.
#include "fc_status.h"                  // Used in Prune() to determine what recombinations can be removed.

#ifdef DMALLOC_FUNC_CHECK
#include "/usr/local/include/dmalloc.h"
#endif

using namespace std;

// Note that some functions/variables defined here have "_Tr" (for "Tree") as a suffix in their name.
// This is to distinguish them from functions/variables in class Range (or RecRange) of the same name, which use "_Rg".

//------------------------------------------------------------------------------------

typedef boost::shared_ptr<RBranch> RBranch_ptr;

//------------------------------------------------------------------------------------

RecTree::RecTree()
    : Tree(),
      m_pMovingLocusVec(),
      m_protoMovingCells()
{
    m_numTargetLinks_Tr = 0;
    m_numNewTargetLinks_Tr = 0;
}

//------------------------------------------------------------------------------------

RecTree::RecTree(const RecTree & tree, bool makestump)
    : Tree(tree, makestump),
      m_pMovingLocusVec(tree.m_pMovingLocusVec),
      m_protoMovingCells(tree.m_protoMovingCells)

{
    m_numTargetLinks_Tr = 0;
    m_numNewTargetLinks_Tr = 0;
}

//------------------------------------------------------------------------------------

Tree * RecTree::Clone() const
{
    RecTree * tree = new RecTree(*this, false);
    return tree;
} // RecTree::Clone

//------------------------------------------------------------------------------------

Tree * RecTree::MakeStump() const
{
    RecTree * tree = new RecTree(*this, true);
    return tree;
}

//------------------------------------------------------------------------------------

void RecTree::Clear()
{
    Tree::Clear();
    m_numTargetLinks_Tr = 0;
    m_numNewTargetLinks_Tr = 0;
}

//------------------------------------------------------------------------------------

void RecTree::CopyTips(const Tree * tree)
{
    Tree::CopyTips(tree);
    m_numTargetLinks_Tr = 0;
    m_numNewTargetLinks_Tr = 0;
    for(vector<Locus>::const_iterator locus = m_pMovingLocusVec->begin(); locus != m_pMovingLocusVec->end(); ++locus)
    {
        m_aliases.push_back(locus->GetDLCalc()->RecalculateAliases(*this, *locus));
    }
}

//------------------------------------------------------------------------------------

void RecTree::CopyBody(const Tree * tree)
{
    Tree::CopyBody(tree);
    m_numTargetLinks_Tr = 0;
    m_numNewTargetLinks_Tr = 0;
}

//------------------------------------------------------------------------------------

void RecTree::CopyPartialBody(const Tree * tree)
{
    Tree::CopyPartialBody(tree);
    m_numTargetLinks_Tr = 0;
    m_numNewTargetLinks_Tr = 0;
}

//------------------------------------------------------------------------------------

void RecTree::Break(Branch_ptr pBranch)
{
    Branch_ptr pParent = pBranch->Parent(0);

    if (pParent->CanRemove(pBranch))
    {
        Break(pParent);
        m_timeList.Remove(pParent);

        pParent = pBranch->Parent(1);
        if (pParent)
        {
            Break(pParent);
            m_timeList.Remove(pParent);
        }
    }
}

//------------------------------------------------------------------------------------

const vector<LocusCell> & RecTree::CollectMovingCells()
{
    if (m_protoMovingCells.empty())
    {
        unsigned long i;
        for (i = 0; i < m_pMovingLocusVec->size(); ++i)
        {
            m_protoMovingCells.push_back((*m_pMovingLocusVec)[i].GetProtoCell());
        }
    }
    return m_protoMovingCells;
} // CollectCells

//------------------------------------------------------------------------------------

vector<Branch_ptr> RecTree::ActivateTips(Tree * othertree)
{
    vector<Branch_ptr> tips = Tree::ActivateTips(othertree);
    m_numTargetLinks_Tr = m_timeList.GetNTips() * (m_totalSites - 1); // -1 because there is no link after last site.
    m_numNewTargetLinks_Tr = 0;
    return tips;
}

//------------------------------------------------------------------------------------

Branch_ptr RecTree::ActivateBranch(Tree * othertree)
{
    Branch_ptr pActive = Tree::ActivateBranch(othertree);

    m_numTargetLinks_Tr += pActive->m_rangePtr->NumTargetLinks_Rg();
    return pActive;
}

//------------------------------------------------------------------------------------

Branch_ptr RecTree::ActivateRoot(FC_Status & fcstatus)
{
    Branch_ptr pRoot = Tree::ActivateRoot(fcstatus);
    m_numTargetLinks_Tr += pRoot->m_rangePtr->NumTargetLinks_Rg();
    m_numNewTargetLinks_Tr = 0;
    return pRoot;
}

//------------------------------------------------------------------------------------

void RecTree::AttachBase(Branch_ptr newroot)
{
    Tree::AttachBase(newroot);
    m_numTargetLinks_Tr = 0;
}

//------------------------------------------------------------------------------------

vector<Branch_ptr> RecTree::FirstInterval(double eventT)
{
    Branch_ptr pBranch;

    vector<Branch_ptr> newinactives = Tree::FirstInterval(eventT);

    m_numNewTargetLinks_Tr = 0;         // Initialize the count of open links.

    return newinactives;
}

//------------------------------------------------------------------------------------

void RecTree::NextInterval(Branch_ptr pBranch)
{
    Branch_ptr pParent = pBranch->Parent(0);

    m_numNewTargetLinks_Tr -= pParent->Child(0)->m_rangePtr->NumNewTargetLinks_Rg();

    // If the branch has a sibling, we must remove it as well (it's a coalescence event replaced by its parent).
    if (pParent->Child(1))
    {
        m_numNewTargetLinks_Tr -= pParent->Child(1)->m_rangePtr->NumNewTargetLinks_Rg();
    }
    else if (pBranch->Parent(1))        // If the branch has a second parent, insert it.
    {
        m_numNewTargetLinks_Tr += pBranch->Parent(1)->m_rangePtr->NumNewTargetLinks_Rg();
    }

    m_numNewTargetLinks_Tr += pParent->m_rangePtr->NumNewTargetLinks_Rg();
}

//------------------------------------------------------------------------------------

void RecTree::Prune()
{
    // This routine DEMANDS that all current-state Range information in the tree is correct!  Thus:
    assert(m_timeList.RevalidateAllRanges());

    // This routine locates recombinations which are no longer relevant due to occuring at links that
    // are no longer targettable.  It removes such recombinations and everything logically dependent
    // on them.  It also makes all "old state" information in the Ranges equal to the current
    // state, in preparation for the next cycle.

    m_numNewTargetLinks_Tr = 0;
    Branch_ptr pBranch, pChild;
    // The original cut branch never needs to be updated or removed by Prune().
    Branchiter start(m_timeList.NextNonTimeTiedBranch(m_firstInvalid));

    FC_Status fcstatus;
    vector<Branch_ptr> branches(FindBranchesImmediatelyTipwardOf(start));
    vector<Branch_ptr>::iterator branch;

#if LAMARC_QA_FC_ON
    for(branch = branches.begin(); branch != branches.end(); ++branch)
    {
        fcstatus.Increment_FC_Counts((*branch)->GetLiveSites_Br());
    }
#endif

    // Loop over pBranch starting from cut point.
    Branchiter brit;
    for (brit = start; brit != m_timeList.EndBranch(); /* increments inside loop */ )
    {
        pBranch = *brit;

#if LAMARC_QA_FC_ON
        // If we're a fully functional, two-legged coalescence ...
        if (pBranch->Event() == btypeCoal && !pBranch->m_marked)
        {
            rangeset coalesced_sites =
                Intersection(pBranch->Child(0)->GetLiveSites_Br(), pBranch->Child(1)->GetLiveSites_Br());
            fcstatus.Decrement_FC_Counts(coalesced_sites);
        }
#endif

        rangeset fcsites;
#if LAMARC_QA_FC_ON
        fcsites = fcstatus.Coalesced_Sites();
#endif

        // Update this branch's oldtargetsites info, which is no longer needed in this cycle,
        // so that it will be right for next cycle.  Done here because it requires FC information.
        pBranch->ResetOldTargetSites_Br(fcsites);

        if (pBranch->IsRemovableRecombinationLeg(fcsites))
        {
#ifndef NDEBUG
            // validation check:
            Range * childrangeptr = pBranch->Child(0)->m_rangePtr;
            rangeset childtargetsites(Union(childrangeptr->GetDiseaseSites_Rg(),
                                            RemoveRangeFromRange(fcsites, childrangeptr->GetLiveSites_Rg())));
            assert(TargetLinksFrom(childtargetsites) == childrangeptr->GetTargetLinks_Rg());
#endif
            Branch_ptr pSpouse = pBranch->GetRecPartner();
            Branch_ptr pRemove(pBranch), pRemain(pSpouse);
            // Which one to remove completely?
            if (pSpouse->IsRemovableRecombinationLeg(fcsites))
            {
                // They are both removable based on targetable sites.
                if (pBranch->m_rangePtr->GetLiveSites_Rg().empty() && pSpouse->m_rangePtr->GetLiveSites_Rg().empty())
                {
                    // ... and we can't break the tie on live sites, so we randomize.
                    bool removespouse = registry.GetRandom().Bool();
                    if (removespouse)
                    {
                        pRemove = pSpouse;
                        pRemain = pBranch;
                    }
                }
                else
                {
                    // ... and we break the tie on live sites.
                    if (pSpouse->m_rangePtr->GetLiveSites_Rg().empty())
                    {
                        pRemove = pSpouse;
                        pRemain = pBranch;
                    }
                }
            }

            // Unhook branch which will be completely removed.
            Branch_ptr pChild = pRemove->Child(0);
            if (pChild->Parent(0) == pRemove)
            {
                pChild->SetParent(0, pChild->Parent(1));
            }
            pChild->SetParent(1, Branch::NONBRANCH);

            // Remove descending material below it.
            Break(pRemove);

            // Get an interator to the next still-valid interval.
            brit = m_timeList.NextBody(brit);
            if (*brit == pSpouse) brit = m_timeList.NextBody(brit);

            // Splice out branch which remains in tree.
            Branch_ptr pParent0 = pRemain->Parent(0);
            pParent0->ReplaceChild(pRemain, pRemain->Child(0));
            m_timeList.SetUpdateDLs(pParent0);
            Branch_ptr pParent1 = pRemain->Parent(1);
            if (pParent1)
            {
                pParent1->ReplaceChild(pRemain, pRemain->Child(0));
                m_timeList.SetUpdateDLs(pParent1);
            }

            // Remove the dead stuff.
            m_timeList.Remove(pRemove);
            m_timeList.Remove(pRemain);
        }
        else
        {
#if LAMARC_QA_FC_ON
            pBranch->UpdateBranchRange(fcsites, true);
#else
            pBranch->UpdateBranchRange(fcsites, false);
#endif
            brit = m_timeList.NextBody(brit);
        }
    }

    Tree::Prune();
    // MDEBUG validating the ranges
    assert(m_timeList.RevalidateAllRanges());

} // RecTree::Prune

//------------------------------------------------------------------------------------

void RecTree::ReassignDLsFor(string lname, long marker, long ind)
{
    // Find out which locus we're dealing with.  We have its name.  If the locus in question
    // is in m_pLocusVec, we call through to Tree::ReassignDLsFor.
    long locus = FLAGLONG;
    for (unsigned long lnum = 0; lnum < m_pMovingLocusVec->size(); lnum++)
    {
        if ((*m_pMovingLocusVec)[lnum].GetName() == lname)
        {
            locus = lnum;
        }
    }
    if (locus == FLAGLONG)
    {
        // The locus must be in m_pLocusVec--call the base function.
        return Tree::ReassignDLsFor(lname, marker, ind);
    }

    vector<Branch_ptr> haps = m_individuals[ind].GetAllTips();
    vector<LocusCell> cells = m_individuals[ind].GetLocusCellsFor(lname, marker);

    for (unsigned long tip = 0; tip < haps.size(); tip++)
    {
        Cell_ptr origcell = haps[tip]->GetDLCell(locus, markerCell, true);
        Cell_ptr newcell  = cells[tip][0];
        origcell->SetSiteDLs(marker, newcell->GetSiteDLs(marker));
        MarkForDLRecalc(haps[tip]);
    }

    // Don't recalculate aliases, since there are no aliases for the moving locus.
} // RecTree::ReassignDLsFor

//------------------------------------------------------------------------------------

Branch_ptr RecTree::CoalesceActive(double eventT, Branch_ptr active1, Branch_ptr active2, const rangeset & fcsites)
{
    Branch_ptr pBranch = Tree::CoalesceActive(eventT, active1, active2, fcsites);

    pBranch->SetMovingDLCells(CollectMovingCells());

    m_numTargetLinks_Tr -= active1->m_rangePtr->NumTargetLinks_Rg();
    m_numTargetLinks_Tr -= active2->m_rangePtr->NumTargetLinks_Rg();
    m_numTargetLinks_Tr += pBranch->m_rangePtr->NumTargetLinks_Rg();

    return pBranch;
}

//------------------------------------------------------------------------------------

Branch_ptr RecTree::CoalesceInactive(double eventT, Branch_ptr active, Branch_ptr inactive, const rangeset & fcsites)
{
    Branch_ptr pBranch = Tree::CoalesceInactive(eventT, active, inactive, fcsites);

    pBranch->SetMovingDLCells(CollectMovingCells());

    m_numNewTargetLinks_Tr -= inactive->m_rangePtr->NumNewTargetLinks_Rg();
    m_numTargetLinks_Tr    -= active->m_rangePtr->NumTargetLinks_Rg();
    m_numNewTargetLinks_Tr += pBranch->m_rangePtr->NumNewTargetLinks_Rg();
    return pBranch;
}

//------------------------------------------------------------------------------------

Branch_ptr RecTree::Migrate(double eventT, long topop, long maxEvents, Branch_ptr active)
{
    Branch_ptr pBranch = Tree::Migrate(eventT, topop, maxEvents, active);

    //    pBranch->m_rangePtr->SetMRange(pBranch->Child(0)->m_rangePtr);  // PRUNE?
    // A migration does not change the active or newly active links.
    return pBranch;
}

//------------------------------------------------------------------------------------

Branch_ptr RecTree::DiseaseMutate(double eventT, long endstatus, long maxEvents, Branch_ptr active)
{
    Branch_ptr pBranch = Tree::DiseaseMutate(eventT, endstatus, maxEvents, active);

    //    pBranch->m_rangePtr->SetMRange(pBranch->Child(0)->m_rangePtr);  // PRUNE?
    // A disease mutation does not change the active or newly active links.
    return pBranch;
}

//------------------------------------------------------------------------------------

branchpair RecTree::RecombineActive(double eventT, long maxEvents, FPartMap fparts,
                                    Branch_ptr pActive, long rec, const rangeset & fcsites)
{
    if (NumberOfRecombinations() >= maxEvents)
    {
        rec_overrun e;
        throw e;
    }

    // Create the left parent of the active branch, saving partition handling until end.

    bool newbranchisinactive(false);
    rangeset transmittedsites;
    rangepair transmit(0L, rec + 1);
    transmittedsites = AddPairToRange(transmit, transmittedsites);
    RBranch_ptr pParentA = RBranch_ptr(new RBranch(pActive->m_rangePtr, newbranchisinactive, transmittedsites, fcsites));
    pParentA->m_eventTime  = eventT;    // Set the event time.
    pParentA->SetChild(0, pActive);
    pActive->SetParent(0, pParentA);

    // Create the right parent of the active branch, saving partition handling until end.
#if 0 // JDEBUG, apparently can't do this in stl
    transmittedsites.begin()->first = rec + 1;
    transmittedsites.begin()->second = m_totalSites;
#endif
    // Can we just change transmit?
    rangeset transmittedsites2;
    rangepair transmit2(rec + 1, m_totalSites);
    transmittedsites2 = AddPairToRange(transmit2, transmittedsites2);
    RBranch_ptr pParentB = RBranch_ptr(new RBranch(pActive->m_rangePtr, newbranchisinactive, transmittedsites2, fcsites));
    pParentB->m_eventTime  = eventT;    // Set the event time.
    pParentB->SetChild(0, pActive);
    pActive->SetParent(1, pParentB);

    // Now deal with partitions.
    if (fparts.empty())
    {
        pParentA->CopyPartitionsFrom(pActive);
        pParentB->CopyPartitionsFrom(pActive);
    }
    else
    {
        pParentA->RecCopyPartitionsFrom(pActive, fparts, true);
        pParentB->RecCopyPartitionsFrom(pActive, fparts, false);
    }

    m_numTargetLinks_Tr -= pActive->m_rangePtr->NumTargetLinks_Rg();
    m_numTargetLinks_Tr += pParentA->m_rangePtr->NumTargetLinks_Rg();
    m_numTargetLinks_Tr += pParentB->m_rangePtr->NumTargetLinks_Rg();

    m_timeList.Collate(pParentA, pParentB);
    return make_pair(pParentA, pParentB);
}

//------------------------------------------------------------------------------------

branchpair RecTree::RecombineInactive(double eventT, long maxEvents, FPartMap fparts,
                                      Branch_ptr pBranch, long rec, const rangeset & fcsites)
{
    if (NumberOfRecombinations() >= maxEvents)
    {
        rec_overrun e;
        throw e;
    }

    // This is the branch already in the tree.
    bool inactive_is_low = pBranch->m_rangePtr->AreLowSitesOnInactiveBranch_Rg(rec);

    // All this might eventually be one constructor call.
    bool newbranchisinactive(true);
    rangepair transmit;
    if (inactive_is_low)
    {
        transmit.first = 0L;
        transmit.second = rec + 1;
    }
    else
    {
        transmit.first = rec + 1;
        transmit.second = m_totalSites;
    }
    rangeset transmittedsitesInac;
    transmittedsitesInac = AddPairToRange(transmit, transmittedsitesInac);
    RBranch_ptr pParentInac = RBranch_ptr(new RBranch(pBranch->m_rangePtr, newbranchisinactive, transmittedsitesInac, fcsites));
    pParentInac->m_eventTime  = eventT;
    pParentInac->CopyPartitionsFrom(pBranch);

    pParentInac->SetChild(0, pBranch);
    pBranch->Parent(0)->ReplaceChild(pBranch, pParentInac);
    pBranch->SetParent(0, pParentInac);
    if (pBranch->Parent(1))
    {
        pBranch->Parent(1)->ReplaceChild(pBranch, pParentInac);
        pBranch->SetParent(1, Branch::NONBRANCH);
    }

    // This is the new, active branch.
    newbranchisinactive = false;
    if (!inactive_is_low)
    {
        transmit.first = 0L;
        transmit.second = rec + 1;
    }
    else
    {
        transmit.first = rec + 1;
        transmit.second = m_totalSites;
    }
    rangeset transmittedsitesAc;
    transmittedsitesAc = AddPairToRange(transmit, transmittedsitesAc);
    assert(Intersection(transmittedsitesAc, pBranch->m_rangePtr->GetOldTargetSites_Rg()).empty());
    RBranch_ptr pParentAc = RBranch_ptr(new RBranch(pBranch->m_rangePtr, newbranchisinactive, transmittedsitesAc, fcsites));
    pParentAc->m_eventTime  = eventT;

    pParentAc->SetChild(0, pBranch);
    pBranch->SetParent(1, pParentAc);

    if (fparts.empty())
    {
        pParentAc->CopyPartitionsFrom(pBranch);
    }
    else
    {
        pParentAc->RecCopyPartitionsFrom(pBranch, fparts, !inactive_is_low);
    }

    // Now put both branches in place
    m_timeList.Collate(pParentInac, pParentAc);

    m_numNewTargetLinks_Tr -= pBranch->m_rangePtr->NumNewTargetLinks_Rg();
    m_numNewTargetLinks_Tr += pParentInac->m_rangePtr->NumNewTargetLinks_Rg();
    m_numTargetLinks_Tr    += pParentAc->m_rangePtr->NumTargetLinks_Rg();

    return make_pair(pParentInac, pParentAc);
}

//------------------------------------------------------------------------------------

rangevector RecTree::GetLocusSubtrees(rangepair span) const
{
    // Get recombination breakpoints
    set<long> breakpt(GetBreakpoints());

    // Push rangepairs representing each breakpoint into subtree vector.  The (pt + 1) avoids creating a rangepair
    // starting at the last site.  The strange iterator maneuver is because (pt + 1) is not allowed in g++.  We do
    // not risk an empty set; we know it contains at least 2 elements.

    rangevector subtrees;
    set<long>::const_iterator pt = breakpt.begin();
    set<long>::const_iterator nextpt = pt;
    ++nextpt;
    long begin = span.first;
    long end = span.second;

    for ( ; nextpt != breakpt.end(); ++pt, ++nextpt)
    {
        // No overlap so we don't add it.
        if (*nextpt <= begin || *pt >= end) continue;
        // Some overlap, but ends may need adjustment.
        long first = *pt;
        long last = *nextpt;
        if (first < begin) first = begin;
        if (last > end) last = end;
        subtrees.push_back(rangepair(first, last));
    }
    return subtrees;
} // GetLocusSubtrees

//------------------------------------------------------------------------------------

void RecTree::SetMovingLocusVec(vector<Locus> * loc)
{
    m_pMovingLocusVec = loc;
    // We don't have to worry about the range--the moving loci move around within the range of the fixed loci.
} // SetMovingLocusVec

//------------------------------------------------------------------------------------

void RecTree::SetMovingMapposition(long mloc, long site)
{
    assert(mloc < static_cast<long>(m_pMovingLocusVec->size()));
    (*m_pMovingLocusVec)[mloc].SetRegionalMapposition(site);
}

//------------------------------------------------------------------------------------

TBranch_ptr RecTree::CreateTip(const TipData & tipdata, const vector<LocusCell> & cells,
                               const vector<LocusCell> & movingcells, const rangeset & diseasesites)
{
    TBranch_ptr pTip = m_timeList.CreateTip(tipdata, cells, movingcells, m_totalSites, diseasesites);
    return pTip;
}

//------------------------------------------------------------------------------------

TBranch_ptr RecTree::CreateTip(const TipData & tipdata, const vector<LocusCell> & cells,
                               const vector<LocusCell> & movingcells, const rangeset & diseasesites,
                               const vector<Locus> & loci)
{
    TBranch_ptr pTip = m_timeList.CreateTip(tipdata, cells, movingcells, m_totalSites, diseasesites, loci);
    return pTip;
}

//------------------------------------------------------------------------------------

bool RecTree::DoesThisLocusJump(long mloc) const
{
    assert(mloc < static_cast<long>(m_pMovingLocusVec->size()));
    return ((*m_pMovingLocusVec)[mloc].GetAnalysisType() == mloc_mapjump);
}

//------------------------------------------------------------------------------------

bool RecTree::AnyRelativeHaplotypes() const
{
    for (IndVec::const_iterator ind = m_individuals.begin(); ind != m_individuals.end(); ind++)
    {
        if (ind->MultipleTraitHaplotypes()) return true;
    }
    return false;
}

//------------------------------------------------------------------------------------

void RecTree::CalculateDataLikes()
{
    m_overallDL = CalculateDataLikesForFixedLoci();
    //The base function accumulates a likelihood into m_overallDL.
    for (unsigned long loc = 0; loc < m_pMovingLocusVec->size(); ++loc)
    {
        if ((*m_pMovingLocusVec)[loc].GetAnalysisType() == mloc_mapjump)
        {
            m_overallDL += CalculateDataLikesForMovingLocus(loc);
        }
    }
    m_timeList.ClearUpdateDLs();        // Reset the updating flags.
}

//------------------------------------------------------------------------------------

double RecTree::CalculateDataLikesForMovingLocus(long loc)
{
    double likelihood;

    const Locus & locus = (*m_pMovingLocusVec)[loc];
    try                                 // Check for need to switch to normalization.
    {
        likelihood = locus.CalculateDataLikelihood(*this, true);
    }

    catch (datalikenorm_error & ex)     // Normalization is set by thrower.
    {
        m_timeList.SetAllUpdateDLs();   // All subsequent loci will recalculate the entire tree.
        RunReport & runreport = registry.GetRunReport();
        runreport.ReportChat("\n", 0);
        runreport.ReportChat("Subtree of likelihood 0.0 found:  Turning on normalization and re-calculating.");

        likelihood = locus.CalculateDataLikelihood(*this, true);
    }

    return likelihood;
} // CalculateDataLikesForMovingLocus

//------------------------------------------------------------------------------------

set<long> RecTree::GetBreakpoints() const
{
    // Create a (sorted) set containing all recombination breakpoints.
    set<long> breakpt;
    breakpt.insert(0);
    breakpt.insert(m_totalSites);

    Branchconstiter brit;
    long site;
    for (brit = m_timeList.FirstBody(); brit != m_timeList.EndBranch(); brit = m_timeList.NextBody(brit))
    {
        site = (*brit)->GetRecSite_Br();
        if (site != FLAGLONG)           // FLAGLONG means there is no recombination here
            breakpt.insert(site + 1);   // inserting site + 1 because of open interval conventions.
    }

    return breakpt;
} // GetBreakpoints

//------------------------------------------------------------------------------------

DoubleVec2d RecTree::GetMapSummary()
{
    DoubleVec2d retvec;
    DoubleVec1d zeroes(m_totalSites, 0.0);
    for (unsigned long mloc = 0; mloc < m_pMovingLocusVec->size(); mloc++)
    {
        DoubleVec1d mlikes = zeroes;
        long currentsite;
        switch((*m_pMovingLocusVec)[mloc].GetAnalysisType())
        {
            case mloc_mapjump:
                //Find out where it is, make the likelihood at that value 1.
                currentsite = (*m_pMovingLocusVec)[mloc].GetRegionalMapposition();
                mlikes[currentsite] = 1.0;
                break;
            case mloc_mapfloat:
                //Calculate the likelihoods over all subtrees, return vectors with it.
                //mlikes = CalculateDataLikesWithRandomHaplotypesForFloatingLocus(mloc);
                CalculateDataLikesForAllHaplotypesForFloatingLocus(mloc, mlikes);
                break;
            case mloc_data:
            case mloc_partition:
                assert(false);
                throw implementation_error("We seem to want to collect mapping data for a segment without that type"
                                           " of analysis.  This is our fault; e-mail us at lamarc@gs.washington.edu");
        }
        retvec.push_back(mlikes);
    }

    return retvec;
}

//------------------------------------------------------------------------------------

DoubleVec1d RecTree::CalculateDataLikesWithRandomHaplotypesForFloatingLocus(long mloc)
{
    RandomizeMovingHaplotypes(mloc);
    return CalculateDataLikesForFloatingLocus(mloc);
}

//------------------------------------------------------------------------------------

void RecTree::CalculateDataLikesForAllHaplotypesForFloatingLocus(long mloc, DoubleVec1d & mlikes)
{
    long ind = 0;
    //We need to change the vector of mlikes from zeroes to EXPMINS.
    mlikes = SafeLog(mlikes);
    UpdateDataLikesForIndividualsFrom(ind, mloc, mlikes);
    //LS DEBUG
    //PrintTipData(mloc, 0);
    //LS TEST
    //cout << "Sum:  " << endl << ToString(mlikes, 5) << endl;
}

//------------------------------------------------------------------------------------

bool RecTree::UpdateDataLikesForIndividualsFrom(long ind, long mloc, DoubleVec1d & mlikes)
{
    if (static_cast<unsigned long>(ind) == m_individuals.size()) return true;
    string lname = (*m_pMovingLocusVec)[mloc].GetName();
    //LS NOTE: If this ASSERTs, we are mapping something with more than one marker.
    // Might be OK, but should be checked.
    assert ((*m_pMovingLocusVec)[mloc].GetNmarkers() == 1);
    for (long marker = 0; marker < (*m_pMovingLocusVec)[mloc].GetNmarkers(); marker++)
    {
        bool newhaps = true;
        for (m_individuals[ind].ChooseFirstHaplotypeFor(lname, marker);
             newhaps;
             newhaps = m_individuals[ind].ChooseNextHaplotypeFor(lname, marker))
        {
            ReassignDLsFor(lname, marker, ind);
            if (UpdateDataLikesForIndividualsFrom(ind + 1, mloc, mlikes))
            {
                //We're on the last one--update the data likelihood.
                DoubleVec1d newlikes = CalculateDataLikesForFloatingLocus(mloc);
                //LS TEST:
                //cout << ToString(newlikes, 5) << endl;
                mlikes = AddValsOfLogs(mlikes, newlikes);
            }
        }
    }
    return false;
}

//------------------------------------------------------------------------------------

DoubleVec1d RecTree::CalculateDataLikesForFloatingLocus(long mloc)
{
    m_timeList.SetAllUpdateDLs();

    // We always need these on because there's only enough storage for a single subtree,
    // and the last subtree we used was at the other end of the range.
    DoubleVec1d datalikes(m_totalSites);
    rangeset allowedranges = (*m_pMovingLocusVec)[mloc].GetAllowedRange();
    set<long> breakpoints = GetBreakpoints();
    breakpoints = IgnoreDisallowedSubTrees(breakpoints, allowedranges);
    set<long>::iterator left(breakpoints.begin()), right(breakpoints.begin());

    right++;

    for ( ; right != breakpoints.end(); left++, right++ )
    {
        (*m_pMovingLocusVec)[mloc].SetRegionalMapposition(*left);
        double datalike = CalculateDataLikesForMovingLocus(mloc);
        for (long siteindex = *left; siteindex<*right; siteindex++)
        {
            datalikes[siteindex] = datalike;
        }
    }

    datalikes = ZeroDisallowedSites(datalikes, allowedranges);
    m_timeList.ClearUpdateDLs();

    //LS NOTE: We could theoretically not call ClearUpdateDLs if there was only one subtree,
    // but only if the *next* call to this function also had the exact same subtree, which
    // we can't enforce.  So, always call ClearUpdateDLs
    return datalikes;
}

//------------------------------------------------------------------------------------

set<long> RecTree::IgnoreDisallowedSubTrees(set<long> breakpoints, rangeset allowedranges)
{
    set<long> newbreakpoints;
    newbreakpoints.insert(allowedranges.begin()->first);
    newbreakpoints.insert(m_totalSites); // To close out the final range.
    for (rangeset::iterator range = allowedranges.begin(); range != allowedranges.end(); range++)
    {
        for (long site = range->first; site < range->second; site++)
        {
            if (breakpoints.find(site) != breakpoints.end())
            {
                newbreakpoints.insert(site);
            }
        }
    }
    return newbreakpoints;
}

//------------------------------------------------------------------------------------

DoubleVec1d RecTree::ZeroDisallowedSites(DoubleVec1d datalikes, rangeset allowedranges)
{
    DoubleVec1d newdatalikes(datalikes.size(), -DBL_BIG);
    for (rangeset::iterator range = allowedranges.begin(); range != allowedranges.end(); range++)
    {
        assert(range->first >= 0);
        assert(static_cast<unsigned long int>(range->second) <= datalikes.size());
        for (long int site = range->first; site < range->second; site++)
        {
            newdatalikes[site] = datalikes[site];
        }
    }
    return newdatalikes;
}

//------------------------------------------------------------------------------------

void RecTree::RandomizeMovingHaplotypes(long mlocus)
{
    //This function loops over all individuals in m_individuals and randomizes the haplotypes.
    string lname = (*m_pMovingLocusVec)[mlocus].GetName();
    long nmarkers = (*m_pMovingLocusVec)[mlocus].GetNmarkers();
    for (long marker = 0; marker < nmarkers; marker++)
    {
        for (unsigned long ind = 0; ind < m_individuals.size(); ind++)
        {
            if (m_individuals[ind].ChooseRandomHaplotypesFor(lname, marker))
            {
                //The chosen one is different from last time.
                ReassignDLsFor(lname, marker, ind);
            }
        }
    }
}

//------------------------------------------------------------------------------------

bool RecTree::SimulateDataIfNeeded()
{
    if (Tree::SimulateDataIfNeeded())
    {
        //LS DEBUG: Current version uses same data for all moving loci.  If we want to simulate data individually,
        // uncomment the following lines and replace all the subsequent loops over lnums to locus.whatever calls.
        //
        //  for (unsigned long mloc = 0; mloc < m_pMovingLocusVec->size(); mloc++)
        //  {
        //      Locus & locus = (*m_pMovingLocusVec)[mloc];
        if ((*m_pMovingLocusVec).size() > 0)
        {
            Locus & locus = (*m_pMovingLocusVec)[0];
            if (locus.GetShouldSimulate())
            {
                for (size_t lnum = 0; lnum < (*m_pMovingLocusVec).size(); lnum++)
                {
                    (*m_pMovingLocusVec)[lnum].ClearVariability();
                    (*m_pMovingLocusVec)[lnum].SetNewMappositionIfMoving();
                }

                // Check to see if this locus's site was already simulated.
                Locus * simlocus = NULL; //to prevent compiler warnings
                for (unsigned long loc = 0; loc < m_pLocusVec->size(); loc++)
                {
                    Locus & oldlocus = (*m_pLocusVec)[loc];
                    if (oldlocus.GetShouldSimulate())
                    {
                        if (oldlocus.SiteInLocus(locus.GetRegionalMapposition()))
                        {
                            simlocus = & oldlocus;
                        }
                    }
                }
                if (simlocus != NULL)
                {
                    //If there was at least one site in the original that overlapped with the range on the moving
                    // locus, we take that chances as the chance that we hit the original at all, but we want
                    // to choose a variable site from that range.  This introduces an ascertainment bias that I'm
                    // not exactly sure how to compensate for, or if it needs to be compensated for.  But anyway.
                    long newsite = simlocus->ChooseVariableSiteFrom(locus.GetAllowedRange());
                    if (newsite != FLAGLONG)
                    {
                        for (size_t lnum = 0; lnum < (*m_pMovingLocusVec).size(); lnum++)
                        {
                            (*m_pMovingLocusVec)[lnum].SetRegionalMapposition(newsite);
                            (*m_pMovingLocusVec)[lnum].SetTrueSite(newsite);
                            //LS TEST
                            //simlocus->PrintOnesAndZeroesForVariableSites();
                            (*m_pMovingLocusVec)[lnum].SetVariableRange(simlocus->CalculateVariableRange());
                        }
                        //LS DEBUG:  Not really for release.
                        //simlocus->CalculateDisEqFor(newsite);
                        //If there were no variable sites in the original, we keep the originally chosen site.
                        // Our trait isn't going to be variable, but them's the breaks.
                    }
                    for (size_t lnum = 0; lnum < (*m_pMovingLocusVec).size(); lnum++)
                    {
                        (*m_pMovingLocusVec)[lnum].CopyDataFrom(*simlocus, *this);
                    }
                    if (locus.IsNighInvariant())
                    {
                        registry.GetRunReport().ReportNormal("All simulated data was nigh invariant for the"
                                                             " source segment, giving us nigh invariant data"
                                                             " for the simulated segment "
                                                             + locus.GetName() + " as well.");
                    }
                    //LS DEBUG SIM:  This would delete information again, for simulations where
                    // the data exists but you never sequenced it.
                    // simlocus->RandomizeHalf(*this);
                }
                else
                {
                    locus.SimulateData(*this, m_totalSites);
                    // Now copy this data into any other moving loci.
                    long newsite = locus.GetRegionalMapposition();
                    for (size_t lnum = 1; lnum < (*m_pMovingLocusVec).size(); lnum++)
                    {
                        (*m_pMovingLocusVec)[lnum].SetRegionalMapposition(newsite);
                        (*m_pMovingLocusVec)[lnum].SetTrueSite(newsite);
                        (*m_pMovingLocusVec)[lnum].CopyDataFrom(locus, *this);
                    }
                }
                for (size_t lnum = 0; lnum < (*m_pMovingLocusVec).size(); lnum++)
                {
                    (*m_pMovingLocusVec)[lnum].MakePhenotypesFor(m_individuals);
                }
            }
        }
        //Redo the aliases, in case we randomized the data in the original locus.
        SetupAliases(*m_pLocusVec);
        return true;
    }
    return false;
}

//------------------------------------------------------------------------------------

long RecTree::NumberOfRecombinations()
{
    return (m_timeList.HowMany(btypeRec) / 2);
}

//------------------------------------------------------------------------------------

void RecTree::SetTargetLinksFrom(const BranchBuffer & brbuffer)
{
    m_numTargetLinks_Tr = 0;

    vector<Branch_ptr> branches(brbuffer.ExtractConstBranches());
    vector<Branch_ptr>::const_iterator branch;
    for(branch = branches.begin(); branch != branches.end(); ++branch)
    {
        m_numTargetLinks_Tr += (*branch)->m_rangePtr->NumTargetLinks_Rg();
    }
}

//------------------------------------------------------------------------------------
// JDEBUG--this is a Debugging function.

void RecTree::SetNewTargetLinksFrom(const BranchBuffer & brbuffer)
{
    m_numNewTargetLinks_Tr = 0;

    vector<Branch_ptr> branches(brbuffer.ExtractConstBranches());
    vector<Branch_ptr>::const_iterator branch;
    for(branch = branches.begin(); branch != branches.end(); ++branch)
    {
        m_numNewTargetLinks_Tr += (*branch)->m_rangePtr->NumNewTargetLinks_Rg();
    }
}

//------------------------------------------------------------------------------------
// Debugging function.

void RecTree::PrintTipData(long mloc, long marker)
{
    cout << endl << "Tip data:" << endl;
    for (Branchconstiter tip = m_timeList.FirstTip(); tip != m_timeList.FirstBody(); tip++)
    {
        cout << " " << (*m_pMovingLocusVec)[mloc].GetDataModel()->CellToData((*tip)->GetDLCell(mloc, marker, true), marker);
    }
    cout << endl;
}

//------------------------------------------------------------------------------------
// Debugging function.

void RecTree::PrintRangeSetCount(const rangesetcount & rsc)
{
    for (RSCcIter rsci = rsc.begin(); rsci != rsc.end(); rsci++)
    {
        cout << (*rsci).first << ": " << ToString((*rsci).second) << endl;
    }
}

//------------------------------------------------------------------------------------
// Debugging function.

rangesetcount RecTree::RemoveEmpty(const rangesetcount & rsc)
{
    rangesetcount result;
    for (RSCcIter rsci = rsc.begin(); rsci != rsc.end(); rsci++)
    {
        if (!(*rsci).second.empty())
        {
            result.insert(make_pair((*rsci).first, (*rsci).second));
        }
    }
    return result;
}

//____________________________________________________________________________________