src/tree/branch.cpp

// $Id: branch.cpp,v 1.100 2011/12/15 20:44:02 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 <algorithm>
#include <numeric>                      // for std::accumulate
#include <iostream>                     // for debug cout

#include "local_build.h"

#include "branch.h"
#include "treesum.h"
#include "summary.h"
#include "datapack.h"
#include "branchbuffer.h"

#include "force.h"                      // for RBranch::RecCopyPartitionsFrom Branch::IsAMember
#include "locus.h"                      // for TipData stuff used in constructor

#include "tinyxml.h"

#ifdef DMALLOC_FUNC_CHECK
#include <dmalloc.h>
#endif

// This turns on the detailed print out of the creation and analysis of each coalescence tree.
// JRM 4/10
//#define PRINT_TREE_DETAILS

// Print details of what was used in the summary scoring.
// JRM 4/10
//#define PRINT_SUMMARY_DETAILS

using namespace std;

// Note that some functions/variables defined here have "_Br" (for "Branch") 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".

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

// Initialization of static variable; needs to be in .cpp file.  It is being initialized to the result
// of calling the implicit default constructor for Branch_ptr (ie, boost::shared_ptr<Branch>).
Branch_ptr Branch::NONBRANCH;

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

string ToString(branch_type btype)
{
    switch(btype)
    {
        case btypeBase:   return "Base";
        case btypeTip:    return "Tip";
        case btypeCoal:   return "Coalescence";
        case btypeMig:    return "Migration";
        case btypeDivMig: return "DivMigration";
        case btypeDisease:return "Disease";
        case btypeRec:    return "Recombination";
        case btypeEpoch:  return "Epochboundary";
    }

    assert(false);
    return "";
}

//------------------------------------------------------------------------------------
// Initializes Branch's Range pointer object to pointer value of argument, not to pointer to copy of that object.
// Thus, be sure that caller of this constructor either constructs a new Range or constructs a new copy of one.
// THIS CTOR does a SHALLOW COPY of its RANGE-OBJECT POINTER.
//
// This constructor is the only one which behaves this way.  For all classes derived from Branch, constructors which
// take as a single argument or as one or more of several arguments a const Range * const value must be given a
// pointer to a freshly-allocated or copied (by a deep-copying constructor) object.  ALL OTHER CTORS (for classes
// derived from Branch) do a DEEP COPY of their RANGE-OBJECT POINTER-valued argument(s).  They usually do so by
// allocating (using CreateRange) or copying (using a deep-copying copy constructor) and passing the result
// to THIS constructor.

Branch::Branch(Range * newrangeptr)
    : boost::enable_shared_from_this<Branch>(),
      m_parents(NELEM, Branch::NONBRANCH),
      m_children(NELEM, Branch::NONBRANCH),
      m_ID(),
      m_equivBranch(Branch::NONBRANCH),
      m_partitions(registry.GetDataPack().GetNPartitionForces(), FLAGLONG),
      // Initialize m_rangePtr to pointer to original object (having been allocated by caller), not to pointer
      // to COPY of pointee.  Ie, this is a shallow copy.  All callers of this constructor must freshly-allocate
      // or deep-copy the Range object whose pointer they pass to this function.
      m_rangePtr(newrangeptr)
{
    m_eventTime     = 0.0;
    m_updateDL      = false;
    m_marked        = false;
    m_wasCoalCalced = false;
    m_isSample      = 1;

} // Branch constructor

//------------------------------------------------------------------------------------
// Must delete the Range object held by pointer and always allocated on the heap.
// All classes derived from Branch have default (ie, empty) destructors, but the
// runtime systems calls this one for the base class (since all Branch classes
// derive from Branch), and this single destructor deallocates the Range object.

Branch::~Branch()
{
    // We need to delete the contained Range member, since we own it.
    delete m_rangePtr;

} // Branch destructor

//------------------------------------------------------------------------------------
// DEEP-COPYING copy constructor.

Branch::Branch(const Branch & src)
    : boost::enable_shared_from_this<Branch>(),
      m_parents(NELEM, Branch::NONBRANCH),
      m_children(NELEM, Branch::NONBRANCH),
      m_equivBranch(Branch::NONBRANCH)
{
    CopyAllMembers(src);                // Does a deep copy of objects (namely Range) held by pointer.
} // Branch copy constructor

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

void Branch::CopyAllMembers(const Branch & src)
{
    fill(m_parents.begin(), m_parents.end(), Branch::NONBRANCH);
    fill(m_children.begin(), m_children.end(), Branch::NONBRANCH);

    m_ID               = src.m_ID;
    m_eventTime        = src.m_eventTime;
    m_partitions       = src.m_partitions;
    m_updateDL         = src.m_updateDL;
    m_marked           = src.m_marked;
    m_rangePtr         = src.m_rangePtr->Clone(); // Here's the deep copy mentioned above.
    m_DLcells          = src.m_DLcells;
    m_movingDLcells    = src.m_movingDLcells;
    m_isSample         = src.m_isSample;
    m_wasCoalCalced    = src.m_wasCoalCalced;

} // Branch::CopyAllMembers

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

bool Branch::IsAMember(const Force & force, const LongVec1d & membership) const
{
    return force.IsAMember(m_partitions, membership);
} // Branch::IsAMember

//------------------------------------------------------------------------------------
// Checks if two branches are functionally the same.
// If you want full equivalency use operator==().

bool Branch::IsEquivalentTo(const Branch_ptr twin)   const
{
    return m_ID == twin->m_ID;

#if 0  // Useful for debugging, perhaps.
    if (Event() != twin->Event())
    {
        // cout << "Different events" << endl;
        return false;
    }

    if (m_partitions != twin->m_partitions)
    {
        // cout << "Different partitions" << endl;
        return false;
    }

    if (m_rangePtr->GetLiveSites_Rg() != twin->m_rangePtr->GetLiveSites_Rg())
    {
        // cout << "Different Active sites" << endl;
        return false;
    }

    if (m_rangePtr->GetRecSite_Rg() != twin->m_rangePtr->GetRecSite_Rg())
    {
        // cout << "Different recombination site" << endl;
        return false;
    }

    if (m_eventTime != twin->m_eventTime)
    {
        // cout << "Different event times" << endl;
        return false;
    }

    return true;
#endif

} // IsEquivalentTo

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

bool Branch::HasSamePartitionsAs(const Branch_ptr twin)   const
{
    if (twin)
    {
        if (m_partitions != twin->m_partitions)
        {
            return false;
        }
    }

    return true;

} //HasSamePartitionsAs

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

bool Branch::PartitionsConsistentWith(const Branch_ptr child)   const
{
    assert(Event() == btypeRec);
    if (child)
    {
        const ForceSummary & fs = registry.GetForceSummary();
        if (fs.CheckForce(force_DISEASE))
        {
            if (m_rangePtr->IsDiseaseSiteTransmitted_Rg())
            {
                long diseaseIndex = fs.GetPartIndex(force_DISEASE);
                long diseaseState = child->GetPartition(force_DISEASE);
                if (diseaseState != m_partitions[diseaseIndex])
                {
                    return false;
                }
            }
        }
    }

    return true;

} // PartitionsConsistentWith

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

void Branch::ResetBuffersForNextRearrangement()
{
    m_rangePtr->SetOldInfoToCurrent_Rg();
    m_rangePtr->ClearNewTargetLinks_Rg();

} // ResetBuffersForNextRearrangement

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

LongVec1d Branch::GetLocalPartitions() const
{
    LongVec1d lppartitions;
    LongVec1d lpindex(registry.GetForceSummary().GetLocalPartitionIndexes());
    LongVec1d::iterator lpforce;

    for(lpforce = lpindex.begin(); lpforce != lpindex.end(); ++lpforce)
    {
        lppartitions.push_back(m_partitions[*lpforce]);
    }

    return lppartitions;

} // GetLocalPartitions

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

