1 // $Id: branch.cpp,v 1.100 2011/12/15 20:44:02 bobgian Exp $
2
3 /*
4 Copyright 2002 Peter Beerli, Mary Kuhner, Jon Yamato and Joseph Felsenstein
5
6 This software is distributed free of charge for non-commercial use
7 and is copyrighted. Of course, we do not guarantee that the software
8 works, and are not responsible for any damage you may cause or have.
9 */
10
11 #include <cassert>
12 #include <algorithm>
13 #include <numeric> // for std::accumulate
14 #include <iostream> // for debug cout
15
16 #include "local_build.h"
17
18 #include "branch.h"
19 #include "treesum.h"
20 #include "summary.h"
21 #include "datapack.h"
22 #include "branchbuffer.h"
23
24 #include "force.h" // for RBranch::RecCopyPartitionsFrom Branch::IsAMember
25 #include "locus.h" // for TipData stuff used in constructor
26
27 #include "tinyxml.h"
28
29 #ifdef DMALLOC_FUNC_CHECK
30 #include <dmalloc.h>
31 #endif
32
33 // This turns on the detailed print out of the creation and analysis of each coalescence tree.
34 // JRM 4/10
35 //#define PRINT_TREE_DETAILS
36
37 // Print details of what was used in the summary scoring.
38 // JRM 4/10
39 //#define PRINT_SUMMARY_DETAILS
40
41 using namespace std;
42
43 // Note that some functions/variables defined here have "_Br" (for "Branch") as a suffix in their name.
44 // This is to distinguish them from functions/variables in class Range (or RecRange) of the same name, which use "_Rg".
45
46 //------------------------------------------------------------------------------------
47
48 // Initialization of static variable; needs to be in .cpp file. It is being initialized to the result
49 // of calling the implicit default constructor for Branch_ptr (ie, boost::shared_ptr<Branch>).
50 Branch_ptr Branch::NONBRANCH;
51
52 //------------------------------------------------------------------------------------
53
ToString(branch_type btype)54 string ToString(branch_type btype)
55 {
56 switch(btype)
57 {
58 case btypeBase: return "Base";
59 case btypeTip: return "Tip";
60 case btypeCoal: return "Coalescence";
61 case btypeMig: return "Migration";
62 case btypeDivMig: return "DivMigration";
63 case btypeDisease:return "Disease";
64 case btypeRec: return "Recombination";
65 case btypeEpoch: return "Epochboundary";
66 }
67
68 assert(false);
69 return "";
70 }
71
72 //------------------------------------------------------------------------------------
73 // Initializes Branch's Range pointer object to pointer value of argument, not to pointer to copy of that object.
74 // Thus, be sure that caller of this constructor either constructs a new Range or constructs a new copy of one.
75 // THIS CTOR does a SHALLOW COPY of its RANGE-OBJECT POINTER.
76 //
77 // This constructor is the only one which behaves this way. For all classes derived from Branch, constructors which
78 // take as a single argument or as one or more of several arguments a const Range * const value must be given a
79 // pointer to a freshly-allocated or copied (by a deep-copying constructor) object. ALL OTHER CTORS (for classes
80 // derived from Branch) do a DEEP COPY of their RANGE-OBJECT POINTER-valued argument(s). They usually do so by
81 // allocating (using CreateRange) or copying (using a deep-copying copy constructor) and passing the result
82 // to THIS constructor.
83
Branch(Range * newrangeptr)84 Branch::Branch(Range * newrangeptr)
85 : boost::enable_shared_from_this<Branch>(),
86 m_parents(NELEM, Branch::NONBRANCH),
87 m_children(NELEM, Branch::NONBRANCH),
88 m_ID(),
89 m_equivBranch(Branch::NONBRANCH),
90 m_partitions(registry.GetDataPack().GetNPartitionForces(), FLAGLONG),
91 // Initialize m_rangePtr to pointer to original object (having been allocated by caller), not to pointer
92 // to COPY of pointee. Ie, this is a shallow copy. All callers of this constructor must freshly-allocate
93 // or deep-copy the Range object whose pointer they pass to this function.
94 m_rangePtr(newrangeptr)
95 {
96 m_eventTime = 0.0;
97 m_updateDL = false;
98 m_marked = false;
99 m_wasCoalCalced = false;
100 m_isSample = 1;
101
102 } // Branch constructor
103
104 //------------------------------------------------------------------------------------
105 // Must delete the Range object held by pointer and always allocated on the heap.
106 // All classes derived from Branch have default (ie, empty) destructors, but the
107 // runtime systems calls this one for the base class (since all Branch classes
108 // derive from Branch), and this single destructor deallocates the Range object.
109
~Branch()110 Branch::~Branch()
111 {
112 // We need to delete the contained Range member, since we own it.
113 delete m_rangePtr;
114
115 } // Branch destructor
116
117 //------------------------------------------------------------------------------------
118 // DEEP-COPYING copy constructor.
119
Branch(const Branch & src)120 Branch::Branch(const Branch & src)
121 : boost::enable_shared_from_this<Branch>(),
122 m_parents(NELEM, Branch::NONBRANCH),
123 m_children(NELEM, Branch::NONBRANCH),
124 m_equivBranch(Branch::NONBRANCH)
125 {
126 CopyAllMembers(src); // Does a deep copy of objects (namely Range) held by pointer.
127 } // Branch copy constructor
128
129 //------------------------------------------------------------------------------------
130
CopyAllMembers(const Branch & src)131 void Branch::CopyAllMembers(const Branch & src)
132 {
133 fill(m_parents.begin(), m_parents.end(), Branch::NONBRANCH);
134 fill(m_children.begin(), m_children.end(), Branch::NONBRANCH);
135
136 m_ID = src.m_ID;
137 m_eventTime = src.m_eventTime;
138 m_partitions = src.m_partitions;
139 m_updateDL = src.m_updateDL;
140 m_marked = src.m_marked;
141 m_rangePtr = src.m_rangePtr->Clone(); // Here's the deep copy mentioned above.
142 m_DLcells = src.m_DLcells;
143 m_movingDLcells = src.m_movingDLcells;
144 m_isSample = src.m_isSample;
145 m_wasCoalCalced = src.m_wasCoalCalced;
146
147 } // Branch::CopyAllMembers
148
149 //------------------------------------------------------------------------------------
150
IsAMember(const Force & force,const LongVec1d & membership) const151 bool Branch::IsAMember(const Force & force, const LongVec1d & membership) const
152 {
153 return force.IsAMember(m_partitions, membership);
154 } // Branch::IsAMember
155
156 //------------------------------------------------------------------------------------
157 // Checks if two branches are functionally the same.
158 // If you want full equivalency use operator==().
159
IsEquivalentTo(const Branch_ptr twin) const160 bool Branch::IsEquivalentTo(const Branch_ptr twin) const
161 {
162 return m_ID == twin->m_ID;
163
164 #if 0 // Useful for debugging, perhaps.
165 if (Event() != twin->Event())
166 {
167 // cout << "Different events" << endl;
168 return false;
169 }
170
171 if (m_partitions != twin->m_partitions)
172 {
173 // cout << "Different partitions" << endl;
174 return false;
175 }
176
177 if (m_rangePtr->GetLiveSites_Rg() != twin->m_rangePtr->GetLiveSites_Rg())
178 {
179 // cout << "Different Active sites" << endl;
180 return false;
181 }
182
183 if (m_rangePtr->GetRecSite_Rg() != twin->m_rangePtr->GetRecSite_Rg())
184 {
185 // cout << "Different recombination site" << endl;
186 return false;
187 }
188
189 if (m_eventTime != twin->m_eventTime)
190 {
191 // cout << "Different event times" << endl;
192 return false;
193 }
194
195 return true;
196 #endif
197
198 } // IsEquivalentTo
199
200 //------------------------------------------------------------------------------------
201
HasSamePartitionsAs(const Branch_ptr twin) const202 bool Branch::HasSamePartitionsAs(const Branch_ptr twin) const
203 {
204 if (twin)
205 {
206 if (m_partitions != twin->m_partitions)
207 {
208 return false;
209 }
210 }
211
212 return true;
213
214 } //HasSamePartitionsAs
215
216 //------------------------------------------------------------------------------------
217
PartitionsConsistentWith(const Branch_ptr child) const218 bool Branch::PartitionsConsistentWith(const Branch_ptr child) const
219 {
220 assert(Event() == btypeRec);
221 if (child)
222 {
223 const ForceSummary & fs = registry.GetForceSummary();
224 if (fs.CheckForce(force_DISEASE))
225 {
226 if (m_rangePtr->IsDiseaseSiteTransmitted_Rg())
227 {
228 long diseaseIndex = fs.GetPartIndex(force_DISEASE);
229 long diseaseState = child->GetPartition(force_DISEASE);
230 if (diseaseState != m_partitions[diseaseIndex])
231 {
232 return false;
233 }
234 }
235 }
236 }
237
238 return true;
239
240 } // PartitionsConsistentWith
241
242 //------------------------------------------------------------------------------------
243
ResetBuffersForNextRearrangement()244 void Branch::ResetBuffersForNextRearrangement()
245 {
246 m_rangePtr->SetOldInfoToCurrent_Rg();
247 m_rangePtr->ClearNewTargetLinks_Rg();
248
249 } // ResetBuffersForNextRearrangement
250
251 //------------------------------------------------------------------------------------
252
GetLocalPartitions() const253 LongVec1d Branch::GetLocalPartitions() const
254 {
255 LongVec1d lppartitions;
256 LongVec1d lpindex(registry.GetForceSummary().GetLocalPartitionIndexes());
257 LongVec1d::iterator lpforce;
258
259 for(lpforce = lpindex.begin(); lpforce != lpindex.end(); ++lpforce)
260 {
261 lppartitions.push_back(m_partitions[*lpforce]);
262 }
263
264 return lppartitions;
265
266 } // GetLocalPartitions
267
268 //------------------------------------------------------------------------------------
269
GetID() const270 long Branch::GetID() const
271 {
272 return m_ID.ID();
273 }
274
275 //------------------------------------------------------------------------------------
276
GetEquivBranch() const277 weakBranch_ptr Branch::GetEquivBranch() const
278 {
279 return m_equivBranch;
280 }
281
282 //------------------------------------------------------------------------------------
283
SetEquivBranch(Branch_ptr twin)284 void Branch::SetEquivBranch(Branch_ptr twin)
285 {
286 m_equivBranch = twin;
287 }
288
289 //------------------------------------------------------------------------------------
290
GetPartition(force_type force) const291 long Branch::GetPartition(force_type force) const
292 {
293 return m_partitions[registry.GetForceSummary().GetPartIndex(force)];
294 } // Branch::GetPartition
295
296 //------------------------------------------------------------------------------------
297
SetPartition(force_type force,long val)298 void Branch::SetPartition(force_type force, long val)
299 {
300 m_partitions[registry.GetForceSummary().GetPartIndex(force)] = val;
301 }
302
303 //------------------------------------------------------------------------------------
304
CopyPartitionsFrom(Branch_ptr src)305 void Branch::CopyPartitionsFrom(Branch_ptr src)
306 {
307 m_partitions = src->m_partitions;
308 } // Branch::CopyPartitionsFrom
309
310 //------------------------------------------------------------------------------------
311
MarkParentsForDLCalc()312 void Branch::MarkParentsForDLCalc()
313 {
314 // Evil kludge necessary because we force parents/children to always have two spaces even when empty!
315 if (!m_parents[0].expired())
316 {
317 Branch_ptr parent0(m_parents[0]);
318 if (!parent0->m_updateDL)
319 {
320 parent0->SetUpdateDL();
321 parent0->MarkParentsForDLCalc();
322 }
323 }
324
325 if (!m_parents[1].expired())
326 {
327 Branch_ptr parent1(m_parents[1]);
328 if (!parent1->m_updateDL)
329 {
330 parent1->SetUpdateDL();
331 parent1->MarkParentsForDLCalc();
332 }
333 }
334
335 } // Branch::MarkParentsForDLCalc
336
337 //------------------------------------------------------------------------------------
338
ReplaceChild(Branch_ptr,Branch_ptr newchild)339 void Branch::ReplaceChild(Branch_ptr, Branch_ptr newchild)
340 {
341 // This version is used by MBranch and DBranch; there are overrides for other branches.
342 newchild->m_parents[0] = shared_from_this();
343 m_children[0] = newchild;
344
345 } // Branch::ReplaceChild
346
347 //------------------------------------------------------------------------------------
348
HasSameActive(const Branch & br)349 bool Branch::HasSameActive(const Branch & br)
350 {
351 return m_rangePtr->SameLiveSites_Rg(br.m_rangePtr->GetLiveSites_Rg());
352 } // Branch::HasSameActive
353
354 //------------------------------------------------------------------------------------
355
GetDLCell(long loc,long ind,bool moving) const356 const Cell_ptr Branch::GetDLCell(long loc, long ind, bool moving) const
357 {
358 if (moving)
359 {
360 assert(loc < static_cast<long>(m_movingDLcells.size()));
361 assert(ind < static_cast<long>(m_movingDLcells[loc].size()));
362 return m_movingDLcells[loc][ind];
363 }
364 else
365 {
366 assert(loc < static_cast<long>(m_DLcells.size()));
367 assert(ind < static_cast<long>(m_DLcells[loc].size()));
368 return m_DLcells[loc][ind];
369 }
370 }
371
372 //------------------------------------------------------------------------------------
373
GetDLCell(long loc,long ind,bool moving)374 Cell_ptr Branch::GetDLCell(long loc, long ind, bool moving)
375 {
376 if (moving)
377 {
378 assert(loc < static_cast<long>(m_movingDLcells.size()));
379 assert(ind < static_cast<long>(m_movingDLcells[loc].size()));
380 return m_movingDLcells[loc][ind];
381 }
382 else
383 {
384 assert(loc < static_cast<long>(m_DLcells.size()));
385 assert(ind < static_cast<long>(m_DLcells[loc].size()));
386 return m_DLcells[loc][ind];
387 }
388 }
389
390 //------------------------------------------------------------------------------------
391
HowFarTo(const Branch & br) const392 double Branch::HowFarTo(const Branch & br) const
393 {
394 return fabs(br.m_eventTime - m_eventTime);
395 } // Branch::HowFarTo
396
397 //------------------------------------------------------------------------------------
398
GetValidChild(Branch_ptr br,long whichpos)399 Branch_ptr Branch::GetValidChild(Branch_ptr br, long whichpos)
400 {
401 #ifdef PRINT_TREE_DETAILS
402 cout << " In GetValidChild whichpos: " << whichpos << " br: " << br << endl;
403 #endif
404
405 Branch_ptr pChild = br->GetActiveChild(whichpos);
406 if (pChild)
407 br = GetValidChild(pChild, whichpos);
408
409 #ifdef PRINT_TREE_DETAILS
410 cout << " return from GetValidChild br: " << br << endl;
411 #endif
412
413 return br;
414
415 } // Branch::GetValidChild
416
417 //------------------------------------------------------------------------------------
418
GetValidPanelChild(Branch_ptr br,long whichpos)419 Branch_ptr Branch::GetValidPanelChild(Branch_ptr br, long whichpos)
420 {
421 #ifdef PRINT_TREE_DETAILS
422 cout << " In GetValidChild whichpos: " << whichpos << " br: " << br << endl;
423 #endif
424
425 Branch_ptr pChild = br->GetActivePanelChild(whichpos);
426 if (pChild)
427 br = GetValidPanelChild(pChild, whichpos);
428
429 #ifdef PRINT_TREE_DETAILS
430 cout << " return from GetValidChild br: " << br << endl;
431 #endif
432
433 return br;
434
435 } // Branch::GetValidPanelChild
436
437 //------------------------------------------------------------------------------------
438
GetValidParent(long whichpos)439 Branch_ptr Branch::GetValidParent(long whichpos)
440 {
441 // We are looking for an ancestor of our branch at which the site of interest coalesces.
442 Branch_ptr pParent(Parent(0));
443 if (Parent(1))
444 {
445 // This branch has two parents; which one is needed?
446 if (Parent(1)->m_rangePtr->IsSiteLive_Rg(whichpos)) pParent = Parent(1);
447 }
448
449 // If we have reached the bottom of the tree, we give up--there is no parent.
450 if (pParent->Event() == btypeBase) return Branch::NONBRANCH;
451
452 // Otherwise, check if this could be the parent we need.
453 if (pParent->Event() == btypeCoal)
454 {
455 // This could be the coalescence we're looking for, but does our site actually coalesce here?
456 if (pParent->Child(0)->m_rangePtr->IsSiteLive_Rg(whichpos) &&
457 pParent->Child(1)->m_rangePtr->IsSiteLive_Rg(whichpos))
458 {
459 return pParent;
460 }
461 }
462
463 // It's not a coalescence (or the right coalescence), so keep going downward.
464 return pParent->GetValidParent(whichpos);
465
466 } // GetValidParent
467
468 //------------------------------------------------------------------------------------
469
DiffersInDLFrom(Branch_ptr branch,long locus,long marker) const470 bool Branch::DiffersInDLFrom(Branch_ptr branch, long locus, long marker) const
471 {
472 return m_DLcells[locus].DiffersFrom(branch->m_DLcells[locus], marker);
473 } // Branch::DiffersInDLFrom
474
475 //------------------------------------------------------------------------------------
476 // Debugging function.
477
CheckInvariant() const478 bool Branch::CheckInvariant() const
479 {
480 // Check correct time relationships among parents and offspring.
481 long index;
482 for (index = 0; index < NELEM; ++index)
483 {
484 // If the child exists it must be earlier.
485 if (Child(index))
486 {
487 if (Child(index)->m_eventTime > m_eventTime)
488 {
489 return false;
490 }
491 }
492
493 // If the parent exists it must be later.
494 if (Parent(index))
495 {
496 if (Parent(index)->m_eventTime < m_eventTime)
497 {
498 return false;
499 }
500 }
501 }
502
503 return true;
504
505 } // CheckInvariant
506
507 //------------------------------------------------------------------------------------
508
operator ==(const Branch & src) const509 bool Branch::operator==(const Branch & src) const
510 {
511 if (Event() != src.Event()) return false;
512 if (m_partitions != src.m_partitions) return false;
513 if (m_eventTime != src.m_eventTime) return false;
514
515 return true;
516
517 } // operator==
518
519 //------------------------------------------------------------------------------------
520 // Debugging function.
521
DLCheck(const Branch & other) const522 string Branch::DLCheck(const Branch & other) const
523 {
524 unsigned long locus;
525 string problems;
526
527 if (m_DLcells.size() != other.m_DLcells.size())
528 {
529 return string(" The DLCells are different sizes--error\n");
530 }
531
532 for (locus = 0; locus < m_DLcells.size(); ++locus)
533 {
534 unsigned long ncells = m_DLcells[locus].size();
535 if (ncells != other.m_DLcells[locus].size())
536 {
537 return string(" Bad branch comparison--error\n");
538 }
539
540 unsigned long ind;
541 for (ind = 0; ind < ncells; ++ind)
542 {
543 long badmarker = m_DLcells[locus][ind]->DiffersFrom(other.m_DLcells[locus][ind]);
544 if (badmarker == FLAGLONG) continue;
545
546 problems += " Branch " + ToString(m_ID.ID());
547 problems += " differs from other branch " + ToString(other.m_ID.ID());
548 problems += " at marker " + ToString(badmarker);
549 problems += "\n";
550
551 return problems;
552 }
553 }
554
555 return problems;
556
557 } // Branch::DLCheck
558
559 //------------------------------------------------------------------------------------
560 // Debugging function.
561
PrintInfo_Br() const562 void Branch::PrintInfo_Br() const
563 {
564 cout << "ID: " << m_ID.ID() << " )" << ToString(Event()) << ")" << endl;
565 cout << "Partitions: ";
566
567 for (unsigned long i = 0; i < m_partitions.size(); i++)
568 {
569 cout << m_partitions[i] << " ";
570 }
571
572 cout << endl;
573 cout << "Event time: " << m_eventTime << endl;
574 m_rangePtr->PrintInfo_Rg();
575 }
576
577 //------------------------------------------------------------------------------------
578 // Debugging function.
579
GetBranchChildren()580 vector<Branch_ptr> Branch::GetBranchChildren()
581 {
582 vector<Branch_ptr> newkids;
583 vector<weakBranch_ptr>::iterator kid;
584
585 for(kid = m_children.begin(); kid != m_children.end(); ++kid)
586 {
587 newkids.push_back(kid->lock());
588 }
589
590 return newkids;
591 }
592
593 //------------------------------------------------------------------------------------
594 // This function is non-const because we need boost::shared_from_this()!
595
ConnectedTo(const Branch_ptr family)596 bool Branch::ConnectedTo(const Branch_ptr family)
597 {
598 if (Child(0) && Child(0) == family) return true;
599 if (Child(1) && Child(1) == family) return true;
600 if (Parent(0) && Parent(0) == family) return true;
601 if (Parent(1) && Parent(1) == family) return true;
602
603 cout << endl << "Branch " << family->m_ID.ID();
604 cout << " is not connected to branch " << m_ID.ID() << " ";
605 cout << "even though it thinks it is via it's ";
606
607 Branch_ptr me(shared_from_this());
608 if (family->Child(0) && family->Child(0) == me) cout << "child0";
609 if (family->Child(1) && family->Child(1) == me) cout << "child1";
610 if (family->Parent(0) && family->Parent(0) == me) cout << "parent0";
611 if (family->Parent(1) && family->Parent(1) == me) cout << "parent1";
612
613 cout << endl;
614 return false;
615 }
616
617 //------------------------------------------------------------------------------------
618 // Debugging function.
619
IsSameExceptForTimes(const Branch_ptr other) const620 bool Branch::IsSameExceptForTimes(const Branch_ptr other) const
621 {
622 if (Event() != other->Event()) return false;
623 if (m_partitions != other->m_partitions) return false;
624
625 return true;
626 }
627
628 //------------------------------------------------------------------------------------
629 // Debugging function.
630
RevalidateRange(FC_Status &) const631 bool Branch::RevalidateRange(FC_Status &) const
632 {
633 return m_rangePtr->SameAsChild_Rg(Child(0)->m_rangePtr);
634 } // Branch::RevalidateRange
635
636 //------------------------------------------------------------------------------------
637
AddGraphML(TiXmlElement * elem) const638 void Branch::AddGraphML(TiXmlElement * elem) const
639 {
640 ////////////////////////////////////////////////////////////
641 // node itself
642 long myID = GetID();
643 long canonicalID = GetCanonicalID();
644 if (canonicalID == myID)
645 {
646 TiXmlElement * node = new TiXmlElement("node");
647 node->SetAttribute("id",ToString(canonicalID));
648 elem->LinkEndChild(node);
649
650 // type of node
651 TiXmlElement * typeInfo = new TiXmlElement("data");
652 node->LinkEndChild(typeInfo);
653 typeInfo->SetAttribute("key","node_type");
654 TiXmlText * nTypeText = new TiXmlText(GetGraphMLNodeType());
655 typeInfo->LinkEndChild(nTypeText);
656
657 // time of node
658 TiXmlElement * timeInfo = new TiXmlElement("data");
659 node->LinkEndChild(timeInfo);
660 timeInfo->SetAttribute("key","node_time");
661 TiXmlText * nTimeText = new TiXmlText(ToString(m_eventTime));
662 timeInfo->LinkEndChild(nTimeText);
663
664 AddNodeInfo(node);
665 }
666
667 ////////////////////////////////////////////////////////////
668 // branch(es) themselves
669 TiXmlElement * branchElem = new TiXmlElement("edge");
670 elem->LinkEndChild(branchElem);
671 branchElem->SetAttribute("source",ToString(GetCanonicalParentID()));
672 branchElem->SetAttribute("target",ToString(canonicalID));
673
674 // testing partitions
675 TiXmlElement * partElem = new TiXmlElement("data");
676 branchElem->LinkEndChild(partElem);
677 partElem->SetAttribute("key","partitions");
678 TiXmlText * pTypeText = new TiXmlText(ToString(m_partitions));
679 partElem->LinkEndChild(pTypeText);
680
681 // testing range info
682 TiXmlElement * rangeElem = new TiXmlElement("data");
683 branchElem->LinkEndChild(rangeElem);
684 rangeElem->SetAttribute("key","active_sites");
685 rangeset rset = m_rangePtr->GetLiveSites_Rg();
686 TiXmlText * rTypeText = new TiXmlText(ToString(rset));
687 rangeElem->LinkEndChild(rTypeText);
688
689 // AddBranchInfo(branchElem);
690 }
691
692 //------------------------------------------------------------------------------------
693
GetParentIDs() const694 string Branch::GetParentIDs() const
695 {
696 string outstr = "";
697 long pCount = NParents();
698 for(long count = 0 ; count < pCount; count++)
699 {
700 const Branch_ptr p = Parent(count);
701 if (p != Branch::NONBRANCH)
702 {
703 outstr = outstr + " " + ToString(p->GetID());
704 }
705 }
706 return outstr;
707 }
708
709 //------------------------------------------------------------------------------------
710
GetChildIDs() const711 string Branch::GetChildIDs() const
712 {
713 string outstr = "";
714 long pCount = NChildren();
715 for(long count = 0 ; count < pCount; count++)
716 {
717 const Branch_ptr p = Child(count);
718 if (p != Branch::NONBRANCH)
719 {
720 outstr = outstr + " " + ToString(p->GetID());
721 }
722 }
723 return outstr;
724 }
725
726 //------------------------------------------------------------------------------------
727
GetCanonicalID() const728 long Branch::GetCanonicalID() const
729 {
730 long myID = GetID();
731 const Branch_ptr p = GetRecPartner();
732 if (p != Branch::NONBRANCH)
733 {
734 long otherID = p->GetID();
735 if (otherID < myID)
736 {
737 return otherID;
738 }
739 }
740 return myID;
741 }
742
743 //------------------------------------------------------------------------------------
744
GetCanonicalParentID() const745 long Branch::GetCanonicalParentID() const
746 {
747 long downID = -1;
748 long pCount = NParents();
749 long trueCount = 0;
750 for(long count = 0 ; count < pCount; count++)
751 {
752 const Branch_ptr p = Parent(count);
753 if (p != Branch::NONBRANCH)
754 {
755 long canonicalParent = p->GetCanonicalID();
756 if (trueCount == 0)
757 {
758 downID = canonicalParent;
759 }
760 else
761 {
762 if (downID != canonicalParent)
763 {
764 assert(false);
765 return -2;
766 }
767 }
768 }
769 }
770 return downID;
771 }
772
773 //------------------------------------------------------------------------------------
774 //------------------------------------------------------------------------------------
775 // This constructor should only be called by the TimeList constructor.
776 // Passes a newly-allocated Range (by pointer) to the Branch constructor, which does a shallow copy.
777
BBranch()778 BBranch::BBranch()
779 : Branch(CreateRange())
780 {
781 // Deliberately blank.
782 } // BBranch constructor
783
784 //------------------------------------------------------------------------------------
785 // Returns a pointer to a newly-heap-allocated Range or RecRange object.
786
CreateRange() const787 Range * BBranch::CreateRange() const
788 {
789 // JNOTE: We would probably need to change the way tree construction is handled
790 // to do this properly (i.e. get rid of the stump logic; see tree constructors).
791 //
792 // So we create a fake Range instead, watching out for validators that will likely fail on this branch!
793 long int nsites(BASEBRANCH_NSITES);
794
795 if (registry.GetForceSummary().CheckForce(force_REC)) // This force includes recombination.
796 {
797 rangeset diseasesites; // There are no disease sites.
798 rangeset alllinks(MakeRangeset(0, nsites - 1)); // All links were and are targetable.
799 rangeset allsites(MakeRangeset(0, nsites)); // All sites are both live and transmitted.
800 return new RecRange(nsites, diseasesites, allsites, allsites, alllinks, alllinks, allsites, allsites);
801 }
802 else
803 {
804 return new Range(nsites);
805 }
806
807 } // BBranch::CreateRange
808
809 //------------------------------------------------------------------------------------
810 // Returns a pointer to a newly-allocated Branch, which contains a deep-copied Range
811 // (thanks to the deep-copying copy constructor this class inherits from Branch).
812
Clone() const813 Branch_ptr BBranch::Clone() const
814 {
815 return Branch_ptr(new BBranch(*this));
816 } // BBranch::Clone
817
818 //------------------------------------------------------------------------------------
819
ScoreEvent(TreeSummary &,BranchBuffer &) const820 void BBranch::ScoreEvent(TreeSummary &, BranchBuffer &) const
821 {
822 assert(false); // I don't think this should ever be called.
823 } // BBranch::ScoreEvent
824
825 //------------------------------------------------------------------------------------
826 // Third arg is ref for consistency with same function in other classes, which return data therein.
827
ScoreEvent(TreeSummary &,BranchBuffer &,long &) const828 void BBranch::ScoreEvent(TreeSummary &, BranchBuffer &, long &) const
829 {
830 assert(false); // I don't think this should ever be called.
831 } // BBranch::ScoreEvent
832
833 //------------------------------------------------------------------------------------
834 // Debugging function.
835
CheckInvariant() const836 bool BBranch::CheckInvariant() const
837 {
838 // Base branches have no parents...
839 long index;
840 for (index = 0; index < NELEM; ++index)
841 {
842 if (Parent(index)) return false;
843 }
844
845 // ...and one child
846 if (!Child(0)) return false;
847 if (Child(1)) return false;
848
849 if (!Branch::CheckInvariant()) return false;
850
851 return true;
852
853 } // CheckInvariant
854
855 //------------------------------------------------------------------------------------
856
GetGraphMLNodeType() const857 string BBranch::GetGraphMLNodeType() const
858 {
859 assert(false);
860 string outstr = "Error BBranch";
861 return outstr;
862 }
863
864 //------------------------------------------------------------------------------------
865 //------------------------------------------------------------------------------------
866 // Creates a new TBranch object containing a newly-allocated Range object (held by pointer),
867 // which was allocated by CreateRange and shallow-copied by the Branch constructor.
868
TBranch(const TipData & tipdata,long int nsites,const rangeset & diseasesites)869 TBranch::TBranch(const TipData & tipdata, long int nsites, const rangeset & diseasesites)
870 : Branch(CreateRange(nsites, diseasesites))
871 {
872 m_partitions = tipdata.GetBranchPartitions();
873 m_label = tipdata.label;
874
875 // This is used for SNP panel corrections.
876 if (tipdata.m_source == dsource_panel)
877 {
878 m_isSample = 0;
879 }
880
881 #if 0
882 cout <<"Tip: " << tipdata.label;
883 if (tipdata.m_source == dsource_panel)
884 cout <<" is a panel member ";
885 else
886 cout <<" is a study member ";
887 cout << endl;
888 #endif
889
890 // WARNING: must be changed for sequential-sampling case
891 m_eventTime = 0.0;
892
893 } // TBranch constructor
894
895 //------------------------------------------------------------------------------------
896 // Returns a pointer to a newly-heap-allocated Range or RecRange object.
897
CreateRange(long int nsites,const rangeset & diseasesites) const898 Range * TBranch::CreateRange(long int nsites, const rangeset & diseasesites) const
899 {
900 if (registry.GetForceSummary().CheckForce(force_REC)) // This force includes recombination.
901 {
902 rangeset alllinks(MakeRangeset(0, nsites - 1)); // All links were and are targetable.
903 rangeset allsites(MakeRangeset(0, nsites));
904 return new RecRange(nsites, diseasesites, allsites, allsites, alllinks, alllinks, allsites, allsites);
905 }
906 else
907 {
908 return new Range(nsites);
909 }
910
911 } // TBranch::CreateRange
912
913 //------------------------------------------------------------------------------------
914 // Returns a pointer to a newly-allocated Branch, which contains a deep-copied Range
915 // (thanks to the deep-copying copy constructor this class inherits from Branch).
916
Clone() const917 Branch_ptr TBranch::Clone() const
918 {
919 return Branch_ptr(new TBranch(*this));
920 } // TBranch::Clone
921
922 //------------------------------------------------------------------------------------
923 // Debugging function.
924
CheckInvariant() const925 bool TBranch::CheckInvariant() const
926 {
927 // Tip branches have no children...
928 long index;
929 for (index = 0; index < NELEM; ++index)
930 {
931 if (Child(index)) return false;
932 }
933
934 // ...and at least one parent.
935 if (!Parent(0)) return false;
936
937 if (!Branch::CheckInvariant()) return false;
938
939 return true;
940
941 } // CheckInvariant
942
943 //------------------------------------------------------------------------------------
944
operator ==(const Branch & src) const945 bool TBranch::operator==(const Branch & src) const
946 {
947 return ((Branch::operator==(src)) && (m_label == dynamic_cast<const TBranch &>(src).m_label));
948 } // operator==
949
950 //------------------------------------------------------------------------------------
951 // Debugging function.
952
IsSameExceptForTimes(const Branch_ptr src) const953 bool TBranch::IsSameExceptForTimes(const Branch_ptr src) const
954 {
955 return ((Branch::IsSameExceptForTimes(src)) && (m_label == boost::dynamic_pointer_cast<const TBranch>(src)->m_label));
956 } // IsSameExceptForTimes
957
958 //------------------------------------------------------------------------------------
959
ScoreEvent(TreeSummary &,BranchBuffer &) const960 void TBranch::ScoreEvent(TreeSummary &, BranchBuffer &) const
961 {
962 assert(false); // I don't think this should ever be called.
963 } // TBranch::ScoreEvent
964
965 //------------------------------------------------------------------------------------
966 // Third arg is ref for consistency with same function in other classes, which return data therein.
967
ScoreEvent(TreeSummary &,BranchBuffer &,long &) const968 void TBranch::ScoreEvent(TreeSummary &, BranchBuffer &, long &) const
969 {
970 assert(false); // I don't think this should ever be called.
971 } // TBranch::ScoreEvent
972
973 //------------------------------------------------------------------------------------
974 // Debugging function.
975
RevalidateRange(FC_Status &) const976 bool TBranch::RevalidateRange(FC_Status &) const
977 {
978 if (m_rangePtr->LiveSitesOnly_Rg()) return true;
979
980 long int nsites(m_rangePtr->NumTotalSites_Rg());
981 rangeset diseasesites(m_rangePtr->GetDiseaseSites_Rg()); // Rec-only?
982 // MDEBUG following is not good form nor exception-safe; consider
983 // fixing
984 Range * newrangeptr(CreateRange(nsites, diseasesites));
985 bool matches = (*newrangeptr == *m_rangePtr);
986 delete newrangeptr;
987
988 return matches;
989
990 } // TBranch::RevalidateRange
991
992 //------------------------------------------------------------------------------------
993 // Debugging function.
994
PrintInfo_Br() const995 void TBranch::PrintInfo_Br() const
996 {
997 cout << "Event type: " << ToString(Event()) << endl;
998 cout << "Label: " << m_label << endl;
999 cout << "ID: " << m_ID.ID() << endl;
1000 cout << "Partitions: ";
1001
1002 for (unsigned long i = 0; i < m_partitions.size(); i++)
1003 {
1004 cout << m_partitions[i] << " ";
1005 }
1006
1007 cout << endl;
1008 cout << "Event time: " << m_eventTime << endl;
1009 m_rangePtr->PrintInfo_Rg();
1010
1011 } // TBranch::PrintInfo_Br
1012
1013 //------------------------------------------------------------------------------------
1014
GetGraphMLNodeType() const1015 string TBranch::GetGraphMLNodeType() const
1016 {
1017 return "Tip";
1018 }
1019
1020 //------------------------------------------------------------------------------------
1021
AddNodeInfo(TiXmlElement * nodeElem) const1022 void TBranch::AddNodeInfo(TiXmlElement * nodeElem) const
1023 {
1024 // type of node
1025 TiXmlElement * labelInfo = new TiXmlElement("data");
1026 nodeElem->LinkEndChild(labelInfo);
1027 labelInfo->SetAttribute("key","node_label");
1028 TiXmlText * labelText = new TiXmlText(m_label);
1029 labelInfo->LinkEndChild(labelText);
1030 }
1031
1032 //------------------------------------------------------------------------------------
1033 //------------------------------------------------------------------------------------
1034 // Creates a new CBranch object containing a newly-allocated Range object (held by pointer),
1035 // which was allocated by CreateRange and shallow-copied by the Branch constructor.
1036
CBranch(const Range * const child1rangeptr,const Range * const child2rangeptr,bool newbranchisinactive,const rangeset & fcsites)1037 CBranch::CBranch(const Range * const child1rangeptr, const Range * const child2rangeptr,
1038 bool newbranchisinactive, const rangeset & fcsites)
1039 : Branch(CreateRange(child1rangeptr, child2rangeptr, newbranchisinactive, fcsites))
1040 {
1041 // Deliberately blank.
1042 } // CBranch constructor
1043
1044 //------------------------------------------------------------------------------------
1045 // Returns a pointer to a newly-heap-allocated Range or RecRange object.
1046
CreateRange(const Range * const child1rangeptr,const Range * const child2rangeptr,bool newbranchisinactive,const rangeset & fcsites) const1047 Range * CBranch::CreateRange(const Range * const child1rangeptr, const Range * const child2rangeptr,
1048 bool newbranchisinactive, const rangeset & fcsites) const
1049 {
1050 // These we can grab from either child as they should be identical in these fields.
1051 long int nsites(child1rangeptr->NumTotalSites_Rg());
1052
1053 if (registry.GetForceSummary().CheckForce(force_REC)) // This force includes recombination.
1054 {
1055 rangeset diseasesites(child1rangeptr->GetDiseaseSites_Rg());
1056 rangeset transmittedsites(MakeRangeset(0, nsites)); // All sites are transmitted.
1057 rangeset livesites(Union(child1rangeptr->GetLiveSites_Rg(), child2rangeptr->GetLiveSites_Rg()));
1058 rangeset targetsites(RemoveRangeFromRange(fcsites, livesites));
1059 targetsites = Union(diseasesites, targetsites);
1060 rangeset targetlinks(TargetLinksFrom(targetsites));
1061
1062 rangeset oldtargetlinks, oldtargetsites, oldlivesites;
1063 if (newbranchisinactive) // Copy the inactive child oldtargetlinks and oldtargetsites.
1064 {
1065 oldtargetlinks = child1rangeptr->GetOldTargetLinks_Rg();
1066 oldtargetsites = child1rangeptr->GetOldTargetSites_Rg();
1067 oldlivesites = child1rangeptr->GetOldLiveSites_Rg();
1068 }
1069 else
1070 {
1071 oldtargetlinks = targetlinks;
1072 oldtargetsites = targetsites;
1073 oldlivesites = livesites;
1074 }
1075
1076 return new RecRange(nsites, diseasesites, transmittedsites, livesites,
1077 targetlinks, oldtargetlinks, oldtargetsites, oldlivesites);
1078 }
1079 else
1080 {
1081 return new Range(nsites);
1082 }
1083
1084 } // CBranch::CreateRange
1085
1086 //------------------------------------------------------------------------------------
1087 // Returns a pointer to a newly-allocated Branch, which contains a deep-copied Range
1088 // (thanks to the deep-copying copy constructor this class inherits from Branch).
1089
Clone() const1090 Branch_ptr CBranch::Clone() const
1091 {
1092 return Branch_ptr(new CBranch(*this));
1093 } // CBranch::Clone
1094
1095 //------------------------------------------------------------------------------------
1096
CanRemove(Branch_ptr checkchild)1097 bool CBranch::CanRemove(Branch_ptr checkchild)
1098 {
1099 if (m_marked) return true;
1100
1101 m_marked = true;
1102 if (Child(0) == checkchild) SetChild(0, Child(1));
1103 SetChild(1, Branch::NONBRANCH);
1104
1105 return false;
1106
1107 } // CBranch::CanRemove
1108
1109 //------------------------------------------------------------------------------------
1110
UpdateBranchRange(const rangeset & fcsites,bool dofc)1111 void CBranch::UpdateBranchRange(const rangeset & fcsites, bool dofc)
1112 {
1113 if (m_marked)
1114 {
1115 m_rangePtr->UpdateOneLeggedCRange_Rg(Child(0)->m_rangePtr);
1116 }
1117 else
1118 {
1119 m_rangePtr->UpdateCRange_Rg(Child(0)->m_rangePtr, Child(1)->m_rangePtr, fcsites, dofc);
1120 }
1121
1122 } // CBranch::UpdateBranchRange
1123
1124 //------------------------------------------------------------------------------------
1125
UpdateRootBranchRange(const rangeset & fcsites,bool dofc)1126 void CBranch::UpdateRootBranchRange(const rangeset & fcsites, bool dofc)
1127 {
1128 if (m_marked)
1129 {
1130 m_rangePtr->UpdateOneLeggedRootRange_Rg(Child(0)->m_rangePtr);
1131 }
1132 else
1133 {
1134 m_rangePtr->UpdateRootRange_Rg(Child(0)->m_rangePtr, Child(1)->m_rangePtr, fcsites, dofc);
1135 }
1136
1137 } // CBranch::UpdateRootBranchRange
1138
1139 //------------------------------------------------------------------------------------
1140
ReplaceChild(Branch_ptr pOldChild,Branch_ptr pNewChild)1141 void CBranch::ReplaceChild(Branch_ptr pOldChild, Branch_ptr pNewChild)
1142 {
1143 if (Child(0) == pOldChild)
1144 SetChild(0, pNewChild);
1145 else
1146 SetChild(1, pNewChild);
1147
1148 pNewChild->SetParent(0, shared_from_this());
1149
1150 } // CBranch::ReplaceChild
1151
1152 //------------------------------------------------------------------------------------
1153
OtherChild(Branch_ptr badchild)1154 Branch_ptr CBranch::OtherChild(Branch_ptr badchild)
1155 {
1156 if (Child(0) == badchild) return Child(1);
1157 return Child(0);
1158
1159 } // CBranch::OtherChild
1160
1161 //------------------------------------------------------------------------------------
1162 // Both children must be both active and included in the DLcalc for this to return true.
1163
CanCalcDL(long site) const1164 bool CBranch::CanCalcDL(long site) const
1165 {
1166 if ((Child(0)->m_rangePtr->IsSiteLive_Rg(site)) && (Child(1)->m_rangePtr->IsSiteLive_Rg(site)))
1167 {
1168 return true;
1169 }
1170 else
1171 {
1172 return false;
1173 }
1174 }
1175
1176 //------------------------------------------------------------------------------------
1177 // If both children are active return Branch::NONBRANCH to signal a stop; otherwise return the active child.
1178
GetActiveChild(long site) const1179 const Branch_ptr CBranch::GetActiveChild(long site) const
1180 {
1181 if (Child(0)->m_rangePtr->IsSiteLive_Rg(site))
1182 {
1183 if (Child(1)->m_rangePtr->IsSiteLive_Rg(site)) return Branch::NONBRANCH;
1184 return Child(0);
1185 }
1186 else
1187 {
1188 return Child(1);
1189 }
1190
1191 } // CBranch::GetActiveChild
1192
1193 //------------------------------------------------------------------------------------
1194
ScoreEvent(TreeSummary & summary,BranchBuffer & ks) const1195 void CBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks) const
1196 {
1197 long i;
1198 Summary * csum = summary.GetSummary(force_COAL);
1199
1200 // Forces have become inconsistent!
1201 assert(dynamic_cast<CoalSummary *>(csum));
1202
1203 long myxpart = registry.GetDataPack().GetCrossPartitionIndex(m_partitions);
1204 LongVec1d emptyvec;
1205
1206 // Score the event.
1207 csum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
1208 FLAGLONG, myxpart, FLAGLONG, FLAGLONG, emptyvec, force_COAL);
1209
1210 // Adjust the branch counts.
1211 for(i = 0; i < NELEM; ++i)
1212 ks.UpdateBranchCounts(Child(i)->m_partitions, false);
1213
1214 ks.UpdateBranchCounts(m_partitions);
1215
1216 } // CBranch::ScoreEvent
1217
1218 //------------------------------------------------------------------------------------
1219 // Third arg is a reference because this function returns results therein.
1220
ScoreEvent(TreeSummary & summary,BranchBuffer & ks,long & s) const1221 void CBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks, long & s) const
1222 {
1223 long i;
1224 Summary * csum = summary.GetSummary(force_COAL);
1225
1226 // Forces have become inconsistent!
1227 assert(dynamic_cast<CoalSummary *>(csum));
1228
1229 long myxpart = registry.GetDataPack().
1230 GetCrossPartitionIndex(m_partitions);
1231 LongVec1d emptyvec;
1232
1233 // Score the event.
1234 csum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
1235 s, myxpart, FLAGLONG, FLAGLONG, emptyvec, force_COAL);
1236
1237 // Rest for recominant branches only? - BOBGIAN
1238 // Adjust the branch and link counts.
1239 for (i = 0; i < NELEM; ++i)
1240 {
1241 ks.UpdateBranchCounts(Child(i)->m_partitions, false);
1242 s -= Child(i)->m_rangePtr->NumTargetLinks_Rg();
1243 }
1244
1245 ks.UpdateBranchCounts(m_partitions);
1246 s += m_rangePtr->NumTargetLinks_Rg();
1247
1248 } // CBranch::ScoreEvent
1249
1250 //------------------------------------------------------------------------------------
1251 // Debugging function.
1252
CheckInvariant() const1253 bool CBranch::CheckInvariant() const
1254 {
1255 // Coalescent branches have two children...
1256 if (!Child(0)) return false;
1257 if (!Child(1)) return false;
1258
1259 //...and at least one parent.
1260 if (!Parent(0)) return false;
1261
1262 if (!Branch::CheckInvariant()) return false;
1263
1264 return true;
1265
1266 } // CheckInvariant
1267
1268 //------------------------------------------------------------------------------------
1269 // Debugging function.
1270
RevalidateRange(FC_Status & fcstatus) const1271 bool CBranch::RevalidateRange(FC_Status & fcstatus) const
1272 {
1273 rangeset livesites;
1274
1275 // Are we the same moeity and population as our children?
1276 if (Child(0)->m_partitions != m_partitions)
1277 {
1278 cout << "Branch changed color in CBranch" << endl;
1279 return false;
1280 }
1281
1282 if (Child(1) && Child(1)->m_partitions != m_partitions)
1283 {
1284 cout << "Branch changed color in CBranch" << endl;
1285 return false;
1286 }
1287
1288 if (!Child(1)) // We're in a one-legger.
1289 {
1290 const Range * const childrangeptr(Child(0)->m_rangePtr);
1291 livesites = childrangeptr->GetLiveSites_Rg();
1292 if (!(m_rangePtr->SameLiveSites_Rg(livesites)))
1293 {
1294 return false;
1295 }
1296 // We don't need to update fc stuff for a one-legger.
1297 }
1298 else
1299 {
1300 // We're in a normal coalescence.
1301 const Range * const child0rangeptr(Child(0)->m_rangePtr);
1302 const Range * const child1rangeptr(Child(1)->m_rangePtr);
1303
1304 livesites = Union(child0rangeptr->GetLiveSites_Rg(), child1rangeptr->GetLiveSites_Rg());
1305 if (!(m_rangePtr->SameLiveSites_Rg(livesites)))
1306 {
1307 return false;
1308 }
1309
1310 #if LAMARC_QA_FC_ON
1311 rangeset fcsites(Intersection(child0rangeptr->GetLiveSites_Rg(), child1rangeptr->GetLiveSites_Rg()));
1312 fcstatus.Decrement_FC_Counts(fcsites);
1313 #endif
1314 }
1315
1316 if (m_rangePtr->LiveSitesOnly_Rg()) return true; // We're done checking.
1317
1318 // Will rest of this function only execute on recombinant branches?
1319 #if LAMARC_QA_FC_ON
1320 rangeset target(TargetLinksFrom(Union(m_rangePtr->GetDiseaseSites_Rg(),
1321 RemoveRangeFromRange(fcstatus.Coalesced_Sites(), livesites))));
1322 #else
1323 rangeset target(TargetLinksFrom(Union(m_rangePtr->GetDiseaseSites_Rg(), livesites)));
1324 #endif
1325
1326 if (m_rangePtr->DifferentTargetLinks_Rg(target))
1327 {
1328 cout << "my range; ";
1329 m_rangePtr->PrintInfo_Rg();
1330 cout << "target: " << ToString(target) << endl;
1331 return false;
1332 }
1333
1334 // Is this a rec-only branch? Otherwise, we may not want GetOldTargetLinks_Rg() to ASSERT.
1335 rangeset newtarget(RemoveRangeFromRange(m_rangePtr->GetOldTargetLinks_Rg(), target));
1336 if (m_rangePtr->DifferentNewTargetLinks_Rg(newtarget))
1337 {
1338 return false;
1339 }
1340
1341 return true;
1342
1343 } // CBranch::RevalidateRange
1344
1345 //------------------------------------------------------------------------------------
1346
GetGraphMLNodeType() const1347 string CBranch::GetGraphMLNodeType() const
1348 {
1349 return "Coal";
1350 }
1351
1352 //------------------------------------------------------------------------------------
1353 //------------------------------------------------------------------------------------
1354 // Creates a new PartitionBranch object containing a newly-allocated Range object
1355 // (held by pointer), which was allocated by CreateRange and shallow-copied by the Branch constructor.
1356
PartitionBranch(const Range * const childrangeptr)1357 PartitionBranch::PartitionBranch(const Range * const childrangeptr)
1358 : Branch(CreateRange(childrangeptr))
1359 {
1360 // Deliberately blank.
1361 } // PartitionBranch constructor
1362
1363 //------------------------------------------------------------------------------------
1364 // Returns a pointer to a newly-heap-allocated Range or RecRange object.
1365
CreateRange(const Range * const childrangeptr) const1366 Range * PartitionBranch::CreateRange(const Range * const childrangeptr) const
1367 {
1368 // NOT copy-constructed from the child, because a PartitionBranch always transmits all sites, and the child might not.
1369 long int nsites(childrangeptr->NumTotalSites_Rg());
1370
1371 if (registry.GetForceSummary().CheckForce(force_REC)) // This force includes recombination.
1372 {
1373 rangeset diseasesites(childrangeptr->GetDiseaseSites_Rg());
1374 rangeset transmittedsites(MakeRangeset(0, nsites));
1375 rangeset livesites(childrangeptr->GetLiveSites_Rg());
1376 rangeset targetlinks(childrangeptr->GetTargetLinks_Rg());
1377 rangeset oldtargetlinks(childrangeptr->GetOldTargetLinks_Rg());
1378 rangeset oldtargetsites(childrangeptr->GetOldTargetSites_Rg());
1379 rangeset oldlivesites(childrangeptr->GetOldLiveSites_Rg());
1380
1381 return new RecRange(nsites, diseasesites, transmittedsites, livesites,
1382 targetlinks, oldtargetlinks, oldtargetsites, oldlivesites);
1383 }
1384 else
1385 {
1386 return new Range(nsites);
1387 }
1388
1389 } // PartitionBranch::CreateRange
1390
1391 //------------------------------------------------------------------------------------
1392
UpdateBranchRange(const rangeset &,bool)1393 void PartitionBranch::UpdateBranchRange(const rangeset &, bool)
1394 {
1395 m_rangePtr->UpdateMRange_Rg(Child(0)->m_rangePtr);
1396 } // PartitionBranch::UpdateBranchRange
1397
1398 //------------------------------------------------------------------------------------
1399 // Debugging function.
1400
CheckInvariant() const1401 bool PartitionBranch::CheckInvariant() const
1402 {
1403 // Partition branches have one child...
1404 if (!Child(0)) return false;
1405 if (Child(1)) return false;
1406
1407 //...and at least one parent.
1408 if (!Parent(0)) return false;
1409
1410 if (!Branch::CheckInvariant()) return false;
1411
1412 return true;
1413
1414 } // CheckInvariant
1415
1416 //------------------------------------------------------------------------------------
1417 //------------------------------------------------------------------------------------
1418 // Creates a new MigLikeBranch object containing a newly-allocated Range object (held by pointer),
1419 // which was allocated by CreateRange in the PartitionBranch constructor.
1420
MigLikeBranch(const Range * const protorangeptr)1421 MigLikeBranch::MigLikeBranch(const Range * const protorangeptr)
1422 : PartitionBranch(protorangeptr)
1423 {
1424 // Deliberately blank.
1425 } // MigLikeBranch constructor
1426
1427 //------------------------------------------------------------------------------------
1428 // Creates a new MigLikeBranch object containing a newly-allocated Range object (held by pointer), which
1429 // was allocated by CreateRange and shallow-copied by the Branch constructor that PartitionBranch inherits.
1430
MigLikeBranch(const MigLikeBranch & src)1431 MigLikeBranch::MigLikeBranch(const MigLikeBranch & src)
1432 : PartitionBranch(src)
1433 {
1434 // Deliberately blank.
1435 } // MigLikeBranch copy constructor
1436
1437 //------------------------------------------------------------------------------------
1438 //------------------------------------------------------------------------------------
1439 // Creates a new MBranch object containing a newly-allocated Range object (held by pointer), which
1440 // was allocated by CreateRange and shallow-copied by the Branch constructor that MigLikeBranch inherits.
1441
MBranch(const Range * const protorangeptr)1442 MBranch::MBranch(const Range * const protorangeptr)
1443 : MigLikeBranch(protorangeptr)
1444 {
1445 // Deliberately blank.
1446 } // MBranch constructor
1447
1448 //------------------------------------------------------------------------------------
1449 // Creates a new MBranch object containing a newly-allocated Range object (held by pointer), which
1450 // was allocated by CreateRange and shallow-copied by the Branch constructor that MigLikeBranch inherits.
1451
MBranch(const MBranch & src)1452 MBranch::MBranch(const MBranch & src)
1453 : MigLikeBranch(src)
1454 {
1455 // Deliberately blank.
1456 } // MBranch copy constructor
1457
1458 //------------------------------------------------------------------------------------
1459 // Creates a new MBranch object containing a newly-allocated Range object (held by pointer), which
1460 // was deep-copied by the PartitionBranch copy constructor that MigLikeBranch inherits.
1461
Clone() const1462 Branch_ptr MBranch::Clone() const
1463 {
1464 return Branch_ptr(new MBranch(*this));
1465 } // MBranch::Clone
1466
1467 //------------------------------------------------------------------------------------
1468
ScoreEvent(TreeSummary & summary,BranchBuffer & ks) const1469 void MBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks) const
1470 {
1471 Summary * msum = summary.GetSummary(force_MIG);
1472
1473 // MDEBUG is correct with divergence? Forces have become inconsistent!
1474 assert(dynamic_cast<MigSummary *>(msum));
1475
1476 long mypop = GetPartition(force_MIG);
1477 long chpop = Child(0)->GetPartition(force_MIG);
1478 LongVec1d emptyvec;
1479
1480 // Score the event.
1481 msum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
1482 FLAGLONG, mypop, chpop, FLAGLONG, emptyvec, force_MIG);
1483
1484 // Adjust the branch counts.
1485 assert(!Child(1)); // too many children??
1486 ks.UpdateBranchCounts(m_partitions);
1487 ks.UpdateBranchCounts(Child(0)->m_partitions, false);
1488
1489 // We do not adjust active sites; they cannot possibly change.
1490
1491 } // MBranch::ScoreEvent
1492
1493 //------------------------------------------------------------------------------------
1494 // Third arg is ref for consistency with same function in other classes, which return data therein.
1495
ScoreEvent(TreeSummary & summary,BranchBuffer & ks,long & s) const1496 void MBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks, long & s) const
1497 {
1498 Summary * msum = summary.GetSummary(force_MIG);
1499
1500 // Forces have become inconsistent!
1501 assert(dynamic_cast<MigSummary *>(msum));
1502
1503 long mypop = GetPartition(force_MIG);
1504 long chpop = Child(0)->GetPartition(force_MIG);
1505 LongVec1d emptyvec;
1506
1507 // Score the event.
1508 msum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
1509 s, mypop, chpop, FLAGLONG, emptyvec, force_MIG);
1510
1511 // Adjust the branch counts.
1512 assert(!Child(1)); // too many children??
1513 ks.UpdateBranchCounts(m_partitions);
1514 ks.UpdateBranchCounts(Child(0)->m_partitions, false);
1515
1516 // We do not adjust active sites; they cannot possibly change.
1517
1518 } // MBranch::ScoreEvent
1519
1520 //------------------------------------------------------------------------------------
1521
GetGraphMLNodeType() const1522 string MBranch::GetGraphMLNodeType() const
1523 {
1524 string outstr = "Mig";
1525 return outstr;
1526 }
1527
1528 //------------------------------------------------------------------------------------
1529 //------------------------------------------------------------------------------------
1530 // Divergence migration branches. Practically identical to MBranch.
1531
DivMigBranch(const Range * const protorangeptr)1532 DivMigBranch::DivMigBranch(const Range * const protorangeptr)
1533 : MigLikeBranch(protorangeptr)
1534 {
1535 // Deliberately blank.
1536 } // DivMigBranch constructor
1537
1538 //------------------------------------------------------------------------------------
1539
DivMigBranch(const DivMigBranch & src)1540 DivMigBranch::DivMigBranch(const DivMigBranch & src)
1541 : MigLikeBranch(src)
1542 {
1543 // Deliberately blank.
1544 } // DivMigBranch copy constructor
1545
1546 //------------------------------------------------------------------------------------
1547
Clone() const1548 Branch_ptr DivMigBranch::Clone() const
1549 {
1550 return Branch_ptr(new DivMigBranch(*this));
1551 } // DivMigBranch::Clone
1552
1553 //------------------------------------------------------------------------------------
1554
ScoreEvent(TreeSummary & summary,BranchBuffer & ks) const1555 void DivMigBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks) const
1556 {
1557 Summary * msum = summary.GetSummary(force_DIVMIG);
1558
1559 // MDEBUG is correct with divergence? Forces have become inconsistent!
1560 assert(dynamic_cast<DivMigSummary *>(msum));
1561
1562 long mypop = GetPartition(force_DIVMIG);
1563 long chpop = Child(0)->GetPartition(force_DIVMIG);
1564 LongVec1d emptyvec;
1565
1566 // Score the event.
1567 msum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
1568 FLAGLONG, mypop, chpop, FLAGLONG, emptyvec, force_DIVMIG);
1569
1570 // Adjust the branch counts.
1571 assert(!Child(1)); // too many children??
1572 ks.UpdateBranchCounts(m_partitions);
1573 ks.UpdateBranchCounts(Child(0)->m_partitions, false);
1574
1575 // We do not adjust active sites; they cannot possibly change.
1576
1577 } // DivMigBranch::ScoreEvent
1578
1579 //------------------------------------------------------------------------------------
1580 // Third arg is ref for consistency with same function in other classes, which return data therein.
1581
ScoreEvent(TreeSummary & summary,BranchBuffer & ks,long & s) const1582 void DivMigBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks, long & s) const
1583 {
1584 Summary * msum = summary.GetSummary(force_DIVMIG);
1585
1586 // Forces have become inconsistent!
1587 assert(dynamic_cast<DivMigSummary *>(msum));
1588
1589 long mypop = GetPartition(force_DIVMIG);
1590 long chpop = Child(0)->GetPartition(force_DIVMIG);
1591 LongVec1d emptyvec;
1592
1593 // Score the event.
1594 msum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
1595 s, mypop, chpop, FLAGLONG, emptyvec, force_DIVMIG);
1596
1597 // Adjust the branch counts.
1598 assert(!Child(1)); // too many children??
1599 ks.UpdateBranchCounts(m_partitions);
1600 ks.UpdateBranchCounts(Child(0)->m_partitions, false);
1601
1602 // We do not adjust active sites; they cannot possibly change.
1603
1604 } // DivMigBranch::ScoreEvent
1605
1606 //------------------------------------------------------------------------------------
1607
GetGraphMLNodeType() const1608 string DivMigBranch::GetGraphMLNodeType() const
1609 {
1610 string outstr = "DivMig";
1611 return outstr;
1612 }
1613
1614 //------------------------------------------------------------------------------------
1615 //------------------------------------------------------------------------------------
1616 // Disease mutation branches. NOT related to Divergence!
1617
DBranch(const Range * const protorangeptr)1618 DBranch::DBranch(const Range * const protorangeptr)
1619 : PartitionBranch(protorangeptr)
1620 {
1621 // Deliberately blank.
1622 } // DBranch constructor
1623
1624 //------------------------------------------------------------------------------------
1625
DBranch(const DBranch & src)1626 DBranch::DBranch(const DBranch & src)
1627 : PartitionBranch(src)
1628 {
1629 // Deliberately blank.
1630 } // DBranch copy constructor
1631
1632 //------------------------------------------------------------------------------------
1633
Clone() const1634 Branch_ptr DBranch::Clone() const
1635 {
1636 return Branch_ptr(new DBranch(*this));
1637 } // DBranch::Clone
1638
1639 //------------------------------------------------------------------------------------
1640
ScoreEvent(TreeSummary & summary,BranchBuffer & ks) const1641 void DBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks) const
1642 {
1643 Summary * dissum = summary.GetSummary(force_DISEASE);
1644
1645 // Forces have become inconsistent!
1646 assert(dynamic_cast<DiseaseSummary *>(dissum));
1647
1648 long mydis = GetPartition(force_DISEASE);
1649 long chdis = Child(0)->GetPartition(force_DISEASE);
1650 LongVec1d emptyvec;
1651
1652 // Score the event.
1653 dissum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
1654 FLAGLONG, mydis, chdis, FLAGLONG, emptyvec, force_DISEASE);
1655
1656 // Adjust the branch counts.
1657 assert(!Child(1)); // too many children??
1658 ks.UpdateBranchCounts(m_partitions);
1659 ks.UpdateBranchCounts(Child(0)->m_partitions, false);
1660
1661 // We do not adjust active sites; they cannot possibly change.
1662
1663 } // DBranch::ScoreEvent
1664
1665 //------------------------------------------------------------------------------------
1666 // Third arg is ref for consistency with same function in other classes, which return data therein.
1667
ScoreEvent(TreeSummary & summary,BranchBuffer & ks,long & s) const1668 void DBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks, long & s) const
1669 {
1670 Summary * dissum = summary.GetSummary(force_DISEASE);
1671
1672 // Forces have become inconsistent!
1673 assert(dynamic_cast<DiseaseSummary *>(dissum));
1674
1675 long mydis = GetPartition(force_DISEASE);
1676 long chdis = Child(0)->GetPartition(force_DISEASE);
1677 LongVec1d emptyvec;
1678
1679 // Score the event.
1680 dissum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
1681 s, mydis, chdis, FLAGLONG, emptyvec, force_DISEASE);
1682
1683 // Adjust the branch counts.
1684 assert(!Child(1)); // too many children??
1685 ks.UpdateBranchCounts(m_partitions);
1686 ks.UpdateBranchCounts(Child(0)->m_partitions, false);
1687
1688 // We do not adjust active sites; they cannot possibly change.
1689
1690 } // DBranch::ScoreEvent
1691
1692 //------------------------------------------------------------------------------------
1693
GetGraphMLNodeType() const1694 string DBranch::GetGraphMLNodeType() const
1695 {
1696 string outstr = "Disease";
1697 return outstr;
1698 }
1699
1700 //------------------------------------------------------------------------------------
1701 //------------------------------------------------------------------------------------
1702
Clone() const1703 Branch_ptr EBranch::Clone() const
1704 {
1705 return Branch_ptr(new EBranch(*this));
1706 } // EBranch::Clone
1707
1708 //------------------------------------------------------------------------------------
1709
ScoreEvent(TreeSummary & summary,BranchBuffer & ks) const1710 void EBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks) const
1711 {
1712 Summary * esum = summary.GetSummary(force_DIVERGENCE);
1713
1714 // Forces have become inconsistent!
1715 assert(dynamic_cast<EpochSummary *>(esum));
1716
1717 // NB While this is not a partition force, it cooperates with DIVMIG here.
1718 long mypop = GetPartition(force_DIVMIG);
1719 long chpop = Child(0)->GetPartition(force_DIVMIG);
1720 LongVec1d emptyvec;
1721
1722 // Score the event.
1723 esum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
1724 FLAGLONG, mypop, chpop, FLAGLONG, emptyvec, force_DIVERGENCE);
1725
1726 // Adjust the branch counts.
1727 assert(!Child(1)); // too many children??
1728 ks.UpdateBranchCounts(m_partitions);
1729 ks.UpdateBranchCounts(Child(0)->m_partitions, false);
1730
1731 // We do not adjust active sites; they cannot possibly change.
1732
1733 } // EBranch::ScoreEvent
1734
1735 //------------------------------------------------------------------------------------
1736 // Third arg is ref for consistency with same function in other classes, which return data therein.
1737
ScoreEvent(TreeSummary & summary,BranchBuffer & ks,long & s) const1738 void EBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks, long & s) const
1739 {
1740 Summary * esum = summary.GetSummary(force_DIVERGENCE);
1741
1742 // Forces have become inconsistent!
1743 assert(dynamic_cast<EpochSummary *>(esum));
1744
1745 // NB While this is not a partition force, it cooperates with DIVMIG here.
1746 long mypop = GetPartition(force_DIVMIG);
1747 long chpop = Child(0)->GetPartition(force_DIVMIG);
1748 LongVec1d emptyvec;
1749
1750 // Score the event.
1751 esum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(),
1752 s, mypop, chpop, FLAGLONG, emptyvec, force_DIVERGENCE);
1753
1754 // Adjust the branch counts.
1755 assert(!Child(1)); // too many children??
1756 ks.UpdateBranchCounts(m_partitions);
1757 ks.UpdateBranchCounts(Child(0)->m_partitions, false);
1758
1759 // We do not adjust active sites; they cannot possibly change.
1760
1761 } // EBranch::ScoreEvent
1762
1763 //------------------------------------------------------------------------------------
1764
GetGraphMLNodeType() const1765 string EBranch::GetGraphMLNodeType() const
1766 {
1767 string outstr = "Epoch";
1768 return outstr;
1769 }
1770
1771 //------------------------------------------------------------------------------------
1772 //------------------------------------------------------------------------------------
1773 // Creates a new RBranch object containing a newly-allocated Range object (held by pointer),
1774 // which was allocated by CreateRange and shallow-copied by the Branch constructor.
1775
RBranch(const Range * const childrangeptr,bool newbranchisinactive,const rangeset & transmittedsites,const rangeset & fcsites)1776 RBranch::RBranch(const Range * const childrangeptr, bool newbranchisinactive,
1777 const rangeset & transmittedsites, const rangeset & fcsites)
1778 : Branch(CreateRange(childrangeptr, newbranchisinactive, transmittedsites, fcsites))
1779 {
1780 // Deliberately blank.
1781 } // RBranch constructor
1782
1783 //------------------------------------------------------------------------------------
1784
Clone() const1785 Branch_ptr RBranch::Clone() const
1786 {
1787 return Branch_ptr(new RBranch(*this));
1788 } // RBranch::Clone
1789
1790 //------------------------------------------------------------------------------------
1791 // Returns a pointer to a newly-heap-allocated RecRange object.
1792 // ASSERTs if attempt is made to create a Range object.
1793
CreateRange(const Range * const childrangeptr,bool newbranchisinactive,const rangeset & transmittedsites,const rangeset & fcsites) const1794 Range * RBranch::CreateRange(const Range * const childrangeptr, bool newbranchisinactive,
1795 const rangeset & transmittedsites, const rangeset & fcsites) const
1796 {
1797 if (registry.GetForceSummary().CheckForce(force_REC)) // This force includes recombination.
1798 {
1799 // The following are identical with our child.
1800 long int nsites(childrangeptr->NumTotalSites_Rg());
1801 rangeset diseasesites(childrangeptr->GetDiseaseSites_Rg());
1802 rangeset livesites(Intersection(childrangeptr->GetLiveSites_Rg(), transmittedsites));
1803
1804 // MNEWCODE I am still not positive the following lines are correct!
1805 // But I can't find the case that would prove them wrong.
1806 rangeset targetsites(RemoveRangeFromRange(fcsites, livesites));
1807 targetsites = Union(diseasesites, targetsites);
1808 rangeset targetlinks(TargetLinksFrom(targetsites));
1809
1810 rangeset oldtargetlinks, oldtargetsites, oldlivesites;
1811 if (newbranchisinactive) // Copy the inactive child oldtargetlinks and oldtargetsites.
1812 {
1813 oldtargetlinks = childrangeptr->GetOldTargetLinks_Rg();
1814 oldtargetsites = childrangeptr->GetOldTargetSites_Rg();
1815 oldlivesites = childrangeptr->GetOldLiveSites_Rg();
1816 }
1817 else
1818 {
1819 oldtargetlinks = targetlinks;
1820 oldtargetsites = targetsites;
1821 oldlivesites = livesites;
1822 }
1823
1824 return new RecRange(nsites, diseasesites, transmittedsites, livesites,
1825 targetlinks, oldtargetlinks, oldtargetsites, oldlivesites);
1826 }
1827 else
1828 {
1829 assert(false); // An RBranch in a non-rec run??
1830 return NULL; // To silence compiler warning.
1831 }
1832
1833 } // RBranch::CreateRange
1834
1835 //------------------------------------------------------------------------------------
1836
CopyPartitionsFrom(Branch_ptr src)1837 void RBranch::CopyPartitionsFrom(Branch_ptr src)
1838 {
1839 m_partitions = src->m_partitions;
1840 } // Branch::CopyPartitionsFrom
1841
1842 //------------------------------------------------------------------------------------
1843
RecCopyPartitionsFrom(Branch_ptr src,FPartMap fparts,bool islow)1844 void RBranch::RecCopyPartitionsFrom(Branch_ptr src, FPartMap fparts, bool islow)
1845 {
1846 m_partitions = src->m_partitions; //fparts may override.
1847
1848 ForceVec forces(registry.GetForceSummary().GetPartitionForces());
1849
1850 for(unsigned long force = 0; force < forces.size(); force++)
1851 {
1852 FPartMapiter thisforcepart = fparts.find(forces[force]->GetTag());
1853 if (thisforcepart != fparts.end())
1854 {
1855 long chosenpart = thisforcepart->second;
1856 long partition = dynamic_cast<PartitionForce *>(forces[force])->
1857 ChoosePartition(src->m_partitions[force], chosenpart, islow, GetRecSite_Br());
1858 m_partitions[force] = partition;
1859 }
1860 }
1861
1862 } // RBranch::RecCopyPartitionsFrom
1863
1864 //------------------------------------------------------------------------------------
1865
UpdateBranchRange(const rangeset & fcsites,bool dofc)1866 void RBranch::UpdateBranchRange(const rangeset & fcsites, bool dofc)
1867 {
1868 m_rangePtr->UpdateRRange_Rg(Child(0)->m_rangePtr, fcsites, dofc);
1869 } // RBranch::UpdateBranchRange
1870
1871 //------------------------------------------------------------------------------------
1872
ReplaceChild(Branch_ptr oldchild,Branch_ptr newchild)1873 void RBranch::ReplaceChild(Branch_ptr oldchild, Branch_ptr newchild)
1874 {
1875 SetChild(0, newchild);
1876
1877 Branch_ptr myself = shared_from_this();
1878
1879 if (oldchild->Parent(0) == myself)
1880 {
1881 newchild->SetParent(0, myself);
1882 }
1883 else
1884 {
1885 newchild->SetParent(1, myself);
1886 }
1887
1888 } // RBranch::ReplaceChild
1889
1890 //------------------------------------------------------------------------------------
1891
IsRemovableRecombinationLeg(const rangeset & fcsites) const1892 bool RBranch::IsRemovableRecombinationLeg(const rangeset & fcsites) const
1893 {
1894 return (!m_rangePtr->AreChildTargetSitesTransmitted_Rg(Child(0)->m_rangePtr, fcsites));
1895 } // RBranch::IsRemovableRecombinationLeg
1896
1897 //------------------------------------------------------------------------------------
1898
operator ==(const Branch & other) const1899 bool RBranch::operator==(const Branch & other) const
1900 {
1901 // Potentially troublesome call to RecRange::operator== here - BOBGIAN.
1902 return (Branch::operator==(other) && (*m_rangePtr == *(other.m_rangePtr)));
1903 }
1904
1905 //------------------------------------------------------------------------------------
1906 // Debugging function.
1907
IsSameExceptForTimes(const Branch_ptr other) const1908 bool RBranch::IsSameExceptForTimes(const Branch_ptr other) const
1909 {
1910 // Potentially troublesome call to RecRange::operator== here - BOBGIAN.
1911 return (Branch::IsSameExceptForTimes(other) && (*m_rangePtr == *(other->m_rangePtr)));
1912 }
1913
1914 //------------------------------------------------------------------------------------
1915 // Recombinant branches only (returns Branch::NONBRANCH in non-recombinant case).
1916
GetRecPartner() const1917 Branch_ptr RBranch::GetRecPartner() const
1918 {
1919 Branch_ptr partner((Child(0)->Parent(0) == shared_from_this()) ?
1920 Child(0)->Parent(1) : Child(0)->Parent(0));
1921
1922 return partner;
1923
1924 } // GetRecPartner
1925
1926 //------------------------------------------------------------------------------------
1927
ScoreEvent(TreeSummary &,BranchBuffer &) const1928 void RBranch::ScoreEvent(TreeSummary &, BranchBuffer &) const
1929 {
1930 assert(false); // This should never be called.
1931 } // Rbranch::ScoreEvent
1932
1933 //------------------------------------------------------------------------------------
1934 // Third arg is a reference because this function returns results therein.
1935
ScoreEvent(TreeSummary & summary,BranchBuffer & ks,long & s) const1936 void RBranch::ScoreEvent(TreeSummary & summary, BranchBuffer & ks, long & s) const
1937 {
1938 // One interval with a recombination at the top involves *two* RBranches. Only one
1939 // will be summarized into the TreeSummary. Thus, this peculiar-looking code only calls
1940 // AddInterval() if the RBranch in question is its Child's *first* Parent, though
1941 // it does clean-up bookkeeping in any case.
1942 if (Child(0)->Parent(0) == shared_from_this())
1943 {
1944 RecSummary * rsum = dynamic_cast<RecSummary *>(summary.GetSummary(force_REC));
1945 RecTreeSummary & rtreesum = dynamic_cast<RecTreeSummary &>(summary);
1946
1947 // Forces have become inconsistent!?
1948 assert(rsum);
1949
1950 // Obtain the partnerpicks information for the branch which had it, which is not necessarily
1951 // this one. MDEBUG: This code assumes there is at most one local partition force and one
1952 // non-local partition force. It will need to be more complex when multiples of either are allowed.
1953
1954 // I assume that either both parents match the child, and it doesn't matter which we use,
1955 // or one parent does not match the child, in which case that parent must be used.
1956
1957 LongVec1d childpartitions = Child(0)->m_partitions;
1958 LongVec1d otherpartitions = Child(0)->Parent(1)->m_partitions; // other branch of rec
1959 LongVec1d partnerpicks;
1960 LongVec1d lpindex(registry.GetForceSummary().GetLocalPartitionIndexes());
1961 LongVec1d::iterator lpforce;
1962 LongVec1d pickedmembership;
1963
1964 // MDEBUG vestigal code assuming multiple lpforces are possible here.
1965 // This is not carried through correctly elsewhere in function!
1966 for(lpforce = lpindex.begin(); lpforce != lpindex.end(); ++lpforce)
1967 {
1968 if (m_partitions[*lpforce] == childpartitions[*lpforce])
1969 {
1970 pickedmembership = otherpartitions;
1971 }
1972 else
1973 {
1974 pickedmembership = m_partitions;
1975 }
1976 partnerpicks.push_back(pickedmembership[*lpforce]);
1977 rsum->AddToRecombinationCounts(pickedmembership);
1978 }
1979
1980 // Score the event.
1981 rsum->AddInterval(m_eventTime, ks.GetBranchParts(), ks.GetBranchXParts(), s,
1982 FLAGLONG, FLAGLONG, GetRecSite_Br(), partnerpicks, force_REC);
1983
1984 // Update list of recombinations by partition.
1985 rtreesum.AddRecToRecsByPart(pickedmembership, rsum->GetLastAdded());
1986
1987 // Adjust the branch counts and activesite counts for removal of the child.
1988 ks.UpdateBranchCounts(Child(0)->m_partitions, false);
1989 s -= Child(0)->m_rangePtr->NumTargetLinks_Rg();
1990 }
1991
1992 // Adjust the branch and activesite counts for addition of this branch.
1993 ks.UpdateBranchCounts(m_partitions);
1994 s += m_rangePtr->NumTargetLinks_Rg();
1995
1996 } // RBranch::ScoreEvent
1997
1998 //------------------------------------------------------------------------------------
1999 // Debugging function.
2000
CheckInvariant() const2001 bool RBranch::CheckInvariant() const
2002 {
2003 // Recombinant branches have one child...
2004 if (!Child(0))
2005 {
2006 return false;
2007 }
2008
2009 if (Child(1))
2010 {
2011 return false;
2012 }
2013
2014 //...and at least one parent.
2015 if (!Parent(0))
2016 {
2017 return false;
2018 }
2019
2020 if (!Branch::CheckInvariant()) return false;
2021 return true;
2022
2023 } // CheckInvariant
2024
2025 //------------------------------------------------------------------------------------
2026 // Debugging function.
2027
RevalidateRange(FC_Status & fcstatus) const2028 bool RBranch::RevalidateRange(FC_Status & fcstatus) const
2029 {
2030 #if 0 // Include this branch if RevalidateRange is only being called after Prune().
2031 long recsite(m_rangePtr->GetRecSite_Rg());
2032 if (!(Child(0)->m_rangePtr->IsLinkTargettable_Rg(recsite)))
2033 {
2034 cout << "Non-targettable recombination site found in RBranch" << endl;
2035 return false;
2036 }
2037 #endif
2038
2039 rangeset livesites(Intersection(Child(0)->m_rangePtr->GetLiveSites_Rg(), m_rangePtr->GetTransmittedSites_Rg()));
2040 if (!(m_rangePtr->SameLiveSites_Rg(livesites)))
2041 {
2042 return false;
2043 }
2044
2045 // Did we change color inappropriately?
2046 if (m_rangePtr->IsDiseaseSiteTransmitted_Rg())
2047 {
2048 if (Child(0)->m_partitions != m_partitions)
2049 {
2050 cout << "Branch changed color in RBranch" << endl;
2051 return false;
2052 }
2053 }
2054
2055 if (m_rangePtr->LiveSitesOnly_Rg()) return true; // We're done.
2056
2057 #if LAMARC_QA_FC_ON
2058 rangeset target(TargetLinksFrom(Union(m_rangePtr->GetDiseaseSites_Rg(),
2059 RemoveRangeFromRange(fcstatus.Coalesced_Sites(), livesites))));
2060 #else
2061 rangeset target(TargetLinksFrom(Union(m_rangePtr->GetDiseaseSites_Rg(), livesites)));
2062 #endif
2063
2064 if (m_rangePtr->DifferentTargetLinks_Rg(target))
2065 {
2066 return false;
2067 }
2068
2069 rangeset newtarget(RemoveRangeFromRange(m_rangePtr->GetOldTargetLinks_Rg(), target));
2070 if (m_rangePtr->DifferentNewTargetLinks_Rg(newtarget))
2071 {
2072 return false;
2073 }
2074
2075 return true;
2076
2077 } // RBranch::RevalidateRange
2078
2079 //------------------------------------------------------------------------------------
2080
GetGraphMLNodeType() const2081 string RBranch::GetGraphMLNodeType() const
2082 {
2083 return "Rec";
2084 }
2085
2086 //------------------------------------------------------------------------------------
2087
AddNodeInfo(TiXmlElement * nodeElem) const2088 void RBranch::AddNodeInfo(TiXmlElement * nodeElem) const
2089 {
2090 // testing recloc info
2091 TiXmlElement * recElem = new TiXmlElement("data");
2092 nodeElem->LinkEndChild(recElem);
2093 recElem->SetAttribute("key","rec_location");
2094 long recLoc = m_rangePtr->GetRecSite_Rg();
2095 TiXmlText * rlTypeText = new TiXmlText(ToString(recLoc));
2096 recElem->LinkEndChild(rlTypeText);
2097 }
2098
2099 //____________________________________________________________________________________
2100