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

// $Id: tree.cpp,v 1.113 2012/05/22 20:25:32 jyamato 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 <algorithm>                    // for std::max and min
#include <cassert>
#include <fstream>
#include <limits>                       // for numeric_limits<double>.epsilon()

#include "local_build.h"

#include "branchbuffer.h"               // for Tree::Groom
#include "datapack.h"
#include "dlcell.h"
#include "errhandling.h"                // for datalikenorm_error in Tree::CalculateDataLikes
#include "locus.h"
#include "mathx.h"                      // for exp()
#include "registry.h"
#include "runreport.h"
#include "timemanager.h"
#include "tree.h"
#include "treesum.h"

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

//#define TEST // erynes test

using namespace std;

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

Tree::Tree()
    : m_overallDL(FLAGDOUBLE),
      m_randomSource(&(registry.GetRandom())),
      m_totalSites(0),
      m_timeManager(NULL)
{
    // Deliberately blank.
}

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

Tree::~Tree()
{
    // Deliberately blank due to shared_ptr auto cleanup.
    delete m_timeManager;

} // Tree destructor

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

// MCHECK Can we do this now using Lucian's new mechanism?  NB This copy constructor does *not* copy
// m_firstInvalid (branch equivalence problem needs to be solved).

Tree::Tree(const Tree & tree, bool makestump)
    : m_protoCells(tree.m_protoCells),
      m_individuals(tree.m_individuals),
      m_overallDL(tree.m_overallDL),
      m_randomSource(tree.m_randomSource),
      m_pLocusVec(tree.m_pLocusVec),
      m_totalSites(tree.m_totalSites)
{
    const ForceSummary & fs(registry.GetForceSummary());
    m_timeManager = fs.CreateTimeManager();

    if (!makestump)
    {
        *m_timeManager = *(tree.m_timeManager); // Make sure the TimeManager state is the same in both trees.
        m_timeList = tree.m_timeList;   // Deep copy the branchlist.

        // Now setup the tip pointers in the individuals correctly, since the
        // individual copy ctor only makes empty containers of tips.
        unsigned long indiv;
        for(indiv = 0; indiv < m_individuals.size(); ++indiv)
        {
            StringVec1d tipnames = tree.m_individuals[indiv].GetAllTipNames();
            vector<Branch_ptr> itips = GetTips(tipnames);
            m_individuals[indiv].SetTips(itips);
        }
    }

}  // Tree copy ctor

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

void Tree::Clear()
{
    m_timeList.Clear();
}

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

const vector<LocusCell> & Tree::CollectCells()
{
    if (m_protoCells.empty())
    {
        unsigned long i;
        for (i = 0; i < m_pLocusVec->size(); ++i)
        {
            m_protoCells.push_back((*m_pLocusVec)[i].GetProtoCell());
        }
    }
    return m_protoCells;

} // CollectCells

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

void Tree::CopyTips(const Tree * tree)
{
    m_timeList.CopyTips(tree->m_timeList);
    m_firstInvalid = m_timeList.BeginBranch();
    m_totalSites = tree->m_totalSites;

    m_overallDL = FLAGDOUBLE;           // It won't be valid anymore.
    m_individuals = tree->m_individuals;

    // Now setup the tip pointers in the individuals correctly.
    unsigned long indiv;
    for(indiv = 0; indiv < m_individuals.size(); ++indiv)
    {
        StringVec1d tipnames = tree->m_individuals[indiv].GetAllTipNames();
        vector<Branch_ptr> itips = GetTips(tipnames);
        m_individuals[indiv].SetTips(itips);
    }

    vector<Locus>::const_iterator locus;
    m_aliases.clear();

    for(locus = m_pLocusVec->begin(); locus != m_pLocusVec->end(); ++locus)
        m_aliases.push_back(locus->GetDLCalc()->RecalculateAliases(*this, *locus));
}

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

void Tree::CopyBody(const Tree * tree)
{
    m_timeList.CopyBody(tree->m_timeList);
    m_overallDL = tree->m_overallDL;
}

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

void Tree::CopyStick(const Tree * tree)
{
    m_timeManager->CopyStick(*(tree->m_timeManager));
}

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

void Tree::CopyPartialBody(const Tree * tree)
{
    // It's wrong to call CopyPartialBody if you don't have a starting point.
    assert((*m_firstInvalid) != Branch::NONBRANCH);

    // We can't do a partial copy if the tree was changed too high up.
    if ((*m_firstInvalid)->BranchGroup() == bgroupTip || (*tree->m_firstInvalid)->BranchGroup() == bgroupTip)
    {
        m_timeList.CopyBody(tree->m_timeList);
    }
    else
    {
        m_timeList.CopyPartialBody(tree->m_timeList, tree->m_firstInvalid, m_firstInvalid);
    }

    m_overallDL = tree->m_overallDL;
}

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

TBranch_ptr Tree::CreateTip(const TipData & tipdata, const vector<LocusCell> & cells,
                            const vector<LocusCell> & movingcells, const rangeset & diseasesites)
{
    // used in sample only case
    TBranch_ptr pTip = m_timeList.CreateTip(tipdata, cells, movingcells, m_totalSites, diseasesites);
    return pTip;
}

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

TBranch_ptr Tree::CreateTip(const TipData & tipdata, const vector<LocusCell> & cells,
                            const vector<LocusCell> & movingcells, const rangeset & diseasesites,
                            const vector<Locus> & loci)
{
    // used in panel case
    TBranch_ptr pTip = m_timeList.CreateTip(tipdata, cells, movingcells, m_totalSites, diseasesites,
                                            loci);
    return pTip;
}
//------------------------------------------------------------------------------------

void Tree::SetTreeTimeManager(TimeManager * tm)
{
    m_timeManager = tm;
}

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

void Tree::SetIndividualsWithTips(const vector<Individual> & indvec)
{
    m_individuals = indvec;
} // SetIndividualsWithTips

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