long Branch::GetID() const
{
    return m_ID.ID();
}

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

weakBranch_ptr Branch::GetEquivBranch() const
{
    return m_equivBranch;
}

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

void Branch::SetEquivBranch(Branch_ptr twin)
{
    m_equivBranch = twin;
}

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

long Branch::GetPartition(force_type force) const
{
    return m_partitions[registry.GetForceSummary().GetPartIndex(force)];
} // Branch::GetPartition

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

void Branch::SetPartition(force_type force, long val)
{
    m_partitions[registry.GetForceSummary().GetPartIndex(force)] = val;
}

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

void Branch::CopyPartitionsFrom(Branch_ptr src)
{
    m_partitions = src->m_partitions;
} // Branch::CopyPartitionsFrom

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

void Branch::MarkParentsForDLCalc()
{
    // Evil kludge necessary because we force parents/children to always have two spaces even when empty!
    if (!m_parents[0].expired())
    {
        Branch_ptr parent0(m_parents[0]);
        if (!parent0->m_updateDL)
        {
            parent0->SetUpdateDL();
            parent0->MarkParentsForDLCalc();
        }
    }

    if (!m_parents[1].expired())
    {
        Branch_ptr parent1(m_parents[1]);
        if (!parent1->m_updateDL)
        {
            parent1->SetUpdateDL();
            parent1->MarkParentsForDLCalc();
        }
    }

} // Branch::MarkParentsForDLCalc

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

void Branch::ReplaceChild(Branch_ptr, Branch_ptr newchild)
{
    // This version is used by MBranch and DBranch; there are overrides for other branches.
    newchild->m_parents[0] = shared_from_this();
    m_children[0]          = newchild;

} // Branch::ReplaceChild

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

bool Branch::HasSameActive(const Branch & br)
{
    return m_rangePtr->SameLiveSites_Rg(br.m_rangePtr->GetLiveSites_Rg());
} // Branch::HasSameActive

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

const Cell_ptr Branch::GetDLCell(long loc, long ind, bool moving) const
{
    if (moving)
    {
        assert(loc < static_cast<long>(m_movingDLcells.size()));
        assert(ind < static_cast<long>(m_movingDLcells[loc].size()));
        return m_movingDLcells[loc][ind];
    }
    else
    {
        assert(loc < static_cast<long>(m_DLcells.size()));
        assert(ind < static_cast<long>(m_DLcells[loc].size()));
        return m_DLcells[loc][ind];
    }
}

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

Cell_ptr Branch::GetDLCell(long loc, long ind, bool moving)
{
    if (moving)
    {
        assert(loc < static_cast<long>(m_movingDLcells.size()));
        assert(ind < static_cast<long>(m_movingDLcells[loc].size()));
        return m_movingDLcells[loc][ind];
    }
    else
    {
        assert(loc < static_cast<long>(m_DLcells.size()));
        assert(ind < static_cast<long>(m_DLcells[loc].size()));
        return m_DLcells[loc][ind];
    }
}

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

double Branch::HowFarTo(const Branch & br) const
{
    return fabs(br.m_eventTime - m_eventTime);
} // Branch::HowFarTo

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

Branch_ptr Branch::GetValidChild(Branch_ptr br, long whichpos)
{
#ifdef PRINT_TREE_DETAILS
    cout << " In GetValidChild whichpos: " << whichpos << "  br: " << br << endl;
#endif

    Branch_ptr pChild = br->GetActiveChild(whichpos);
    if (pChild)
        br = GetValidChild(pChild, whichpos);

#ifdef PRINT_TREE_DETAILS
    cout << " return from GetValidChild br: " << br << endl;
#endif

    return br;

} // Branch::GetValidChild

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

Branch_ptr Branch::GetValidPanelChild(Branch_ptr br, long whichpos)
{
#ifdef PRINT_TREE_DETAILS
    cout << " In GetValidChild whichpos: " << whichpos << "  br: " << br << endl;
#endif

    Branch_ptr pChild = br->GetActivePanelChild(whichpos);
    if (pChild)
        br = GetValidPanelChild(pChild, whichpos);

#ifdef PRINT_TREE_DETAILS
    cout << " return from GetValidChild br: " << br << endl;
#endif

    return br;

} // Branch::GetValidPanelChild

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

Branch_ptr Branch::GetValidParent(long whichpos)
{
    // We are looking for an ancestor of our branch at which the site of interest coalesces.
    Branch_ptr pParent(Parent(0));
    if (Parent(1))
    {
        // This branch has two parents; which one is needed?
        if (Parent(1)->m_rangePtr->IsSiteLive_Rg(whichpos)) pParent = Parent(1);
    }

    // If we have reached the bottom of the tree, we give up--there is no parent.
    if (pParent->Event() == btypeBase) return Branch::NONBRANCH;

    // Otherwise, check if this could be the parent we need.
    if (pParent->Event() == btypeCoal)
    {
        // This could be the coalescence we're looking for, but does our site actually coalesce here?
        if (pParent->Child(0)->m_rangePtr->IsSiteLive_Rg(whichpos) &&
            pParent->Child(1)->m_rangePtr->IsSiteLive_Rg(whichpos))
        {
            return pParent;
        }
    }

    // It's not a coalescence (or the right coalescence), so keep going downward.
    return pParent->GetValidParent(whichpos);

} // GetValidParent

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

bool Branch::DiffersInDLFrom(Branch_ptr branch, long locus, long marker) const
{
    return m_DLcells[locus].DiffersFrom(branch->m_DLcells[locus], marker);
} // Branch::DiffersInDLFrom

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

bool Branch::CheckInvariant() const
{
    // Check correct time relationships among parents and offspring.
    long index;
    for (index = 0; index < NELEM; ++index)
    {
        // If the child exists it must be earlier.
        if (Child(index))
        {
            if (Child(index)->m_eventTime > m_eventTime)
            {
                return false;
            }
        }

        // If the parent exists it must be later.
        if (Parent(index))
        {
            if (Parent(index)->m_eventTime < m_eventTime)
            {
                return false;
            }
        }
    }

    return true;

} // CheckInvariant

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

bool Branch::operator==(const Branch & src) const
{
    if (Event() != src.Event()) return false;
    if (m_partitions != src.m_partitions) return false;
    if (m_eventTime != src.m_eventTime) return false;

    return true;

} // operator==

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

string Branch::DLCheck(const Branch & other) const
{
    unsigned long locus;
    string problems;

    if (m_DLcells.size() != other.m_DLcells.size())
    {
        return string("   The DLCells are different sizes--error\n");
    }

    for (locus = 0; locus < m_DLcells.size(); ++locus)
    {
        unsigned long ncells = m_DLcells[locus].size();
        if (ncells != other.m_DLcells[locus].size())
        {
            return string("   Bad branch comparison--error\n");
        }

        unsigned long ind;
        for (ind = 0; ind < ncells; ++ind)
        {
            long badmarker = m_DLcells[locus][ind]->DiffersFrom(other.m_DLcells[locus][ind]);
            if (badmarker == FLAGLONG) continue;

            problems += "   Branch " + ToString(m_ID.ID());
            problems += " differs from other branch " + ToString(other.m_ID.ID());
            problems += " at marker " + ToString(badmarker);
            problems += "\n";

            return problems;
        }
    }

    return problems;

} // Branch::DLCheck

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

void Branch::PrintInfo_Br() const
{
    cout << "ID: " << m_ID.ID() << " )" << ToString(Event()) << ")" << endl;
    cout << "Partitions:  ";

    for (unsigned long i = 0; i < m_partitions.size(); i++)
    {
        cout << m_partitions[i] << " ";
    }

    cout << endl;
    cout << "Event time:  " << m_eventTime << endl;
    m_rangePtr->PrintInfo_Rg();
}

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

vector<Branch_ptr> Branch::GetBranchChildren()
{
    vector<Branch_ptr> newkids;
    vector<weakBranch_ptr>::iterator kid;

    for(kid = m_children.begin(); kid != m_children.end(); ++kid)
    {
        newkids.push_back(kid->lock());
    }

    return newkids;
}

//------------------------------------------------------------------------------------
// This function is non-const because we need boost::shared_from_this()!

bool Branch::ConnectedTo(const Branch_ptr family)
{
    if (Child(0) && Child(0) == family) return true;
    if (Child(1) && Child(1) == family) return true;
    if (Parent(0) && Parent(0) == family) return true;
    if (Parent(1) && Parent(1) == family) return true;

    cout << endl << "Branch " << family->m_ID.ID();
    cout << " is not connected to branch " << m_ID.ID() << " ";
    cout << "even though it thinks it is via it's ";

    Branch_ptr me(shared_from_this());
    if (family->Child(0) && family->Child(0) == me) cout << "child0";
    if (family->Child(1) && family->Child(1) == me) cout << "child1";
    if (family->Parent(0) && family->Parent(0) == me) cout << "parent0";
    if (family->Parent(1) && family->Parent(1) == me) cout << "parent1";

    cout << endl;
    return false;
}

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

bool Branch::IsSameExceptForTimes(const Branch_ptr other) const
{
    if (Event() != other->Event()) return false;
    if (m_partitions != other->m_partitions) return false;

    return true;
}

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

bool Branch::RevalidateRange(FC_Status &) const
{
    return m_rangePtr->SameAsChild_Rg(Child(0)->m_rangePtr);
} // Branch::RevalidateRange

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

void Branch::AddGraphML(TiXmlElement * elem) const
{
    ////////////////////////////////////////////////////////////
    // node itself
    long myID = GetID();
    long canonicalID = GetCanonicalID();
    if (canonicalID == myID)
    {
        TiXmlElement * node = new TiXmlElement("node");
        node->SetAttribute("id",ToString(canonicalID));
        elem->LinkEndChild(node);

        // type of node
        TiXmlElement * typeInfo = new TiXmlElement("data");
        node->LinkEndChild(typeInfo);
        typeInfo->SetAttribute("key","node_type");
        TiXmlText * nTypeText = new TiXmlText(GetGraphMLNodeType());
        typeInfo->LinkEndChild(nTypeText);

        // time of node
        TiXmlElement * timeInfo = new TiXmlElement("data");
        node->LinkEndChild(timeInfo);
        timeInfo->SetAttribute("key","node_time");
        TiXmlText * nTimeText = new TiXmlText(ToString(m_eventTime));
        timeInfo->LinkEndChild(nTimeText);

        AddNodeInfo(node);
    }

    ////////////////////////////////////////////////////////////
    // branch(es) themselves
    TiXmlElement * branchElem = new TiXmlElement("edge");
    elem->LinkEndChild(branchElem);
    branchElem->SetAttribute("source",ToString(GetCanonicalParentID()));
    branchElem->SetAttribute("target",ToString(canonicalID));

    // testing partitions
    TiXmlElement * partElem = new TiXmlElement("data");
    branchElem->LinkEndChild(partElem);
    partElem->SetAttribute("key","partitions");
    TiXmlText * pTypeText = new TiXmlText(ToString(m_partitions));
    partElem->LinkEndChild(pTypeText);

    // testing range info
    TiXmlElement * rangeElem = new TiXmlElement("data");
    branchElem->LinkEndChild(rangeElem);
    rangeElem->SetAttribute("key","active_sites");
    rangeset rset = m_rangePtr->GetLiveSites_Rg();
    TiXmlText * rTypeText = new TiXmlText(ToString(rset));
    rangeElem->LinkEndChild(rTypeText);

    // AddBranchInfo(branchElem);
}

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

string Branch::GetParentIDs() const
{
    string outstr = "";
    long pCount = NParents();
    for(long count = 0 ; count < pCount; count++)
    {
        const Branch_ptr p = Parent(count);
        if (p != Branch::NONBRANCH)
        {
            outstr = outstr + " " + ToString(p->GetID());
        }
    }
    return outstr;
}

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

string Branch::GetChildIDs() const
{
    string outstr = "";
    long pCount = NChildren();
    for(long count = 0 ; count < pCount; count++)
    {
        const Branch_ptr p = Child(count);
        if (p != Branch::NONBRANCH)
        {
            outstr = outstr + " " + ToString(p->GetID());
        }
    }
    return outstr;
}

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

long Branch::GetCanonicalID() const
{
    long myID = GetID();
    const Branch_ptr p = GetRecPartner();
    if (p != Branch::NONBRANCH)
    {
        long otherID = p->GetID();
        if (otherID < myID)
        {
            return otherID;
        }
    }
    return myID;
}

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

long Branch::GetCanonicalParentID() const
{
    long downID = -1;
    long pCount = NParents();
    long trueCount = 0;
    for(long count = 0 ; count < pCount; count++)
    {
        const Branch_ptr p = Parent(count);
        if (p != Branch::NONBRANCH)
        {
            long canonicalParent = p->GetCanonicalID();
            if (trueCount == 0)
            {
                downID = canonicalParent;
            }
            else
            {
                if (downID != canonicalParent)
                {
                    assert(false);
                    return -2;
                }
            }
        }
    }
    return downID;
}

//------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------
// This constructor should only be called by the TimeList constructor.
// Passes a newly-allocated Range (by pointer) to the Branch constructor, which does a shallow copy.

BBranch::BBranch()
    : Branch(CreateRange())
{
    // Deliberately blank.
} // BBranch constructor

//------------------------------------------------------------------------------------
// Returns a pointer to a newly-heap-allocated Range or RecRange object.

Range * BBranch::CreateRange() const
{
    // JNOTE: We would probably need to change the way tree construction is handled
    // to do this properly (i.e. get rid of the stump logic; see tree constructors).
    //
    // So we create a fake Range instead, watching out for validators that will likely fail on this branch!
    long int nsites(BASEBRANCH_NSITES);

    if (registry.GetForceSummary().CheckForce(force_REC)) // This force includes recombination.
    {
        rangeset diseasesites;                            // There are no disease sites.
        rangeset alllinks(MakeRangeset(0, nsites - 1));   // All links were and are targetable.
        rangeset allsites(MakeRangeset(0, nsites));       // All sites are both live and transmitted.
        return new RecRange(nsites, diseasesites, allsites, allsites, alllinks, alllinks, allsites, allsites);
    }
    else
    {
        return new Range(nsites);
    }

} // BBranch::CreateRange

//------------------------------------------------------------------------------------
// Returns a pointer to a newly-allocated Branch, which contains a deep-copied Range
// (thanks to the deep-copying copy constructor this class inherits from Branch).

Branch_ptr BBranch::Clone() const
{
    return Branch_ptr(new BBranch(*this));
} // BBranch::Clone

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

void BBranch::ScoreEvent(TreeSummary &, BranchBuffer &) const
{
    assert(false);                      // I don't think this should ever be called.
} // BBranch::ScoreEvent

//------------------------------------------------------------------------------------
// Third arg is ref for consistency with same function in other classes, which return data therein.

void BBranch::ScoreEvent(TreeSummary &, BranchBuffer &, long &) const
{
    assert(false);                      // I don't think this should ever be called.
} // BBranch::ScoreEvent

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