void Tree::SetLocusVec(vector<Locus> * loc)
{
    m_pLocusVec = loc;
    m_totalSites = 0;
    vector<Locus>::const_iterator locit;
    for(locit = loc->begin(); locit != loc->end(); ++locit)
    {
        if (locit->GetSiteSpan().second > m_totalSites)
            m_totalSites = locit->GetSiteSpan().second;
    }

} // SetLocusVec

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

TreeSummary * Tree::SummarizeTree() const
{
    TreeSummary * treesum = registry.GetProtoTreeSummary().Clone();
    treesum->Summarize(*this);
    return treesum;

} // SummarizeTree

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

double Tree::CalculateDataLikesForFixedLoci()
{
    unsigned long loc;
#ifndef NDEBUG
    vector<double> likes;
#endif
    double likelihood, overalllike(0);

    for (loc = 0; loc < m_pLocusVec->size(); ++loc)
    {
        const Locus & locus = (*m_pLocusVec)[loc];
        try
        {                               // Check for need to switch to normalization.
            likelihood = locus.CalculateDataLikelihood(*this, false);
        }
        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, false);
        }
#ifndef NDEBUG
        likes.push_back(likelihood);
#endif
        //cout << "short loc: " << loc << " likelihood: " << likelihood << endl;
        overalllike += likelihood;      // Add to accumulator.

#if LIKETRACK
        ofstream of;
        of.open("likes1", ios::app);
        of << GetDLValue() << " ";
        of.close();
#endif
    }

#ifndef NDEBUG
    // This code tests a full likelihood recalculation against the time-saving partial recalculation.  Slow but sure.
    double checkDL = 0.0;
    m_timeList.SetAllUpdateDLs();
    for (loc = 0; loc < m_pLocusVec->size(); ++loc)
    {
        const Locus & locus = (*m_pLocusVec)[loc];
        double likelihood = locus.CalculateDataLikelihood(*this, false);
        //cout << "long loc: " << loc << " likelihood: " << likelihood << endl << endl ;
        checkDL += likelihood;
#ifndef LAMARC_QA_SINGLE_DENOVOS
        // likelihood can be zero in single denovos test
        assert(fabs(likelihood - likes[loc]) / likelihood < EPSILON);
#endif // LAMARC_QA_SINGLE_DENOVOS
    }
#ifndef LAMARC_QA_SINGLE_DENOVOS
    // overalllike can be zero in single denovos test
    assert(fabs(overalllike - checkDL) / overalllike < EPSILON);
#endif // LAMARC_QA_SINGLE_DENOVOS
#endif // NDEBUG

    return overalllike;

} // CalculateDataLikesForFixedLoci

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

void Tree::SetupAliases(const vector<Locus> & loci)
{
    // Set up the aliases; this must be done after filling the tips with data!
    m_aliases.clear();
    vector<Locus>::const_iterator locus;

    for (locus = loci.begin(); locus != loci.end(); ++locus)
    {
        m_aliases.push_back(locus->GetDLCalc()->RecalculateAliases(*this, *locus));
    }

} // SetupAliases

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

vector<Branch_ptr> Tree::ActivateTips(Tree * othertree)
{
    vector<Branch_ptr> tips;
    m_firstInvalid = m_timeList.BeginBranch();
    othertree->SetFirstInvalid();

    m_timeList.ClearBody();

    Branchiter brit;
    for (brit = m_timeList.FirstTip(); brit != m_timeList.EndBranch(); brit = m_timeList.NextTip(brit))
    {
        tips.push_back(*brit);
    }

    return tips;
}

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

Branch_ptr Tree::ActivateBranch(Tree * othertree)
{
    // NCuttable - 1 because the root branch is of a cuttable type but is not, in fact, a cuttable branch.
    // + 1 because it's helpful to count cuttable branches from 1, not 0.
    long rn = m_randomSource->Long(m_timeList.GetNCuttable() - 1) + 1;

    Branchiter brit;
    Branch_ptr pActive = Branch::NONBRANCH;
    assert(m_timeList.IsValidTimeList());

    for (brit = m_timeList.BeginBranch(); rn > 0; ++brit)
    {
        m_firstInvalid = brit;
        pActive = *brit;
        rn -= pActive->Cuttable();
    }

    Branchconstiter constfirstinvalid = m_firstInvalid;
    othertree->SetFirstInvalidFrom(constfirstinvalid);

    assert(pActive != Branch::NONBRANCH); // We didn't find any cuttable branches?

    Break(pActive);
    pActive->SetParent(0, Branch::NONBRANCH);
    pActive->SetParent(1, Branch::NONBRANCH);
    return pActive;
}

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

void Tree::SetFirstInvalidFrom(Branchconstiter & target)
{
    m_firstInvalid = m_timeList.FindEquivBranch(*target);
} // SetFirstInvalidFrom

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

Branch_ptr Tree::ActivateRoot(FC_Status & fcstatus)
{
    Branch_ptr pBase   = m_timeList.Base(); // Get the tree base and root.
    Branch_ptr pRoot   = pBase->Child(0);
    pBase->SetChild(0, Branch::NONBRANCH); // Disconnect them.
    pRoot->SetParent(0, Branch::NONBRANCH);

    // Update the fcstatus with root info.
#if LAMARC_QA_FC_ON
    if (pRoot->Child(1))                // We have two children.
    {
        rangeset sites(Intersection(pRoot->Child(0)->GetLiveSites_Br(), pRoot->Child(1)->GetLiveSites_Br()));
        fcstatus.Decrement_FC_Counts(sites);
    }
#endif

#if LAMARC_QA_FC_ON
    pRoot->UpdateRootBranchRange(fcstatus.Coalesced_Sites(), true); // Update the root's range object.
#else
    rangeset emptyset;
    pRoot->UpdateRootBranchRange(emptyset, false);
#endif
    return pRoot;
}

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