bool BBranch::CheckInvariant() const
{
    // Base branches have no parents...
    long index;
    for (index = 0; index < NELEM; ++index)
    {
        if (Parent(index)) return false;
    }

    // ...and one child
    if (!Child(0)) return false;
    if (Child(1)) return false;

    if (!Branch::CheckInvariant()) return false;

    return true;

} // CheckInvariant

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

string BBranch::GetGraphMLNodeType() const
{
    assert(false);
    string outstr = "Error BBranch";
    return outstr;
}

//------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------
// Creates a new TBranch object containing a newly-allocated Range object (held by pointer),
// which was allocated by CreateRange and shallow-copied by the Branch constructor.

TBranch::TBranch(const TipData & tipdata, long int nsites, const rangeset & diseasesites)
    : Branch(CreateRange(nsites, diseasesites))
{
    m_partitions = tipdata.GetBranchPartitions();
    m_label = tipdata.label;

    // This is used for SNP panel corrections.
    if (tipdata.m_source == dsource_panel)
    {
        m_isSample = 0;
    }

#if 0
    cout <<"Tip: " << tipdata.label;
    if (tipdata.m_source == dsource_panel)
        cout <<" is a panel member ";
    else
        cout <<" is a study member ";
    cout << endl;
#endif

    // WARNING:  must be changed for sequential-sampling case
    m_eventTime = 0.0;

} // TBranch constructor

//------------------------------------------------------------------------------------
// Returns a pointer to a newly-heap-allocated Range or RecRange object.

Range * TBranch::CreateRange(long int nsites, const rangeset & diseasesites) const
{
    if (registry.GetForceSummary().CheckForce(force_REC)) // This force includes recombination.
    {
        rangeset alllinks(MakeRangeset(0, nsites - 1)); // All links were and are targetable.
        rangeset allsites(MakeRangeset(0, nsites));
        return new RecRange(nsites, diseasesites, allsites, allsites, alllinks, alllinks, allsites, allsites);
    }
    else
    {
        return new Range(nsites);
    }

} // TBranch::CreateRange

//------------------------------------------------------------------------------------
// Returns a pointer to a newly-allocated Branch, which contains a deep-copied Range
// (thanks to the deep-copying copy constructor this class inherits from Branch).

Branch_ptr TBranch::Clone() const
{
    return Branch_ptr(new TBranch(*this));
} // TBranch::Clone

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

bool TBranch::CheckInvariant() const
{
    // Tip branches have no children...
    long index;
    for (index = 0; index < NELEM; ++index)
    {
        if (Child(index)) return false;
    }

    // ...and at least one parent.
    if (!Parent(0)) return false;

    if (!Branch::CheckInvariant()) return false;

    return true;

} // CheckInvariant

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

bool TBranch::operator==(const Branch & src) const
{
    return ((Branch::operator==(src)) && (m_label == dynamic_cast<const TBranch &>(src).m_label));
} // operator==

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

bool TBranch::IsSameExceptForTimes(const Branch_ptr src) const
{
    return ((Branch::IsSameExceptForTimes(src)) && (m_label == boost::dynamic_pointer_cast<const TBranch>(src)->m_label));
} // IsSameExceptForTimes

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

void TBranch::ScoreEvent(TreeSummary &, BranchBuffer &) const
{
    assert(false);                      // I don't think this should ever be called.
} // TBranch::ScoreEvent

//------------------------------------------------------------------------------------
// Third arg is ref for consistency with same function in other classes, which return data therein.

void TBranch::ScoreEvent(TreeSummary &, BranchBuffer &, long &) const
{
    assert(false);                      // I don't think this should ever be called.
} // TBranch::ScoreEvent

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

bool TBranch::RevalidateRange(FC_Status &) const
{
    if (m_rangePtr->LiveSitesOnly_Rg()) return true;

    long int nsites(m_rangePtr->NumTotalSites_Rg());
    rangeset diseasesites(m_rangePtr->GetDiseaseSites_Rg()); // Rec-only?
    // MDEBUG following is not good form nor exception-safe; consider
    // fixing
    Range * newrangeptr(CreateRange(nsites, diseasesites));
    bool matches = (*newrangeptr == *m_rangePtr);
    delete newrangeptr;

    return matches;

} // TBranch::RevalidateRange

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

void TBranch::PrintInfo_Br() const
{
    cout << "Event type:  " << ToString(Event()) << endl;
    cout << "Label: " << m_label << endl;
    cout << "ID: " << m_ID.ID() << endl;
    cout << "Partitions:  ";

    for (unsigned long i = 0; i < m_partitions.size(); i++)
    {
        cout << m_partitions[i] << " ";
    }

    cout << endl;
    cout << "Event time:  " << m_eventTime << endl;
    m_rangePtr->PrintInfo_Rg();

} // TBranch::PrintInfo_Br

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

string TBranch::GetGraphMLNodeType() const
{
    return "Tip";
}

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

void TBranch::AddNodeInfo(TiXmlElement * nodeElem) const
{
    // type of node
    TiXmlElement * labelInfo = new TiXmlElement("data");
    nodeElem->LinkEndChild(labelInfo);
    labelInfo->SetAttribute("key","node_label");
    TiXmlText * labelText = new TiXmlText(m_label);
    labelInfo->LinkEndChild(labelText);
}

//------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------
// Creates a new CBranch object containing a newly-allocated Range object (held by pointer),
// which was allocated by CreateRange and shallow-copied by the Branch constructor.

CBranch::CBranch(const Range * const child1rangeptr, const Range * const child2rangeptr,
                 bool newbranchisinactive, const rangeset & fcsites)
    : Branch(CreateRange(child1rangeptr, child2rangeptr, newbranchisinactive, fcsites))
{
    // Deliberately blank.
} // CBranch constructor

//------------------------------------------------------------------------------------
// Returns a pointer to a newly-heap-allocated Range or RecRange object.

Range * CBranch::CreateRange(const Range * const child1rangeptr, const Range * const child2rangeptr,
                             bool newbranchisinactive, const rangeset & fcsites) const
{
    // These we can grab from either child as they should be identical in these fields.
    long int nsites(child1rangeptr->NumTotalSites_Rg());

    if (registry.GetForceSummary().CheckForce(force_REC)) // This force includes recombination.
    {
        rangeset diseasesites(child1rangeptr->GetDiseaseSites_Rg());
        rangeset transmittedsites(MakeRangeset(0, nsites)); // All sites are transmitted.
        rangeset livesites(Union(child1rangeptr->GetLiveSites_Rg(), child2rangeptr->GetLiveSites_Rg()));
        rangeset targetsites(RemoveRangeFromRange(fcsites, livesites));
        targetsites = Union(diseasesites, targetsites);
        rangeset targetlinks(TargetLinksFrom(targetsites));

        rangeset oldtargetlinks, oldtargetsites, oldlivesites;
        if (newbranchisinactive)        // Copy the inactive child oldtargetlinks and oldtargetsites.
        {
            oldtargetlinks = child1rangeptr->GetOldTargetLinks_Rg();
            oldtargetsites = child1rangeptr->GetOldTargetSites_Rg();
            oldlivesites = child1rangeptr->GetOldLiveSites_Rg();
        }
        else
        {
            oldtargetlinks = targetlinks;
            oldtargetsites = targetsites;
            oldlivesites = livesites;
        }

        return new RecRange(nsites, diseasesites, transmittedsites, livesites,
                            targetlinks, oldtargetlinks, oldtargetsites, oldlivesites);
    }
    else
    {
        return new Range(nsites);
    }

} // CBranch::CreateRange

//------------------------------------------------------------------------------------
// Returns a pointer to a newly-allocated Branch, which contains a deep-copied Range
// (thanks to the deep-copying copy constructor this class inherits from Branch).

Branch_ptr CBranch::Clone() const
{
    return Branch_ptr(new CBranch(*this));
} // CBranch::Clone

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