Branch_ptr Tree::ChoosePreferentiallyTowardsRoot(Tree * othertree)
{
    // Choose a random body branch, triangle weighted towards the root.
    long nbranches = m_timeList.GetNBodyBranches();
    long weightedbranches = (nbranches * (nbranches + 1)) / 2;
    // weightedbranches = 0; //JRM remove
    long chosenweight = m_randomSource->Long(weightedbranches) + 1;
    long chosenbranch(0L), weight(0L);
    while(weight < chosenweight)
    {
        ++chosenbranch;
        weight += chosenbranch;
    }

    long thisbranch(1L);
    for(m_firstInvalid = m_timeList.FirstBody();
        m_firstInvalid != m_timeList.EndBranch();
        m_firstInvalid = m_timeList.NextBody(m_firstInvalid), ++thisbranch)
    {
        if (thisbranch == chosenbranch) break;
    }
    assert(m_firstInvalid != m_timeList.EndBranch() && thisbranch == chosenbranch);

    // m_firstInvalid may be set incorrectly when there are multiple branches with the same eventtime.
    m_firstInvalid = m_timeList.ResolveTiedTimes(m_firstInvalid);

    Branchconstiter constfirstinvalid = m_firstInvalid;
    othertree->SetFirstInvalidFrom(constfirstinvalid);

    return *m_firstInvalid;
}

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

Branch_ptr Tree::ChooseFirstBranchInEpoch(double targettime, Tree * othertree)
{
    for(m_firstInvalid = m_timeList.FirstBody();
        m_firstInvalid != m_timeList.EndBranch();
        m_firstInvalid = m_timeList.NextBody(m_firstInvalid))
    {
        if ((*m_firstInvalid)->GetTime() > targettime) break;
    }

    assert(m_firstInvalid != m_timeList.EndBranch());

    // m_firstInvalid may be set incorrectly when there are multiple branches with the same eventtime.
    m_firstInvalid = m_timeList.ResolveTiedTimes(m_firstInvalid);

    Branchconstiter constfirstinvalid = m_firstInvalid;
    othertree->SetFirstInvalidFrom(constfirstinvalid);
    return *m_firstInvalid;
}

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

void Tree::AttachBase(Branch_ptr newroot)
{
    Branch_ptr pBase     = m_timeList.Base();  // Get the tree base.
    pBase->SetChild(0, newroot);               // Link the root to the base.
    newroot->SetParent(0, pBase);
}

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

vector<Branch_ptr> Tree::FindAllBranchesAtTime(double eventT)
{
    Branch_ptr pBranch;
    vector<Branch_ptr> inactives;

    Branchiter brit;
    for (brit = m_timeList.BeginBranch(); brit != m_timeList.EndBranch(); ++brit)
    {
        pBranch = *brit;
        if (eventT < pBranch->m_eventTime)
        {
            break;
        }

        if (pBranch->Parent(0))
        {
            if (eventT < pBranch->Parent(0)->m_eventTime)
            {
                inactives.push_back(pBranch);
            }
        }
        else assert(!pBranch->Parent(1));
    }

    return inactives;
}

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

vector<Branch_ptr> Tree::FirstInterval(double eventT)
{
    return FindAllBranchesAtTime(eventT);
}

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

vector<Branch_ptr> Tree::FindBranchesImmediatelyTipwardOf(Branchiter startbr)
{
    assert((*startbr)->m_eventTime != 0.0);
    Branchiter prevbr(startbr);
    --prevbr;                           // Assumes we have reversible iterator.
    while((*prevbr)->m_eventTime == (*startbr)->m_eventTime)
    {
        assert(prevbr != m_timeList.BeginBranch()); // Too far!
        --prevbr;
    }
    return FindAllBranchesAtTime((*prevbr)->m_eventTime);
}

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

vector<Branch_ptr> Tree::FindBranchesBetween(double startint, double endint)
{
    vector<Branch_ptr> branches;

    Branchiter brit;
    for (brit = m_timeList.BeginBranch(); brit != m_timeList.EndBranch(); ++brit)
    {
        Branch_ptr pBranch = *brit;
        if (startint <= pBranch->m_eventTime && pBranch->m_eventTime < endint)
        {
            branches.push_back(pBranch);
        }
    }

    return branches;
}

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

void Tree::NextInterval(Branch_ptr branch)
{
    // Deliberately blank; needed only in rectree.
}

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

void Tree::Prune()
{
    m_timeList.Prune();
    // This is done last so that the information it overwrites can be used during
    // pruning and for validation--do not move it into the main loop!
    Branchiter brit;
    for (brit = m_timeList.FirstBody(); brit != m_timeList.EndBranch(); brit = m_timeList.NextBody(brit))
    {
        (*brit)->ResetBuffersForNextRearrangement();
    }
    TrimStickToRoot();
}

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

void Tree::TrimStickToRoot()
{
    m_timeManager->ChopOffStickAt(m_timeList.RootTime());
} // TrimStickToRoot

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

Branch_ptr Tree::Coalesce(Branch_ptr child1, Branch_ptr child2, double tevent, const rangeset & fcsites)
{
    assert(child1->m_partitions == child2->m_partitions);

    bool newbranchisinactive(false);
    Branch_ptr parent = Branch_ptr(new CBranch(child1->m_rangePtr, child2->m_rangePtr, newbranchisinactive, fcsites));
    parent->m_eventTime = tevent;
    parent->CopyPartitionsFrom(child1);

    parent->SetChild(0, child1);
    parent->SetChild(1, child2);
    child1->SetParent(0, parent);
    child2->SetParent(0, parent);

    parent->SetDLCells(CollectCells());
    parent->SetUpdateDL();              // Mark this branch for updating.

    m_timeList.Collate(parent);
    assert(child1->m_partitions == parent->m_partitions);

    return parent;

} // Tree::Coalesce

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

Branch_ptr Tree::CoalesceActive(double eventT, Branch_ptr pActive1, Branch_ptr pActive2, const rangeset & fcsites)
{
    // Create the parent of the two branches.
    Branch_ptr pParent = Coalesce(pActive1, pActive2, eventT, fcsites);
    return pParent;
}

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

Branch_ptr Tree::CoalesceInactive(double eventT, Branch_ptr pActive, Branch_ptr pInactive, const rangeset & fcsites)
{
    // Create the parent of the two branches.
    bool newbranchisinactive(true);
    Branch_ptr pParent = Branch_ptr(new CBranch(pInactive->m_rangePtr, pActive->m_rangePtr, newbranchisinactive, fcsites));
    pParent->m_eventTime  = eventT;
    pParent->CopyPartitionsFrom(pActive);

    // Hook children to parent.
    pParent->SetChild(0, pInactive);
    pParent->SetChild(1, pActive);

    // Hook pInactive's old parent up to new parent.
    pInactive->Parent(0)->ReplaceChild(pInactive, pParent);
    pInactive->SetParent(0, pParent);
    if (pInactive->Parent(1))
    {
        pInactive->Parent(1)->ReplaceChild(pInactive, pParent);
        pInactive->SetParent(1, Branch::NONBRANCH);
    }

    pActive->SetParent(0, pParent);

    pParent->SetDLCells(CollectCells());
    m_timeList.SetUpdateDLs(pParent);   // Set the update data likelihood flag.

    m_timeList.Collate(pParent);
    return pParent;
}

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

Branch_ptr Tree::Migrate(double eventT, long pop, long maxEvents, Branch_ptr pActive)
{
    if (m_timeList.HowMany(btypeMig) >= maxEvents)
    {
        mig_overrun e;
        throw e;
    }

    force_type force;
    if (registry.GetForceSummary().CheckForce(force_MIG))
    {
        force = force_MIG;
    }
    else
    {
        force = force_DIVMIG;
    }

    Branch_ptr pParent;
    if (force == force_MIG)
    {
        pParent = Branch_ptr(new MBranch(pActive->m_rangePtr));
    }
    else
    {
        pParent = Branch_ptr(new DivMigBranch(pActive->m_rangePtr));
    }

    pParent->m_eventTime  = eventT;
    pParent->CopyPartitionsFrom(pActive);
    pParent->SetPartition(force, pop);

    pParent->SetChild(0, pActive);
    pActive->SetParent(0, pParent);

    m_timeList.Collate(pParent);
    return pParent;
}

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

Branch_ptr Tree::DiseaseMutate(double eventT, long endstatus, long maxevents, Branch_ptr pActive)
{
    if (m_timeList.HowMany(btypeDisease) >= maxevents)
    {
        dis_overrun e;
        throw e;
    }

    Branch_ptr pParent = Branch_ptr(new DBranch(pActive->m_rangePtr));
    pParent->m_eventTime  = eventT;
    pParent->CopyPartitionsFrom(pActive);
    pParent->SetPartition(force_DISEASE, endstatus);

    pParent->SetChild(0, pActive);
    pActive->SetParent(0, pParent);

    m_timeList.Collate(pParent);

    return pParent;

} // DiseaseMutate

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

Branch_ptr Tree::TransitionEpoch(double eventT, long newpop, long maxevents, Branch_ptr pActive)
{
    if (m_timeList.HowMany(btypeEpoch) >= maxevents)
    {
        epoch_overrun e;
        throw e;
    }

    Branch_ptr pParent = Branch_ptr(new EBranch(pActive->m_rangePtr));
    pParent->m_eventTime = eventT;
    pParent->CopyPartitionsFrom(pActive);
    pParent->SetPartition(force_DIVMIG, newpop);

    pParent->SetChild(0, pActive);
    pActive->SetParent(0, pParent);

    m_timeList.Collate(pParent);
    return pParent;

} // TransitionEpoch

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

vector<Branch_ptr> Tree::FindBranchesStartingOnOpenInterval(
    double starttime, double endtime)
{
    vector<Branch_ptr> branches;
    Branchiter br;
    for(br = m_timeList.BeginBranch(); br != m_timeList.EndBranch(); ++br)
    {
        if ((*br)->m_eventTime <= starttime) continue;
        if ((*br)->m_eventTime >= endtime) break;
        branches.push_back(*br);
    }

    return branches;

} // FindBranchesStartingOnOpenInterval

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

vector<Branch_ptr> Tree::FindEpochBranchesAt(double time)
{
    vector<Branch_ptr> branches;
    Branchiter br;
    for(br = m_timeList.BeginBranch(); br != m_timeList.EndBranch(); ++br)
    {
        if ((*br)->m_eventTime == time && (*br)->Event() == btypeEpoch)
            branches.push_back(*br);
    }

    return branches;

} // FindEpochBranchesAt

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

vector<Branch_ptr> Tree::FindBranchesStartingRootwardOf(double time)
{
    vector<Branch_ptr> branches;
    Branchiter br;
    for(br = m_timeList.BeginBranch(); br != m_timeList.EndBranch(); ++br)
    {
        if ((*br)->m_eventTime > time)
            branches.push_back(*br);
    }

    return branches;

} // FindBranchesStartingRootwardOf

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