bool CBranch::CanRemove(Branch_ptr checkchild)
{
    if (m_marked) return true;

    m_marked = true;
    if (Child(0) == checkchild) SetChild(0, Child(1));
    SetChild(1, Branch::NONBRANCH);

    return false;

} // CBranch::CanRemove

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

void CBranch::UpdateBranchRange(const rangeset & fcsites, bool dofc)
{
    if (m_marked)
    {
        m_rangePtr->UpdateOneLeggedCRange_Rg(Child(0)->m_rangePtr);
    }
    else
    {
        m_rangePtr->UpdateCRange_Rg(Child(0)->m_rangePtr, Child(1)->m_rangePtr, fcsites, dofc);
    }

} // CBranch::UpdateBranchRange

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

void CBranch::UpdateRootBranchRange(const rangeset & fcsites, bool dofc)
{
    if (m_marked)
    {
        m_rangePtr->UpdateOneLeggedRootRange_Rg(Child(0)->m_rangePtr);
    }
    else
    {
        m_rangePtr->UpdateRootRange_Rg(Child(0)->m_rangePtr, Child(1)->m_rangePtr, fcsites, dofc);
    }

} // CBranch::UpdateRootBranchRange

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

void CBranch::ReplaceChild(Branch_ptr pOldChild, Branch_ptr pNewChild)
{
    if (Child(0) == pOldChild)
        SetChild(0, pNewChild);
    else
        SetChild(1, pNewChild);

    pNewChild->SetParent(0, shared_from_this());

} // CBranch::ReplaceChild

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

Branch_ptr CBranch::OtherChild(Branch_ptr badchild)
{
    if (Child(0) == badchild) return Child(1);
    return Child(0);

} // CBranch::OtherChild

//------------------------------------------------------------------------------------
// Both children must be both active and included in the DLcalc for this to return true.

bool CBranch::CanCalcDL(long site) const
{
    if ((Child(0)->m_rangePtr->IsSiteLive_Rg(site)) && (Child(1)->m_rangePtr->IsSiteLive_Rg(site)))
    {
        return true;
    }
    else
    {
        return false;
    }
}

//------------------------------------------------------------------------------------
// If both children are active return Branch::NONBRANCH to signal a stop; otherwise return the active child.

const Branch_ptr CBranch::GetActiveChild(long site)  const
{
    if (Child(0)->m_rangePtr->IsSiteLive_Rg(site))
    {
        if (Child(1)->m_rangePtr->IsSiteLive_Rg(site)) return Branch::NONBRANCH;
        return Child(0);
    }
    else
    {
        return Child(1);
    }

} // CBranch::GetActiveChild

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

void CBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks) const
{
    long i;
    Summary * csum = summary.GetSummary(force_COAL);

    // Forces have become inconsistent!
    assert(dynamic_cast<CoalSummary *>(csum));

    long myxpart = registry.GetDataPack().GetCrossPartitionIndex(m_partitions);
    LongVec1d emptyvec;

    // Score the event.
    csum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
                      FLAGLONG, myxpart, FLAGLONG, FLAGLONG, emptyvec, force_COAL);

    // Adjust the branch counts.
    for(i = 0; i < NELEM; ++i)
        ks.UpdateBranchCounts(Child(i)->m_partitions, false);

    ks.UpdateBranchCounts(m_partitions);

} // CBranch::ScoreEvent

//------------------------------------------------------------------------------------
// Third arg is a reference because this function returns results therein.

void CBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks, long & s) const
{
    long i;
    Summary * csum = summary.GetSummary(force_COAL);

    // Forces have become inconsistent!
    assert(dynamic_cast<CoalSummary *>(csum));

    long myxpart = registry.GetDataPack().
        GetCrossPartitionIndex(m_partitions);
    LongVec1d emptyvec;

    // Score the event.
    csum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
                      s, myxpart, FLAGLONG, FLAGLONG, emptyvec, force_COAL);

    // Rest for recominant branches only? - BOBGIAN
    // Adjust the branch and link counts.
    for (i = 0; i < NELEM; ++i)
    {
        ks.UpdateBranchCounts(Child(i)->m_partitions, false);
        s -= Child(i)->m_rangePtr->NumTargetLinks_Rg();
    }

    ks.UpdateBranchCounts(m_partitions);
    s += m_rangePtr->NumTargetLinks_Rg();

} // CBranch::ScoreEvent

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

bool CBranch::CheckInvariant() const
{
    // Coalescent branches have two children...
    if (!Child(0)) return false;
    if (!Child(1)) return false;

    //...and at least one parent.
    if (!Parent(0)) return false;

    if (!Branch::CheckInvariant()) return false;

    return true;

} // CheckInvariant

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

bool CBranch::RevalidateRange(FC_Status & fcstatus) const
{
    rangeset livesites;

    // Are we the same moeity and population as our children?
    if (Child(0)->m_partitions != m_partitions)
    {
        cout << "Branch changed color in CBranch" << endl;
        return false;
    }

    if (Child(1) && Child(1)->m_partitions != m_partitions)
    {
        cout << "Branch changed color in CBranch" << endl;
        return false;
    }

    if (!Child(1))                      // We're in a one-legger.
    {
        const Range * const childrangeptr(Child(0)->m_rangePtr);
        livesites = childrangeptr->GetLiveSites_Rg();
        if (!(m_rangePtr->SameLiveSites_Rg(livesites)))
        {
            return false;
        }
        // We don't need to update fc stuff for a one-legger.
    }
    else
    {
        // We're in a normal coalescence.
        const Range * const child0rangeptr(Child(0)->m_rangePtr);
        const Range * const child1rangeptr(Child(1)->m_rangePtr);

        livesites = Union(child0rangeptr->GetLiveSites_Rg(), child1rangeptr->GetLiveSites_Rg());
        if (!(m_rangePtr->SameLiveSites_Rg(livesites)))
        {
            return false;
        }

#if LAMARC_QA_FC_ON
        rangeset fcsites(Intersection(child0rangeptr->GetLiveSites_Rg(), child1rangeptr->GetLiveSites_Rg()));
        fcstatus.Decrement_FC_Counts(fcsites);
#endif
    }

    if (m_rangePtr->LiveSitesOnly_Rg()) return true; // We're done checking.

    // Will rest of this function only execute on recombinant branches?
#if LAMARC_QA_FC_ON
    rangeset target(TargetLinksFrom(Union(m_rangePtr->GetDiseaseSites_Rg(),
                                          RemoveRangeFromRange(fcstatus.Coalesced_Sites(), livesites))));
#else
    rangeset target(TargetLinksFrom(Union(m_rangePtr->GetDiseaseSites_Rg(), livesites)));
#endif

    if (m_rangePtr->DifferentTargetLinks_Rg(target))
    {
        cout << "my range; ";
        m_rangePtr->PrintInfo_Rg();
        cout << "target: " << ToString(target) << endl;
        return false;
    }

    // Is this a rec-only branch?  Otherwise, we may not want GetOldTargetLinks_Rg() to ASSERT.
    rangeset newtarget(RemoveRangeFromRange(m_rangePtr->GetOldTargetLinks_Rg(), target));
    if (m_rangePtr->DifferentNewTargetLinks_Rg(newtarget))
    {
        return false;
    }

    return true;

} // CBranch::RevalidateRange

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

string CBranch::GetGraphMLNodeType() const
{
    return "Coal";
}

//------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------
// Creates a new PartitionBranch object containing a newly-allocated Range object
// (held by pointer), which was allocated by CreateRange and shallow-copied by the Branch constructor.

PartitionBranch::PartitionBranch(const Range * const childrangeptr)
    : Branch(CreateRange(childrangeptr))
{
    // Deliberately blank.
} // PartitionBranch constructor