void Tree::SwapSiteDLs()
{
    Individual ind;

    // Pick an individual with phase-unknown sites at random.
    do {
        ind  = m_individuals[m_randomSource->Long(m_individuals.size())];
    } while (!ind.AnyPhaseUnknownSites());

    // Pick a phase marker at random.
    pair<long, long> locusmarker = ind.PickRandomPhaseMarker(*m_randomSource);
    long locus = locusmarker.first;
    long posn = locusmarker.second;

    // Pick two tips in the individual at random.
    vector<Branch_ptr> haps = ind.GetAllTips();
    long size = haps.size();
    long rnd1 = m_randomSource->Long(size);
    long rnd2 = m_randomSource->Long(size - 1);
    if (rnd1 <= rnd2)
        rnd2++;

    Branch_ptr pTip1   = haps[rnd1];
    Branch_ptr pTip2   = haps[rnd2];

    // NB: Swap first DLCell only!

    assert(pTip1->GetNcells(0) == pTip2->GetNcells(0));

    // Get the data likelihood arrays.
    Cell_ptr dlcell1 = pTip1->GetDLCell(locus, markerCell, false);
    Cell_ptr dlcell2 = pTip2->GetDLCell(locus, markerCell, false);

    // Swap the cells at the marker position.
    dlcell1->SwapDLs(dlcell2, posn);

    // Take care of needed data-likelihood corrections.
    const Locus & loc = (*m_pLocusVec)[locus];
    m_aliases[locus] = loc.GetDLCalc()->RecalculateAliases(*this, loc);

    MarkForDLRecalc(pTip1);
    MarkForDLRecalc(pTip2);

} // Tree::SwapSiteDLs

//------------------------------------------------------------------------------------
// PickNewSiteDLs is used when haplotype arranging for non-50/50 haplotype probabilities.
// Otherwise, SwapSiteDLs is used.

void Tree::PickNewSiteDLs()
{
    long ind;

    // Pick an individual with phase-unknown sites at random.
    do {
        ind  = m_randomSource->Long(m_individuals.size());
    } while (!m_individuals[ind].MultipleTraitHaplotypes());

    // Pick a locus/marker at random.
    pair<string, long> locusmarker = m_individuals[ind].PickRandomHaplotypeMarker();
    string lname = locusmarker.first;
    long marker  = locusmarker.second;

    // Have the individual in question pick up a new set of haplotypes.
    m_individuals[ind].ChooseNewHaplotypesFor(lname, marker);

    // Finish this routine in a subroutine, because chances are good that the locus in question
    // is a member of m_pMovingLocusVec, which we only have if we're a RecTree.

    ReassignDLsFor(lname, marker, ind);
    return;

} // PickNewSiteDLs

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

void Tree::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_pMovingLocusVec, we've already dealt with it in RecTree::ReassignDLsFor.
    long locus = FLAGLONG;
    for (unsigned long lnum = 0; lnum < m_pLocusVec->size(); lnum++)
    {
        if ((*m_pLocusVec)[lnum].GetName() == lname)
        {
            locus = lnum;
        }
    }
    if (locus == FLAGLONG)
    {
        throw implementation_error("We have the segment name " + lname
                                   + ", but no segment matches it.  Did something happen"
                                   " to this tree?  To the segment in question?");
    }

    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, false);
        Cell_ptr newcell  = cells[tip][0];
        origcell->SetSiteDLs(marker, newcell->GetSiteDLs(marker));
        MarkForDLRecalc(haps[tip]);
    }

    // take care of needed data-likelihood corrections
    const Locus & loc = (*m_pLocusVec)[locus];
    m_aliases[locus] = loc.GetDLCalc()->RecalculateAliases(*this, loc);

} // Tree::ReassignDLsFor

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

void Tree::MarkForDLRecalc(Branch_ptr markbr)
{
    markbr->SetUpdateDL();
    markbr->MarkParentsForDLCalc();

} // Tree::MarkForDLRecalc

//------------------------------------------------------------------------------------
// This version is used to limit a multi-locus case to just the locus under consideration.

rangevector Tree::GetLocusSubtrees(rangepair span) const
{
    rangevector subtrees;
    subtrees.push_back(span);
    return subtrees;

} // GetLocusSubtrees

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

vector<Branch_ptr> Tree::GetTips(StringVec1d & tipnames) const
{
    vector<Branch_ptr> tips;
    StringVec1d::iterator tname;
    for(tname = tipnames.begin(); tname != tipnames.end(); ++tname)
        tips.push_back(GetTip(*tname));

    return tips;

} // GetTips

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

Branch_ptr Tree::GetTip(const string & name) const
{
    return m_timeList.GetTip(name);
} // GetTip

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

bool Tree::NoPhaseUnknownSites() const
{
    IndVec::const_iterator ind;
    unsigned long locus;
    for (locus = 0; locus < m_pLocusVec->size(); ++locus)
    {
        for(ind = m_individuals.begin(); ind != m_individuals.end(); ++ind)
            if (ind->AnyPhaseUnknownSites()) return false;
    }

    return true;

} // Tree:NoPhaseUnknownSites

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

long Tree::GetNsites() const
{
long nsites(0L);
vector<Locus>::iterator locus;
for(locus = m_pLocusVec->begin() ; locus != m_pLocusVec->end(); ++locus)
{
    nsites += locus->GetNsites();
}

return nsites;

} // Tree::GetNsites

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

bool Tree::IsValidTree() const
{
    // NB:  This can only be called on a completed tree, not one in mid-rearrangement.
    // Each Branch validates itself.
    Branchconstiter brit = m_timeList.BeginBranch();
    long ncuttable = 0;

    for ( ; brit != m_timeList.EndBranch(); ++brit)
    {
        //LS TEST
        //(*brit)->PrintInfo_Br();
        if (!(*brit)->CheckInvariant())
        {
            return false;
        }
        ncuttable += (*brit)->Cuttable();
    }

    if (ncuttable != m_timeList.GetNCuttable())
    {
        return false;
    }
    return m_timeList.IsValidTimeList();

    //  if (FindContradictParts()) return false;
    //  return true;

} // IsValidTree

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

bool Tree::ConsistentWithParameters(const ForceParameters& fp) const
{
  Branchconstiter brit;
  vector<double> epochtimes = fp.GetEpochTimes();
  for (brit=m_timeList.FirstBody(); brit != m_timeList.EndBranch();
    brit = m_timeList.NextBody(brit))
  {
    if ((*brit)->Event() == btypeEpoch) {
      if (find(epochtimes.begin(), epochtimes.end(),
        (*brit)->m_eventTime) == epochtimes.end()) {
         return false; // tree inconsistent with parameters
      }
    }
  }
  return true;
}

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

bool Tree::operator==(const Tree & src) const
{
    if (m_pLocusVec->size() != src.m_pLocusVec->size()) return false;
    unsigned long i;
    for (i = 0; i < m_pLocusVec->size(); ++i)
    {
        if ((*m_pLocusVec)[i].GetNsites() != (*src.m_pLocusVec)[i].GetNsites()) return false;
    }
    if (m_overallDL != src.m_overallDL) return false;
    if (m_timeList != src.m_timeList) return false;

    // WARNING I don't guarantee these are the same tree, especially if
    // we are in mid-rearrangement, but they are at least very similar. --Mary
    return true;

} // operator==

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

bool Tree::SimulateDataIfNeeded()
{
    bool simulated = false;
    for (unsigned long loc = 0; loc < m_pLocusVec->size(); loc++)
    {
        Locus & locus = (*m_pLocusVec)[loc];
        if (locus.GetShouldSimulate())
        {
            simulated = true;
            locus.SimulateData(*this, m_totalSites);
        }
    }

    // Redo the aliases.
    SetupAliases(*m_pLocusVec);

    // Note that we have to calculate all DLCells.
    m_timeList.SetAllUpdateDLs();

    return simulated;
}

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

void Tree::DLCheck(const Tree & other) const
{
    cout << m_timeList.DLCheck(other.GetTimeList()) << endl;
} // DLCheck

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

void Tree::PrintStickThetasToFile(ofstream & of) const
{
    m_timeManager->PrintStickThetasToFile(of);
} // PrintStickThetasToFile

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

void Tree::PrintStickFreqsToFile(ofstream & of) const
{
    m_timeManager->PrintStickFreqsToFile(of);
} // PrintStickFreqsToFile

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

void Tree::PrintStickFreqsToFileAtTime(ofstream & of, double time) const
{
    m_timeManager->PrintStickFreqsToFileAtTime(of, time);
} // PrintStickFreqsToFileAtTime

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

void Tree::PrintStickThetasToFileForJoint300(ofstream & of) const
{
    m_timeManager->PrintStickThetasToFileForJoint300(of);
} // PrintStickThetasToFileForJoint300

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

void Tree::PrintStickToFile(ofstream & of) const
{
    m_timeManager->PrintStickToFile(of);
} // PrintStickAndSummaryToFile

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

void Tree::PrintDirectionalMutationEventCountsToFile(ofstream & of) const
{
    m_timeList.PrintDirectionalMutationEventCountsToFile(of);
} // PrintDirectionalMutationEventCountsToFile

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

void Tree::PrintTimeTilFirstEventToFile(ofstream & of) const
{
    m_timeList.PrintTimeTilFirstEventToFile(of);
} // PrintTimeTilFirstEventToFile

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

void Tree::PrintTraitPhenotypeAtLastCoalescence(ofstream & of) const
{
    m_timeList.PrintTraitPhenotypeAtLastCoalescence(of);
} // PrintTraitPhenotypeAtLastCoalescence

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

void Tree::DestroyStick()
{
    m_timeManager->ClearStick();
} // DestroyStick

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

void Tree::SetStickParams(const ForceParameters & fp)
{
    m_timeManager->SetStickParameters(fp);
} // SetStickParams

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

bool Tree::UsingStick() const
{
    return m_timeManager->UsingStick();
} // UsingStick

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

void Tree::ScoreStick(TreeSummary & treesum) const
{
    m_timeManager->ScoreStick(treesum);
} // ScoreStick

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

DoubleVec1d Tree::XpartThetasAtT(double time, const ForceParameters & fp) const
{
    return m_timeManager->XpartThetasAtT(time, fp);
} // XpartThetasAtT

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

DoubleVec1d Tree::PartitionThetasAtT(double time, const force_type force, const ForceParameters & fp) const
{
    return m_timeManager->PartitionThetasAtT(time, force, fp);
} // XpartThetasAtT

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

void Tree::SetNewTimesFrom(Branchiter startpoint, const DoubleVec1d & newtimes)
{
    assert(IsValidTree());
    assert(!newtimes.empty());

    bool updatefirstinvalid(startpoint == m_firstInvalid);

    // First retime and remove all the branches from startpoint on...
    vector<Branch_ptr> newbranches;
    Branchiter branch(startpoint);
    DoubleVec1d::const_iterator newtime(newtimes.begin());

    for( ; newtime != newtimes.end(); ++newtime)
    {
        assert(branch != m_timeList.EndBranch());
        // No special handling for recombinations needed; assumed handled in calling code that sets up "newtimes".
        (*branch)->m_eventTime = *newtime;
        Branch_ptr newbr(*branch);
        branch = m_timeList.NextBody(branch);
        m_timeList.Remove(newbr);
        newbranches.push_back(newbr);
    }

    // Now put them back in, setting datalikelihood update flags as needed.
    vector<Branch_ptr>::iterator newbranch;
    for(newbranch = newbranches.begin(); newbranch != newbranches.end(); ++newbranch)
    {
        Branchiter brit(m_timeList.Collate(*newbranch));
        if (newbranch == newbranches.begin() && updatefirstinvalid)
        {
            m_firstInvalid = brit;
        }
        m_timeList.SetUpdateDLs(*newbranch);
    }

    assert(IsValidTree());

} // SetNewTimesFrom

//------------------------------------------------------------------------------------
// Method to, if necessary, shrink interval lengths, usually at the root, in an attempt to prevent explosions in the
// next maximization.  See the extensive comment in ChainManager::DoChain(), which is where this method is called.
//
// TEMPERATURE arg used only for debugging.