//------------------------------------------------------------------------------------
// Returns a pointer to a newly-heap-allocated Range or RecRange object.

Range * PartitionBranch::CreateRange(const Range * const childrangeptr) const
{
    // NOT copy-constructed from the child, because a PartitionBranch always transmits all sites, and the child might not.
    long int nsites(childrangeptr->NumTotalSites_Rg());

    if (registry.GetForceSummary().CheckForce(force_REC)) // This force includes recombination.
    {
        rangeset diseasesites(childrangeptr->GetDiseaseSites_Rg());
        rangeset transmittedsites(MakeRangeset(0, nsites));
        rangeset livesites(childrangeptr->GetLiveSites_Rg());
        rangeset targetlinks(childrangeptr->GetTargetLinks_Rg());
        rangeset oldtargetlinks(childrangeptr->GetOldTargetLinks_Rg());
        rangeset oldtargetsites(childrangeptr->GetOldTargetSites_Rg());
        rangeset oldlivesites(childrangeptr->GetOldLiveSites_Rg());

        return new RecRange(nsites, diseasesites, transmittedsites, livesites,
                            targetlinks, oldtargetlinks, oldtargetsites, oldlivesites);
    }
    else
    {
        return new Range(nsites);
    }

} // PartitionBranch::CreateRange

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

void PartitionBranch::UpdateBranchRange(const rangeset &, bool)
{
    m_rangePtr->UpdateMRange_Rg(Child(0)->m_rangePtr);
} // PartitionBranch::UpdateBranchRange

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

bool PartitionBranch::CheckInvariant() const
{
    // Partition branches have one child...
    if (!Child(0)) return false;
    if (Child(1)) return false;

    //...and at least one parent.
    if (!Parent(0)) return false;

    if (!Branch::CheckInvariant()) return false;

    return true;

} // CheckInvariant

//------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------
// Creates a new MigLikeBranch object containing a newly-allocated Range object (held by pointer),
// which was allocated by CreateRange in the PartitionBranch constructor.

MigLikeBranch::MigLikeBranch(const Range * const protorangeptr)
    : PartitionBranch(protorangeptr)
{
    // Deliberately blank.
} // MigLikeBranch constructor

//------------------------------------------------------------------------------------
// Creates a new MigLikeBranch object containing a newly-allocated Range object (held by pointer), which
// was allocated by CreateRange and shallow-copied by the Branch constructor that PartitionBranch inherits.

MigLikeBranch::MigLikeBranch(const MigLikeBranch & src)
    : PartitionBranch(src)
{
    // Deliberately blank.
} // MigLikeBranch copy constructor

//------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------
// Creates a new MBranch object containing a newly-allocated Range object (held by pointer), which
// was allocated by CreateRange and shallow-copied by the Branch constructor that MigLikeBranch inherits.

MBranch::MBranch(const Range * const protorangeptr)
    : MigLikeBranch(protorangeptr)
{
    // Deliberately blank.
} // MBranch constructor

//------------------------------------------------------------------------------------
// Creates a new MBranch object containing a newly-allocated Range object (held by pointer), which
// was allocated by CreateRange and shallow-copied by the Branch constructor that MigLikeBranch inherits.

MBranch::MBranch(const MBranch & src)
    : MigLikeBranch(src)
{
    // Deliberately blank.
} // MBranch copy constructor

//------------------------------------------------------------------------------------
// Creates a new MBranch object containing a newly-allocated Range object (held by pointer), which
// was deep-copied by the PartitionBranch copy constructor that MigLikeBranch inherits.

Branch_ptr MBranch::Clone() const
{
    return Branch_ptr(new MBranch(*this));
} // MBranch::Clone

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

void MBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks) const
{
    Summary * msum = summary.GetSummary(force_MIG);

    // MDEBUG is correct with divergence?  Forces have become inconsistent!
    assert(dynamic_cast<MigSummary *>(msum));

    long mypop = GetPartition(force_MIG);
    long chpop = Child(0)->GetPartition(force_MIG);
    LongVec1d emptyvec;

    // Score the event.
    msum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
                      FLAGLONG, mypop, chpop, FLAGLONG, emptyvec, force_MIG);

    // Adjust the branch counts.
    assert(!Child(1));                  // too many children??
    ks.UpdateBranchCounts(m_partitions);
    ks.UpdateBranchCounts(Child(0)->m_partitions, false);

    // We do not adjust active sites; they cannot possibly change.

} // MBranch::ScoreEvent

//------------------------------------------------------------------------------------
// Third arg is ref for consistency with same function in other classes, which return data therein.

void MBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks, long & s) const
{
    Summary * msum = summary.GetSummary(force_MIG);

    // Forces have become inconsistent!
    assert(dynamic_cast<MigSummary *>(msum));

    long mypop = GetPartition(force_MIG);
    long chpop = Child(0)->GetPartition(force_MIG);
    LongVec1d emptyvec;

    // Score the event.
    msum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
                      s, mypop, chpop, FLAGLONG, emptyvec, force_MIG);

    // Adjust the branch counts.
    assert(!Child(1));                  // too many children??
    ks.UpdateBranchCounts(m_partitions);
    ks.UpdateBranchCounts(Child(0)->m_partitions, false);

    // We do not adjust active sites; they cannot possibly change.

} // MBranch::ScoreEvent

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

string MBranch::GetGraphMLNodeType() const
{
    string outstr = "Mig";
    return outstr;
}

//------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------
// Divergence migration branches.  Practically identical to MBranch.

DivMigBranch::DivMigBranch(const Range * const protorangeptr)
    : MigLikeBranch(protorangeptr)
{
    // Deliberately blank.
} // DivMigBranch constructor

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

DivMigBranch::DivMigBranch(const DivMigBranch & src)
    : MigLikeBranch(src)
{
    // Deliberately blank.
} // DivMigBranch copy constructor

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

Branch_ptr DivMigBranch::Clone() const
{
    return Branch_ptr(new DivMigBranch(*this));
} // DivMigBranch::Clone

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

void DivMigBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks) const
{
    Summary * msum = summary.GetSummary(force_DIVMIG);

    // MDEBUG is correct with divergence?  Forces have become inconsistent!
    assert(dynamic_cast<DivMigSummary *>(msum));

    long mypop = GetPartition(force_DIVMIG);
    long chpop = Child(0)->GetPartition(force_DIVMIG);
    LongVec1d emptyvec;

    // Score the event.
    msum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
                      FLAGLONG, mypop, chpop, FLAGLONG, emptyvec, force_DIVMIG);

    // Adjust the branch counts.
    assert(!Child(1));                  // too many children??
    ks.UpdateBranchCounts(m_partitions);
    ks.UpdateBranchCounts(Child(0)->m_partitions, false);

    // We do not adjust active sites; they cannot possibly change.

} // DivMigBranch::ScoreEvent

//------------------------------------------------------------------------------------
// Third arg is ref for consistency with same function in other classes, which return data therein.

void DivMigBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks, long & s) const
{
    Summary * msum = summary.GetSummary(force_DIVMIG);

    // Forces have become inconsistent!
    assert(dynamic_cast<DivMigSummary *>(msum));

    long mypop = GetPartition(force_DIVMIG);
    long chpop = Child(0)->GetPartition(force_DIVMIG);
    LongVec1d emptyvec;

    // Score the event.
    msum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
                      s, mypop, chpop, FLAGLONG, emptyvec, force_DIVMIG);

    // Adjust the branch counts.
    assert(!Child(1));                  // too many children??
    ks.UpdateBranchCounts(m_partitions);
    ks.UpdateBranchCounts(Child(0)->m_partitions, false);

    // We do not adjust active sites; they cannot possibly change.

} // DivMigBranch::ScoreEvent

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