bool Tree::GroomForGrowth(const DoubleVec1d & thetas, const DoubleVec1d & growths, double temperature)
{
    if (thetas.empty() || growths.empty() || thetas.size() != growths.size())
        throw implementation_error("Tree::GroomForGrowth() received an invalid theta and/or growth vector");

    bool treeWasModified = false;
    double starttime(0.0), cumulativeShift(0.0), prevIntervalLength(DBL_BIG);
    Branchiter brit;
    BranchBuffer lineages(registry.GetDataPack());
    for (brit = m_timeList.FirstTip(); brit != m_timeList.EndBranch(); brit = m_timeList.NextTip(brit))
        lineages.UpdateBranchCounts((*brit)->m_partitions, true);

    // Iterate through the timelist, looking for fatal intervals.
    for (brit = m_timeList.FirstBody(); brit != m_timeList.EndBranch(); brit = m_timeList.NextBody(brit))
    {
        Branch_ptr pBranch = *brit;
        unsigned long nchildren(0UL), worst_pop(0UL), npops = thetas.size();
        double smallestExpectedLength(DBL_BIG); // Smallest E[dt] among the pops.
        double expectedLength;
        LongVec1d xpartitions = lineages.GetBranchXParts();
        string msg;

        // Once we shrink an interval, we need to shift all intervals below it by this amount.
        // The amount accumulates.
        if (cumulativeShift > 0.0)
        {
            pBranch->m_eventTime -= cumulativeShift;
            pBranch->SetUpdateDL();
        }

        // Determine how many child branches we will need to update.
        switch(pBranch->Event())
        {
            case btypeCoal:
                nchildren = static_cast<unsigned long>(NELEM);
                break;
            case btypeMig:
            case btypeDivMig:
            case btypeDisease:
            case btypeEpoch:
                nchildren = 1UL;
                break;
            case btypeRec:
                nchildren = 0UL;        // Not literally true; this is a "magic value".
                break;
            case btypeBase:
            case btypeTip:
            {
                // MDEBUG We assume that we won't find a Tip; this will need to be fixed when serial sampling
                // is added.  Unknown event type; we don't know how many children or parents this event type has.
                assert(false);
                string msg = " Unrecognized branch type (\""
                    +ToString(pBranch->Event())
                    + "\") received by Tree::Groom().";
                throw implementation_error(msg);
                break;
            }
        }
        double intervalLength = pBranch->m_eventTime - starttime;
        if (intervalLength <= 0.0)
        {
            if (0.0 == intervalLength && btypeRec == pBranch->Event() &&
                prevIntervalLength > 0.0)
            {
                // This is perfectly okay--rec branches come in pairs w/equal timestamps.
                lineages.UpdateBranchCounts(pBranch->m_partitions, true);
                // Recombination:  The opposite of a coalescence--it has one child and two parents.
                // Each parent branch comes through separately and gets updated.  We need to avoid
                // updating their shared child twice, so we choose to update the child when the
                // "left" parent comes through.  (nchildren was set to the "magic value" of 0.)
                if (pBranch->Child(0)->Parent(0) == pBranch)
                {
                    lineages.UpdateBranchCounts(pBranch->Child(0)->m_partitions, false);
                }
                starttime = pBranch->m_eventTime;
                prevIntervalLength = intervalLength;
                continue;
            }
            if (0.0 == intervalLength && btypeEpoch == pBranch->Event())
            {
                lineages.UpdateBranchCounts(pBranch->m_partitions, true);
                lineages.UpdateBranchCounts(pBranch->Child(0)->m_partitions, false);
                starttime = pBranch->m_eventTime;
                prevIntervalLength = intervalLength;
                continue;
            }
            string msg = "Tree::GroomForGrowth(), encountered a time interval length of ";
            msg += ToString(intervalLength) + ", followed by an event of type ";
            msg += ToString(pBranch->Event()) + ".";
            throw impossible_error(msg);
        }

        // Calculate and store E[dt] for each population.  Calculate E[dt] for coalescence + growth,
        // regardless of which event occurs at the bottom of the interval, because coalescence
        // + growth has the most extreme/sensitive values for lnWait() and DlnWait().
        for (unsigned long pop = 0; pop < npops; pop++)
        {
            const unsigned long k = xpartitions[pop];
            if (k > 1 && growths[pop] > 0.0)
            {
                expectedLength =
                    (1.0 / growths[pop])
                    *
                    ExpE1(SafeProductWithExp(k * (k - 1) / (growths[pop] * thetas[pop]), growths[pop] * starttime));
            }
            else
            {
                continue;
            }

            if (expectedLength < smallestExpectedLength)
            {
                smallestExpectedLength = expectedLength;
                worst_pop = pop;        // For debugging message only.
            }
        }

        // Note:  Calling an interval fatal if it's 100 times its new expectation value is handwaving. So far,
        // 100 seems to be OK, but this value might need to be reduced to 20 or 10 if growth is very large.
        if (smallestExpectedLength < DBL_BIG && intervalLength >= 100.0 * smallestExpectedLength)
        {
            msg = "\nTemperature " + ToString(temperature) + ", population "
                + ToString(worst_pop) + ", k = " + ToString(xpartitions[worst_pop])
                + ", interval length is " + ToString(intervalLength)
                + ", expected value is " + ToString(smallestExpectedLength)
                + ".  Interval ends with a " + ToString(pBranch->Event())
                + " event.  Time stamps are " + ToString(starttime) + " and "
                + ToString(pBranch->m_eventTime) + "; theta = " + ToString(thetas[worst_pop])
                + " and g = " + ToString(growths[worst_pop]) + ".  ";

            // Shrink this overlong interval to a value that is between 1/2 and 3/2 times the
            // expectation value under the new parameter values for the "worst" population.
            double factor = 0.5 + m_randomSource->Float();
            double newIntervalLength = factor * smallestExpectedLength;

            if (newIntervalLength < starttime * numeric_limits<double>::epsilon())
            {
                msg += "Tried to shrink this interval to a length of "
                    + ToString(newIntervalLength) + ", but this is smaller than the "
                    + "minimum length of "
                    + ToString(starttime * numeric_limits<double>::epsilon())
                    + " that can be added to " + ToString(starttime)
                    + " without being lost to rounding error.  Giving up and "
                    + "copying the cold tree into the tree for temperature "
                    + ToString(temperature) + "....\n";
                registry.GetRunReport().ReportDebug(msg);
                return false;
            }

            pBranch->m_eventTime = starttime + newIntervalLength;
            pBranch->SetUpdateDL();
            treeWasModified = true;

            msg += "Shrinking this interval to a length of "
                + ToString(newIntervalLength) + ", with a new ending time "
                + "stamp of " + ToString(pBranch->m_eventTime) + ".\n";
            registry.GetRunReport().ReportDebug(msg);

            // Once we shrink an interval, we need to shift all intervals below it to reflect this.
            cumulativeShift += intervalLength - newIntervalLength;

        }

        // Prepare for next loop iteration (next branch in the timelist).
        for (unsigned long i = 0; i < nchildren; i++)
            lineages.UpdateBranchCounts(pBranch->Child(i)->m_partitions, false);
        lineages.UpdateBranchCounts(pBranch->m_partitions, true);

        // Recombination:  The opposite of a coalescence--it has one child and two parents.
        // Each parent branch comes through separately and gets updated.  We need to avoid updating
        // their shared child twice, so we choose to update the child when the "left" parent comes
        // through.  (nchildren was set to the "magic value" of 0.)
        if (btypeRec == pBranch->Event() && pBranch->Child(0)->Parent(0) == pBranch)
        {
            lineages.UpdateBranchCounts(pBranch->Child(0)->m_partitions, false);
        }

        starttime = pBranch->m_eventTime;
        prevIntervalLength = intervalLength; // Used to detect successive zero-lengths, if any.
    }

    if (treeWasModified)
        CalculateDataLikes();

    return true;
}