string DivMigBranch::GetGraphMLNodeType() const
{
    string outstr = "DivMig";
    return outstr;
}

//------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------
// Disease mutation branches.  NOT related to Divergence!

DBranch::DBranch(const Range * const protorangeptr)
    : PartitionBranch(protorangeptr)
{
    // Deliberately blank.
} // DBranch constructor

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

DBranch::DBranch(const DBranch & src)
    : PartitionBranch(src)
{
    // Deliberately blank.
} // DBranch copy constructor

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

Branch_ptr DBranch::Clone() const
{
    return Branch_ptr(new DBranch(*this));
} // DBranch::Clone

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

void DBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks) const
{
    Summary * dissum = summary.GetSummary(force_DISEASE);

    // Forces have become inconsistent!
    assert(dynamic_cast<DiseaseSummary *>(dissum));

    long mydis = GetPartition(force_DISEASE);
    long chdis = Child(0)->GetPartition(force_DISEASE);
    LongVec1d emptyvec;

    // Score the event.
    dissum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
                        FLAGLONG, mydis, chdis, FLAGLONG, emptyvec, force_DISEASE);

    // Adjust the branch counts.
    assert(!Child(1));                  // too many children??
    ks.UpdateBranchCounts(m_partitions);
    ks.UpdateBranchCounts(Child(0)->m_partitions, false);

    // We do not adjust active sites; they cannot possibly change.

} // DBranch::ScoreEvent

//------------------------------------------------------------------------------------
// Third arg is ref for consistency with same function in other classes, which return data therein.

void DBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks, long & s) const
{
    Summary * dissum = summary.GetSummary(force_DISEASE);

    // Forces have become inconsistent!
    assert(dynamic_cast<DiseaseSummary *>(dissum));

    long mydis = GetPartition(force_DISEASE);
    long chdis = Child(0)->GetPartition(force_DISEASE);
    LongVec1d emptyvec;

    // Score the event.
    dissum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
                        s, mydis, chdis, FLAGLONG, emptyvec, force_DISEASE);

    // Adjust the branch counts.
    assert(!Child(1));                  // too many children??
    ks.UpdateBranchCounts(m_partitions);
    ks.UpdateBranchCounts(Child(0)->m_partitions, false);

    // We do not adjust active sites; they cannot possibly change.

} // DBranch::ScoreEvent

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

string DBranch::GetGraphMLNodeType() const
{
    string outstr = "Disease";
    return outstr;
}

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

Branch_ptr EBranch::Clone() const
{
    return Branch_ptr(new EBranch(*this));
} // EBranch::Clone

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

void EBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks) const
{
    Summary * esum = summary.GetSummary(force_DIVERGENCE);

    // Forces have become inconsistent!
    assert(dynamic_cast<EpochSummary *>(esum));

    // NB While this is not a partition force, it cooperates with DIVMIG here.
    long mypop = GetPartition(force_DIVMIG);
    long chpop = Child(0)->GetPartition(force_DIVMIG);
    LongVec1d emptyvec;

    // Score the event.
    esum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
                      FLAGLONG, mypop, chpop, FLAGLONG, emptyvec, force_DIVERGENCE);

    // Adjust the branch counts.
    assert(!Child(1));                  // too many children??
    ks.UpdateBranchCounts(m_partitions);
    ks.UpdateBranchCounts(Child(0)->m_partitions, false);

    // We do not adjust active sites; they cannot possibly change.

} // EBranch::ScoreEvent

//------------------------------------------------------------------------------------
// Third arg is ref for consistency with same function in other classes, which return data therein.

void EBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks, long & s) const
{
    Summary * esum = summary.GetSummary(force_DIVERGENCE);

    // Forces have become inconsistent!
    assert(dynamic_cast<EpochSummary *>(esum));

    // NB While this is not a partition force, it cooperates with DIVMIG here.
    long mypop = GetPartition(force_DIVMIG);
    long chpop = Child(0)->GetPartition(force_DIVMIG);
    LongVec1d emptyvec;

    // Score the event.
    esum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
                      s, mypop, chpop, FLAGLONG, emptyvec, force_DIVERGENCE);

    // Adjust the branch counts.
    assert(!Child(1));                  // too many children??
    ks.UpdateBranchCounts(m_partitions);
    ks.UpdateBranchCounts(Child(0)->m_partitions, false);

    // We do not adjust active sites; they cannot possibly change.

} // EBranch::ScoreEvent

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

string EBranch::GetGraphMLNodeType() const
{
    string outstr = "Epoch";
    return outstr;
}

//------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------
// Creates a new RBranch object containing a newly-allocated Range object (held by pointer),
// which was allocated by CreateRange and shallow-copied by the Branch constructor.

RBranch::RBranch(const Range * const childrangeptr, bool newbranchisinactive,
                 const rangeset & transmittedsites, const rangeset & fcsites)
    : Branch(CreateRange(childrangeptr, newbranchisinactive, transmittedsites, fcsites))
{
    // Deliberately blank.
} // RBranch constructor

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

Branch_ptr RBranch::Clone() const
{
    return Branch_ptr(new RBranch(*this));
} // RBranch::Clone

//------------------------------------------------------------------------------------
// Returns a pointer to a newly-heap-allocated RecRange object.
// ASSERTs if attempt is made to create a Range object.

Range * RBranch::CreateRange(const Range * const childrangeptr, bool newbranchisinactive,
                             const rangeset & transmittedsites, const rangeset & fcsites) const
{
    if (registry.GetForceSummary().CheckForce(force_REC)) // This force includes recombination.
    {
        // The following are identical with our child.
        long int nsites(childrangeptr->NumTotalSites_Rg());
        rangeset diseasesites(childrangeptr->GetDiseaseSites_Rg());
        rangeset livesites(Intersection(childrangeptr->GetLiveSites_Rg(), transmittedsites));

        // MNEWCODE I am still not positive the following lines are correct!
        // But I can't find the case that would prove them wrong.
        rangeset targetsites(RemoveRangeFromRange(fcsites, livesites));
        targetsites = Union(diseasesites, targetsites);
        rangeset targetlinks(TargetLinksFrom(targetsites));

        rangeset oldtargetlinks, oldtargetsites, oldlivesites;
        if (newbranchisinactive)        // Copy the inactive child oldtargetlinks and oldtargetsites.
        {
            oldtargetlinks = childrangeptr->GetOldTargetLinks_Rg();
            oldtargetsites = childrangeptr->GetOldTargetSites_Rg();
            oldlivesites = childrangeptr->GetOldLiveSites_Rg();
        }
        else
        {
            oldtargetlinks = targetlinks;
            oldtargetsites = targetsites;
            oldlivesites = livesites;
        }

        return new RecRange(nsites, diseasesites, transmittedsites, livesites,
                            targetlinks, oldtargetlinks, oldtargetsites, oldlivesites);
    }
    else
    {
        assert(false);                  // An RBranch in a non-rec run??
        return NULL;                    // To silence compiler warning.
    }

} // RBranch::CreateRange

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

void RBranch::CopyPartitionsFrom(Branch_ptr src)
{
    m_partitions = src->m_partitions;
} // Branch::CopyPartitionsFrom

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

void RBranch::RecCopyPartitionsFrom(Branch_ptr src, FPartMap fparts, bool islow)
{
    m_partitions = src->m_partitions; //fparts may override.

    ForceVec forces(registry.GetForceSummary().GetPartitionForces());

    for(unsigned long force = 0; force < forces.size(); force++)
    {
        FPartMapiter thisforcepart = fparts.find(forces[force]->GetTag());
        if (thisforcepart != fparts.end())
        {
            long chosenpart = thisforcepart->second;
            long partition = dynamic_cast<PartitionForce *>(forces[force])->
                ChoosePartition(src->m_partitions[force], chosenpart, islow, GetRecSite_Br());
            m_partitions[force] = partition;
        }
    }

} // RBranch::RecCopyPartitionsFrom

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

void RBranch::UpdateBranchRange(const rangeset & fcsites, bool dofc)
{
    m_rangePtr->UpdateRRange_Rg(Child(0)->m_rangePtr, fcsites, dofc);
} // RBranch::UpdateBranchRange

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

void RBranch::ReplaceChild(Branch_ptr oldchild, Branch_ptr newchild)
{
    SetChild(0, newchild);

    Branch_ptr myself = shared_from_this();

    if (oldchild->Parent(0) == myself)
    {
        newchild->SetParent(0, myself);
    }
    else
    {
        newchild->SetParent(1, myself);
    }

} // RBranch::ReplaceChild

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

bool RBranch::IsRemovableRecombinationLeg(const rangeset & fcsites) const
{
    return (!m_rangePtr->AreChildTargetSitesTransmitted_Rg(Child(0)->m_rangePtr, fcsites));
} // RBranch::IsRemovableRecombinationLeg

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

bool RBranch::operator==(const Branch & other) const
{
    // Potentially troublesome call to RecRange::operator== here - BOBGIAN.
    return (Branch::operator==(other) && (*m_rangePtr == *(other.m_rangePtr)));
}

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

bool RBranch::IsSameExceptForTimes(const Branch_ptr other) const
{
    // Potentially troublesome call to RecRange::operator== here - BOBGIAN.
    return (Branch::IsSameExceptForTimes(other) && (*m_rangePtr == *(other->m_rangePtr)));
}

//------------------------------------------------------------------------------------
// Recombinant branches only (returns Branch::NONBRANCH in non-recombinant case).

Branch_ptr RBranch::GetRecPartner() const
{
    Branch_ptr partner((Child(0)->Parent(0) == shared_from_this()) ?
                       Child(0)->Parent(1) : Child(0)->Parent(0));

    return partner;

} // GetRecPartner

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

void RBranch::ScoreEvent(TreeSummary &, BranchBuffer &) const
{
    assert(false);                      // This should never be called.
} // Rbranch::ScoreEvent

//------------------------------------------------------------------------------------
// Third arg is a reference because this function returns results therein.

void RBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks, long & s) const
{
    // One interval with a recombination at the top involves *two* RBranches.  Only one
    // will be summarized into the TreeSummary.  Thus, this peculiar-looking code only calls
    // AddInterval() if the RBranch in question is its Child's *first* Parent, though
    // it does clean-up bookkeeping in any case.
    if (Child(0)->Parent(0) == shared_from_this())
    {
        RecSummary * rsum = dynamic_cast<RecSummary *>(summary.GetSummary(force_REC));
        RecTreeSummary & rtreesum = dynamic_cast<RecTreeSummary &>(summary);

        // Forces have become inconsistent!?
        assert(rsum);

        // Obtain the partnerpicks information for the branch which had it, which is not necessarily
        // this one.  MDEBUG:  This code assumes there is at most one local partition force and one
        // non-local partition force.  It will need to be more complex when multiples of either are allowed.

        // I assume that either both parents match the child, and it doesn't matter which we use,
        // or one parent does not match the child, in which case that parent must be used.

        LongVec1d childpartitions = Child(0)->m_partitions;
        LongVec1d otherpartitions = Child(0)->Parent(1)->m_partitions;  // other branch of rec
        LongVec1d partnerpicks;
        LongVec1d lpindex(registry.GetForceSummary().GetLocalPartitionIndexes());
        LongVec1d::iterator lpforce;
        LongVec1d pickedmembership;

        // MDEBUG vestigal code assuming multiple lpforces are possible here.
        // This is not carried through correctly elsewhere in function!
        for(lpforce = lpindex.begin(); lpforce != lpindex.end(); ++lpforce)
        {
            if (m_partitions[*lpforce] == childpartitions[*lpforce])
            {
                pickedmembership = otherpartitions;
            }
            else
            {
                pickedmembership = m_partitions;
            }
            partnerpicks.push_back(pickedmembership[*lpforce]);
            rsum->AddToRecombinationCounts(pickedmembership);
        }

        // Score the event.
        rsum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(), s,
                          FLAGLONG, FLAGLONG, GetRecSite_Br(), partnerpicks, force_REC);

        // Update list of recombinations by partition.
        rtreesum.AddRecToRecsByPart(pickedmembership, rsum->GetLastAdded());

        // Adjust the branch counts and activesite counts for removal of the child.
        ks.UpdateBranchCounts(Child(0)->m_partitions, false);
        s -= Child(0)->m_rangePtr->NumTargetLinks_Rg();
    }

    // Adjust the branch and activesite counts for addition of this branch.
    ks.UpdateBranchCounts(m_partitions);
    s += m_rangePtr->NumTargetLinks_Rg();

} // RBranch::ScoreEvent

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

bool RBranch::CheckInvariant() const
{
    // Recombinant branches have one child...
    if (!Child(0))
    {
        return false;
    }

    if (Child(1))
    {
        return false;
    }

    //...and at least one parent.
    if (!Parent(0))
    {
        return false;
    }

    if (!Branch::CheckInvariant()) return false;
    return true;

} // CheckInvariant

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

bool RBranch::RevalidateRange(FC_Status & fcstatus) const
{
#if 0 // Include this branch if RevalidateRange is only being called after Prune().
    long recsite(m_rangePtr->GetRecSite_Rg());
    if (!(Child(0)->m_rangePtr->IsLinkTargettable_Rg(recsite)))
    {
        cout << "Non-targettable recombination site found in RBranch" << endl;
        return false;
    }
#endif

    rangeset livesites(Intersection(Child(0)->m_rangePtr->GetLiveSites_Rg(), m_rangePtr->GetTransmittedSites_Rg()));
    if (!(m_rangePtr->SameLiveSites_Rg(livesites)))
    {
        return false;
    }

    // Did we change color inappropriately?
    if (m_rangePtr->IsDiseaseSiteTransmitted_Rg())
    {
        if (Child(0)->m_partitions != m_partitions)
        {
            cout << "Branch changed color in RBranch" << endl;
            return false;
        }
    }

    if (m_rangePtr->LiveSitesOnly_Rg()) return true; // We're done.

#if LAMARC_QA_FC_ON
    rangeset target(TargetLinksFrom(Union(m_rangePtr->GetDiseaseSites_Rg(),
                                          RemoveRangeFromRange(fcstatus.Coalesced_Sites(), livesites))));
#else
    rangeset target(TargetLinksFrom(Union(m_rangePtr->GetDiseaseSites_Rg(), livesites)));
#endif

    if (m_rangePtr->DifferentTargetLinks_Rg(target))
    {
        return false;
    }

    rangeset newtarget(RemoveRangeFromRange(m_rangePtr->GetOldTargetLinks_Rg(), target));
    if (m_rangePtr->DifferentNewTargetLinks_Rg(newtarget))
    {
        return false;
    }

    return true;

} // RBranch::RevalidateRange

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

string RBranch::GetGraphMLNodeType() const
{
    return "Rec";
}

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

void RBranch::AddNodeInfo(TiXmlElement * nodeElem) const
{
    // testing recloc info
    TiXmlElement * recElem = new TiXmlElement("data");
    nodeElem->LinkEndChild(recElem);
    recElem->SetAttribute("key","rec_location");
    long recLoc = m_rangePtr->GetRecSite_Rg();
    TiXmlText * rlTypeText = new TiXmlText(ToString(recLoc));
    recElem->LinkEndChild(rlTypeText);
}

//____________________________________________________________________________________