//------------------------------------------------------------------------------------
// Method to, if necessary, shrink interval lengths, usually at the root, in an attempt to prevent
// explosions in the next maximization. See the extensive comment in ChainManager::DoChain(),
// which is where this method is called.

bool Tree::GroomForLogisticSelection(const DoubleVec1d & thetas,
                                     double s,           // logistic selection coefficient
                                     double temperature) // temperature used only for debugging

{
    if (0.0 == s)
        return true;
    if (2 != thetas.size() || 0.0 == thetas[0] || 0.0 == thetas[1])
    {
        string msg = "Tree:GroomForLogisticSelection(), received invalid Theta vector.";
        throw implementation_error(msg);
    }

    bool treeWasModified = false;
    double starttime(0.0), cumulativeShift(0.0), shrinkageFactor(1.0);
    double max_allowed_starttime(DBL_BIG);
    Branchiter brit, firstBadInterval, lastBadInterval;
    string msg;

    if (s > 0.0)
        max_allowed_starttime = (EXPMAX - log(thetas[1])) / s;
    else
        max_allowed_starttime = (EXPMAX - log(thetas[0])) / (-s);

    firstBadInterval = lastBadInterval = m_timeList.EndBranch();

    // Iterate through the timelist, looking for fatal intervals.
    for (brit = m_timeList.FirstBody(); brit != m_timeList.EndBranch();
         brit = m_timeList.NextBody(brit))
    {
        if ((*brit)->m_eventTime >= max_allowed_starttime)
        {
            Branchiter brit2 = brit;
            firstBadInterval = brit;
            while (brit2 != m_timeList.EndBranch())
            {
                lastBadInterval = brit2;
                brit2 = m_timeList.NextBody(brit2);
            }
            break;
        }
        starttime = (*brit)->m_eventTime;
    }

    if (firstBadInterval == m_timeList.EndBranch())
        return true;                    // nothing to shrink

    if ((*lastBadInterval)->m_eventTime - starttime <= 0.0)
    {
        msg += "Tree::GroomForLogisticSelection(), encountered an interval ";
        msg += "(or sum of intervals) of length ";
        msg += ToString((*lastBadInterval)->m_eventTime - starttime) + ".";
        throw implementation_error(msg);
    }

    shrinkageFactor = (max_allowed_starttime - starttime) /
        ((*lastBadInterval)->m_eventTime - starttime);

    if (shrinkageFactor <= 0.0)
    {
        msg += "Tree::GroomForLogisticSelection(), computed an invalid shrinkage ";
        msg += "factor of " + ToString(shrinkageFactor) + ".";
        throw implementation_error(msg);
    }

    // If there are any fatal intervals, iterate through these, and shrink them.
    for (brit = firstBadInterval; brit != m_timeList.EndBranch();
         brit = m_timeList.NextBody(brit))
    {
        Branch_ptr pBranch = *brit;

        // Once we shrink an interval, we need to shift all intervals below it by this amount.
        // The amount accumulates.
        if (cumulativeShift > 0.0)
        {
            pBranch->m_eventTime -= cumulativeShift;
            pBranch->SetUpdateDL();
        }

        double intervalLength = pBranch->m_eventTime - starttime;
        pBranch->m_eventTime = starttime + intervalLength * shrinkageFactor;
        pBranch->SetUpdateDL();
        treeWasModified = true;

        // Once we shrink an interval, we need to shift all intervals below it to reflect this.
        cumulativeShift += intervalLength - intervalLength * shrinkageFactor;

        starttime = pBranch->m_eventTime;
    }

    if (treeWasModified)
        CalculateDataLikes();

    return true;
}

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

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

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

Tree * PlainTree::MakeStump() const
{
    return new PlainTree(*this, true);
} // PlainTree stump maker

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

void PlainTree::CalculateDataLikes()
{
    m_overallDL = CalculateDataLikesForFixedLoci();
    m_timeList.ClearUpdateDLs();        // Reset the updating flags.
}

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

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

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

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

//____________________________________________________________________________________