1 // $Id: tree.cpp,v 1.113 2012/05/22 20:25:32 jyamato 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 <algorithm> // for std::max and min
12 #include <cassert>
13 #include <fstream>
14 #include <limits> // for numeric_limits<double>.epsilon()
15
16 #include "local_build.h"
17
18 #include "branchbuffer.h" // for Tree::Groom
19 #include "datapack.h"
20 #include "dlcell.h"
21 #include "errhandling.h" // for datalikenorm_error in Tree::CalculateDataLikes
22 #include "locus.h"
23 #include "mathx.h" // for exp()
24 #include "registry.h"
25 #include "runreport.h"
26 #include "timemanager.h"
27 #include "tree.h"
28 #include "treesum.h"
29
30 #ifdef DMALLOC_FUNC_CHECK
31 #include "/usr/local/include/dmalloc.h"
32 #endif
33
34 //#define TEST // erynes test
35
36 using namespace std;
37
38 //------------------------------------------------------------------------------------
39
Tree()40 Tree::Tree()
41 : m_overallDL(FLAGDOUBLE),
42 m_randomSource(&(registry.GetRandom())),
43 m_totalSites(0),
44 m_timeManager(NULL)
45 {
46 // Deliberately blank.
47 }
48
49 //------------------------------------------------------------------------------------
50
~Tree()51 Tree::~Tree()
52 {
53 // Deliberately blank due to shared_ptr auto cleanup.
54 delete m_timeManager;
55
56 } // Tree destructor
57
58 //------------------------------------------------------------------------------------
59
60 // MCHECK Can we do this now using Lucian's new mechanism? NB This copy constructor does *not* copy
61 // m_firstInvalid (branch equivalence problem needs to be solved).
62
Tree(const Tree & tree,bool makestump)63 Tree::Tree(const Tree & tree, bool makestump)
64 : m_protoCells(tree.m_protoCells),
65 m_individuals(tree.m_individuals),
66 m_overallDL(tree.m_overallDL),
67 m_randomSource(tree.m_randomSource),
68 m_pLocusVec(tree.m_pLocusVec),
69 m_totalSites(tree.m_totalSites)
70 {
71 const ForceSummary & fs(registry.GetForceSummary());
72 m_timeManager = fs.CreateTimeManager();
73
74 if (!makestump)
75 {
76 *m_timeManager = *(tree.m_timeManager); // Make sure the TimeManager state is the same in both trees.
77 m_timeList = tree.m_timeList; // Deep copy the branchlist.
78
79 // Now setup the tip pointers in the individuals correctly, since the
80 // individual copy ctor only makes empty containers of tips.
81 unsigned long indiv;
82 for(indiv = 0; indiv < m_individuals.size(); ++indiv)
83 {
84 StringVec1d tipnames = tree.m_individuals[indiv].GetAllTipNames();
85 vector<Branch_ptr> itips = GetTips(tipnames);
86 m_individuals[indiv].SetTips(itips);
87 }
88 }
89
90 } // Tree copy ctor
91
92 //------------------------------------------------------------------------------------
93
Clear()94 void Tree::Clear()
95 {
96 m_timeList.Clear();
97 }
98
99 //------------------------------------------------------------------------------------
100
CollectCells()101 const vector<LocusCell> & Tree::CollectCells()
102 {
103 if (m_protoCells.empty())
104 {
105 unsigned long i;
106 for (i = 0; i < m_pLocusVec->size(); ++i)
107 {
108 m_protoCells.push_back((*m_pLocusVec)[i].GetProtoCell());
109 }
110 }
111 return m_protoCells;
112
113 } // CollectCells
114
115 //------------------------------------------------------------------------------------
116
CopyTips(const Tree * tree)117 void Tree::CopyTips(const Tree * tree)
118 {
119 m_timeList.CopyTips(tree->m_timeList);
120 m_firstInvalid = m_timeList.BeginBranch();
121 m_totalSites = tree->m_totalSites;
122
123 m_overallDL = FLAGDOUBLE; // It won't be valid anymore.
124 m_individuals = tree->m_individuals;
125
126 // Now setup the tip pointers in the individuals correctly.
127 unsigned long indiv;
128 for(indiv = 0; indiv < m_individuals.size(); ++indiv)
129 {
130 StringVec1d tipnames = tree->m_individuals[indiv].GetAllTipNames();
131 vector<Branch_ptr> itips = GetTips(tipnames);
132 m_individuals[indiv].SetTips(itips);
133 }
134
135 vector<Locus>::const_iterator locus;
136 m_aliases.clear();
137
138 for(locus = m_pLocusVec->begin(); locus != m_pLocusVec->end(); ++locus)
139 m_aliases.push_back(locus->GetDLCalc()->RecalculateAliases(*this, *locus));
140 }
141
142 //------------------------------------------------------------------------------------
143
CopyBody(const Tree * tree)144 void Tree::CopyBody(const Tree * tree)
145 {
146 m_timeList.CopyBody(tree->m_timeList);
147 m_overallDL = tree->m_overallDL;
148 }
149
150 //------------------------------------------------------------------------------------
151
CopyStick(const Tree * tree)152 void Tree::CopyStick(const Tree * tree)
153 {
154 m_timeManager->CopyStick(*(tree->m_timeManager));
155 }
156
157 //------------------------------------------------------------------------------------
158
CopyPartialBody(const Tree * tree)159 void Tree::CopyPartialBody(const Tree * tree)
160 {
161 // It's wrong to call CopyPartialBody if you don't have a starting point.
162 assert((*m_firstInvalid) != Branch::NONBRANCH);
163
164 // We can't do a partial copy if the tree was changed too high up.
165 if ((*m_firstInvalid)->BranchGroup() == bgroupTip || (*tree->m_firstInvalid)->BranchGroup() == bgroupTip)
166 {
167 m_timeList.CopyBody(tree->m_timeList);
168 }
169 else
170 {
171 m_timeList.CopyPartialBody(tree->m_timeList, tree->m_firstInvalid, m_firstInvalid);
172 }
173
174 m_overallDL = tree->m_overallDL;
175 }
176
177 //------------------------------------------------------------------------------------
178
CreateTip(const TipData & tipdata,const vector<LocusCell> & cells,const vector<LocusCell> & movingcells,const rangeset & diseasesites)179 TBranch_ptr Tree::CreateTip(const TipData & tipdata, const vector<LocusCell> & cells,
180 const vector<LocusCell> & movingcells, const rangeset & diseasesites)
181 {
182 // used in sample only case
183 TBranch_ptr pTip = m_timeList.CreateTip(tipdata, cells, movingcells, m_totalSites, diseasesites);
184 return pTip;
185 }
186
187 //------------------------------------------------------------------------------------
188
CreateTip(const TipData & tipdata,const vector<LocusCell> & cells,const vector<LocusCell> & movingcells,const rangeset & diseasesites,const vector<Locus> & loci)189 TBranch_ptr Tree::CreateTip(const TipData & tipdata, const vector<LocusCell> & cells,
190 const vector<LocusCell> & movingcells, const rangeset & diseasesites,
191 const vector<Locus> & loci)
192 {
193 // used in panel case
194 TBranch_ptr pTip = m_timeList.CreateTip(tipdata, cells, movingcells, m_totalSites, diseasesites,
195 loci);
196 return pTip;
197 }
198 //------------------------------------------------------------------------------------
199
SetTreeTimeManager(TimeManager * tm)200 void Tree::SetTreeTimeManager(TimeManager * tm)
201 {
202 m_timeManager = tm;
203 }
204
205 //------------------------------------------------------------------------------------
206
SetIndividualsWithTips(const vector<Individual> & indvec)207 void Tree::SetIndividualsWithTips(const vector<Individual> & indvec)
208 {
209 m_individuals = indvec;
210 } // SetIndividualsWithTips
211
212 //------------------------------------------------------------------------------------
213
SetLocusVec(vector<Locus> * loc)214 void Tree::SetLocusVec(vector<Locus> * loc)
215 {
216 m_pLocusVec = loc;
217 m_totalSites = 0;
218 vector<Locus>::const_iterator locit;
219 for(locit = loc->begin(); locit != loc->end(); ++locit)
220 {
221 if (locit->GetSiteSpan().second > m_totalSites)
222 m_totalSites = locit->GetSiteSpan().second;
223 }
224
225 } // SetLocusVec
226
227 //------------------------------------------------------------------------------------
228
SummarizeTree() const229 TreeSummary * Tree::SummarizeTree() const
230 {
231 TreeSummary * treesum = registry.GetProtoTreeSummary().Clone();
232 treesum->Summarize(*this);
233 return treesum;
234
235 } // SummarizeTree
236
237 //------------------------------------------------------------------------------------
238
CalculateDataLikesForFixedLoci()239 double Tree::CalculateDataLikesForFixedLoci()
240 {
241 unsigned long loc;
242 #ifndef NDEBUG
243 vector<double> likes;
244 #endif
245 double likelihood, overalllike(0);
246
247 for (loc = 0; loc < m_pLocusVec->size(); ++loc)
248 {
249 const Locus & locus = (*m_pLocusVec)[loc];
250 try
251 { // Check for need to switch to normalization.
252 likelihood = locus.CalculateDataLikelihood(*this, false);
253 }
254 catch (datalikenorm_error & ex) // Normalization is set by thrower.
255 {
256 m_timeList.SetAllUpdateDLs(); // All subsequent loci will recalculate the entire tree.
257 RunReport & runreport = registry.GetRunReport();
258 runreport.ReportChat("\n", 0);
259 runreport.ReportChat("Subtree of likelihood 0.0 found: Turning on normalization and re-calculating.");
260
261 likelihood = locus.CalculateDataLikelihood(*this, false);
262 }
263 #ifndef NDEBUG
264 likes.push_back(likelihood);
265 #endif
266 //cout << "short loc: " << loc << " likelihood: " << likelihood << endl;
267 overalllike += likelihood; // Add to accumulator.
268
269 #if LIKETRACK
270 ofstream of;
271 of.open("likes1", ios::app);
272 of << GetDLValue() << " ";
273 of.close();
274 #endif
275 }
276
277 #ifndef NDEBUG
278 // This code tests a full likelihood recalculation against the time-saving partial recalculation. Slow but sure.
279 double checkDL = 0.0;
280 m_timeList.SetAllUpdateDLs();
281 for (loc = 0; loc < m_pLocusVec->size(); ++loc)
282 {
283 const Locus & locus = (*m_pLocusVec)[loc];
284 double likelihood = locus.CalculateDataLikelihood(*this, false);
285 //cout << "long loc: " << loc << " likelihood: " << likelihood << endl << endl ;
286 checkDL += likelihood;
287 #ifndef LAMARC_QA_SINGLE_DENOVOS
288 // likelihood can be zero in single denovos test
289 assert(fabs(likelihood - likes[loc]) / likelihood < EPSILON);
290 #endif // LAMARC_QA_SINGLE_DENOVOS
291 }
292 #ifndef LAMARC_QA_SINGLE_DENOVOS
293 // overalllike can be zero in single denovos test
294 assert(fabs(overalllike - checkDL) / overalllike < EPSILON);
295 #endif // LAMARC_QA_SINGLE_DENOVOS
296 #endif // NDEBUG
297
298 return overalllike;
299
300 } // CalculateDataLikesForFixedLoci
301
302 //------------------------------------------------------------------------------------
303
SetupAliases(const vector<Locus> & loci)304 void Tree::SetupAliases(const vector<Locus> & loci)
305 {
306 // Set up the aliases; this must be done after filling the tips with data!
307 m_aliases.clear();
308 vector<Locus>::const_iterator locus;
309
310 for (locus = loci.begin(); locus != loci.end(); ++locus)
311 {
312 m_aliases.push_back(locus->GetDLCalc()->RecalculateAliases(*this, *locus));
313 }
314
315 } // SetupAliases
316
317 //------------------------------------------------------------------------------------
318
ActivateTips(Tree * othertree)319 vector<Branch_ptr> Tree::ActivateTips(Tree * othertree)
320 {
321 vector<Branch_ptr> tips;
322 m_firstInvalid = m_timeList.BeginBranch();
323 othertree->SetFirstInvalid();
324
325 m_timeList.ClearBody();
326
327 Branchiter brit;
328 for (brit = m_timeList.FirstTip(); brit != m_timeList.EndBranch(); brit = m_timeList.NextTip(brit))
329 {
330 tips.push_back(*brit);
331 }
332
333 return tips;
334 }
335
336 //------------------------------------------------------------------------------------
337
ActivateBranch(Tree * othertree)338 Branch_ptr Tree::ActivateBranch(Tree * othertree)
339 {
340 // NCuttable - 1 because the root branch is of a cuttable type but is not, in fact, a cuttable branch.
341 // + 1 because it's helpful to count cuttable branches from 1, not 0.
342 long rn = m_randomSource->Long(m_timeList.GetNCuttable() - 1) + 1;
343
344 Branchiter brit;
345 Branch_ptr pActive = Branch::NONBRANCH;
346 assert(m_timeList.IsValidTimeList());
347
348 for (brit = m_timeList.BeginBranch(); rn > 0; ++brit)
349 {
350 m_firstInvalid = brit;
351 pActive = *brit;
352 rn -= pActive->Cuttable();
353 }
354
355 Branchconstiter constfirstinvalid = m_firstInvalid;
356 othertree->SetFirstInvalidFrom(constfirstinvalid);
357
358 assert(pActive != Branch::NONBRANCH); // We didn't find any cuttable branches?
359
360 Break(pActive);
361 pActive->SetParent(0, Branch::NONBRANCH);
362 pActive->SetParent(1, Branch::NONBRANCH);
363 return pActive;
364 }
365
366 //------------------------------------------------------------------------------------
367
SetFirstInvalidFrom(Branchconstiter & target)368 void Tree::SetFirstInvalidFrom(Branchconstiter & target)
369 {
370 m_firstInvalid = m_timeList.FindEquivBranch(*target);
371 } // SetFirstInvalidFrom
372
373 //------------------------------------------------------------------------------------
374
ActivateRoot(FC_Status & fcstatus)375 Branch_ptr Tree::ActivateRoot(FC_Status & fcstatus)
376 {
377 Branch_ptr pBase = m_timeList.Base(); // Get the tree base and root.
378 Branch_ptr pRoot = pBase->Child(0);
379 pBase->SetChild(0, Branch::NONBRANCH); // Disconnect them.
380 pRoot->SetParent(0, Branch::NONBRANCH);
381
382 // Update the fcstatus with root info.
383 #if LAMARC_QA_FC_ON
384 if (pRoot->Child(1)) // We have two children.
385 {
386 rangeset sites(Intersection(pRoot->Child(0)->GetLiveSites_Br(), pRoot->Child(1)->GetLiveSites_Br()));
387 fcstatus.Decrement_FC_Counts(sites);
388 }
389 #endif
390
391 #if LAMARC_QA_FC_ON
392 pRoot->UpdateRootBranchRange(fcstatus.Coalesced_Sites(), true); // Update the root's range object.
393 #else
394 rangeset emptyset;
395 pRoot->UpdateRootBranchRange(emptyset, false);
396 #endif
397 return pRoot;
398 }
399
400 //------------------------------------------------------------------------------------
401
ChoosePreferentiallyTowardsRoot(Tree * othertree)402 Branch_ptr Tree::ChoosePreferentiallyTowardsRoot(Tree * othertree)
403 {
404 // Choose a random body branch, triangle weighted towards the root.
405 long nbranches = m_timeList.GetNBodyBranches();
406 long weightedbranches = (nbranches * (nbranches + 1)) / 2;
407 // weightedbranches = 0; //JRM remove
408 long chosenweight = m_randomSource->Long(weightedbranches) + 1;
409 long chosenbranch(0L), weight(0L);
410 while(weight < chosenweight)
411 {
412 ++chosenbranch;
413 weight += chosenbranch;
414 }
415
416 long thisbranch(1L);
417 for(m_firstInvalid = m_timeList.FirstBody();
418 m_firstInvalid != m_timeList.EndBranch();
419 m_firstInvalid = m_timeList.NextBody(m_firstInvalid), ++thisbranch)
420 {
421 if (thisbranch == chosenbranch) break;
422 }
423 assert(m_firstInvalid != m_timeList.EndBranch() && thisbranch == chosenbranch);
424
425 // m_firstInvalid may be set incorrectly when there are multiple branches with the same eventtime.
426 m_firstInvalid = m_timeList.ResolveTiedTimes(m_firstInvalid);
427
428 Branchconstiter constfirstinvalid = m_firstInvalid;
429 othertree->SetFirstInvalidFrom(constfirstinvalid);
430
431 return *m_firstInvalid;
432 }
433
434 //------------------------------------------------------------------------------------
435
ChooseFirstBranchInEpoch(double targettime,Tree * othertree)436 Branch_ptr Tree::ChooseFirstBranchInEpoch(double targettime, Tree * othertree)
437 {
438 for(m_firstInvalid = m_timeList.FirstBody();
439 m_firstInvalid != m_timeList.EndBranch();
440 m_firstInvalid = m_timeList.NextBody(m_firstInvalid))
441 {
442 if ((*m_firstInvalid)->GetTime() > targettime) break;
443 }
444
445 assert(m_firstInvalid != m_timeList.EndBranch());
446
447 // m_firstInvalid may be set incorrectly when there are multiple branches with the same eventtime.
448 m_firstInvalid = m_timeList.ResolveTiedTimes(m_firstInvalid);
449
450 Branchconstiter constfirstinvalid = m_firstInvalid;
451 othertree->SetFirstInvalidFrom(constfirstinvalid);
452 return *m_firstInvalid;
453 }
454
455 //------------------------------------------------------------------------------------
456
AttachBase(Branch_ptr newroot)457 void Tree::AttachBase(Branch_ptr newroot)
458 {
459 Branch_ptr pBase = m_timeList.Base(); // Get the tree base.
460 pBase->SetChild(0, newroot); // Link the root to the base.
461 newroot->SetParent(0, pBase);
462 }
463
464 //------------------------------------------------------------------------------------
465
FindAllBranchesAtTime(double eventT)466 vector<Branch_ptr> Tree::FindAllBranchesAtTime(double eventT)
467 {
468 Branch_ptr pBranch;
469 vector<Branch_ptr> inactives;
470
471 Branchiter brit;
472 for (brit = m_timeList.BeginBranch(); brit != m_timeList.EndBranch(); ++brit)
473 {
474 pBranch = *brit;
475 if (eventT < pBranch->m_eventTime)
476 {
477 break;
478 }
479
480 if (pBranch->Parent(0))
481 {
482 if (eventT < pBranch->Parent(0)->m_eventTime)
483 {
484 inactives.push_back(pBranch);
485 }
486 }
487 else assert(!pBranch->Parent(1));
488 }
489
490 return inactives;
491 }
492
493 //------------------------------------------------------------------------------------
494
FirstInterval(double eventT)495 vector<Branch_ptr> Tree::FirstInterval(double eventT)
496 {
497 return FindAllBranchesAtTime(eventT);
498 }
499
500 //------------------------------------------------------------------------------------
501
FindBranchesImmediatelyTipwardOf(Branchiter startbr)502 vector<Branch_ptr> Tree::FindBranchesImmediatelyTipwardOf(Branchiter startbr)
503 {
504 assert((*startbr)->m_eventTime != 0.0);
505 Branchiter prevbr(startbr);
506 --prevbr; // Assumes we have reversible iterator.
507 while((*prevbr)->m_eventTime == (*startbr)->m_eventTime)
508 {
509 assert(prevbr != m_timeList.BeginBranch()); // Too far!
510 --prevbr;
511 }
512 return FindAllBranchesAtTime((*prevbr)->m_eventTime);
513 }
514
515 //------------------------------------------------------------------------------------
516
FindBranchesBetween(double startint,double endint)517 vector<Branch_ptr> Tree::FindBranchesBetween(double startint, double endint)
518 {
519 vector<Branch_ptr> branches;
520
521 Branchiter brit;
522 for (brit = m_timeList.BeginBranch(); brit != m_timeList.EndBranch(); ++brit)
523 {
524 Branch_ptr pBranch = *brit;
525 if (startint <= pBranch->m_eventTime && pBranch->m_eventTime < endint)
526 {
527 branches.push_back(pBranch);
528 }
529 }
530
531 return branches;
532 }
533
534 //------------------------------------------------------------------------------------
535
NextInterval(Branch_ptr branch)536 void Tree::NextInterval(Branch_ptr branch)
537 {
538 // Deliberately blank; needed only in rectree.
539 }
540
541 //------------------------------------------------------------------------------------
542
Prune()543 void Tree::Prune()
544 {
545 m_timeList.Prune();
546 // This is done last so that the information it overwrites can be used during
547 // pruning and for validation--do not move it into the main loop!
548 Branchiter brit;
549 for (brit = m_timeList.FirstBody(); brit != m_timeList.EndBranch(); brit = m_timeList.NextBody(brit))
550 {
551 (*brit)->ResetBuffersForNextRearrangement();
552 }
553 TrimStickToRoot();
554 }
555
556 //------------------------------------------------------------------------------------
557
TrimStickToRoot()558 void Tree::TrimStickToRoot()
559 {
560 m_timeManager->ChopOffStickAt(m_timeList.RootTime());
561 } // TrimStickToRoot
562
563 //------------------------------------------------------------------------------------
564
Coalesce(Branch_ptr child1,Branch_ptr child2,double tevent,const rangeset & fcsites)565 Branch_ptr Tree::Coalesce(Branch_ptr child1, Branch_ptr child2, double tevent, const rangeset & fcsites)
566 {
567 assert(child1->m_partitions == child2->m_partitions);
568
569 bool newbranchisinactive(false);
570 Branch_ptr parent = Branch_ptr(new CBranch(child1->m_rangePtr, child2->m_rangePtr, newbranchisinactive, fcsites));
571 parent->m_eventTime = tevent;
572 parent->CopyPartitionsFrom(child1);
573
574 parent->SetChild(0, child1);
575 parent->SetChild(1, child2);
576 child1->SetParent(0, parent);
577 child2->SetParent(0, parent);
578
579 parent->SetDLCells(CollectCells());
580 parent->SetUpdateDL(); // Mark this branch for updating.
581
582 m_timeList.Collate(parent);
583 assert(child1->m_partitions == parent->m_partitions);
584
585 return parent;
586
587 } // Tree::Coalesce
588
589 //------------------------------------------------------------------------------------
590
CoalesceActive(double eventT,Branch_ptr pActive1,Branch_ptr pActive2,const rangeset & fcsites)591 Branch_ptr Tree::CoalesceActive(double eventT, Branch_ptr pActive1, Branch_ptr pActive2, const rangeset & fcsites)
592 {
593 // Create the parent of the two branches.
594 Branch_ptr pParent = Coalesce(pActive1, pActive2, eventT, fcsites);
595 return pParent;
596 }
597
598 //------------------------------------------------------------------------------------
599
CoalesceInactive(double eventT,Branch_ptr pActive,Branch_ptr pInactive,const rangeset & fcsites)600 Branch_ptr Tree::CoalesceInactive(double eventT, Branch_ptr pActive, Branch_ptr pInactive, const rangeset & fcsites)
601 {
602 // Create the parent of the two branches.
603 bool newbranchisinactive(true);
604 Branch_ptr pParent = Branch_ptr(new CBranch(pInactive->m_rangePtr, pActive->m_rangePtr, newbranchisinactive, fcsites));
605 pParent->m_eventTime = eventT;
606 pParent->CopyPartitionsFrom(pActive);
607
608 // Hook children to parent.
609 pParent->SetChild(0, pInactive);
610 pParent->SetChild(1, pActive);
611
612 // Hook pInactive's old parent up to new parent.
613 pInactive->Parent(0)->ReplaceChild(pInactive, pParent);
614 pInactive->SetParent(0, pParent);
615 if (pInactive->Parent(1))
616 {
617 pInactive->Parent(1)->ReplaceChild(pInactive, pParent);
618 pInactive->SetParent(1, Branch::NONBRANCH);
619 }
620
621 pActive->SetParent(0, pParent);
622
623 pParent->SetDLCells(CollectCells());
624 m_timeList.SetUpdateDLs(pParent); // Set the update data likelihood flag.
625
626 m_timeList.Collate(pParent);
627 return pParent;
628 }
629
630 //------------------------------------------------------------------------------------
631
Migrate(double eventT,long pop,long maxEvents,Branch_ptr pActive)632 Branch_ptr Tree::Migrate(double eventT, long pop, long maxEvents, Branch_ptr pActive)
633 {
634 if (m_timeList.HowMany(btypeMig) >= maxEvents)
635 {
636 mig_overrun e;
637 throw e;
638 }
639
640 force_type force;
641 if (registry.GetForceSummary().CheckForce(force_MIG))
642 {
643 force = force_MIG;
644 }
645 else
646 {
647 force = force_DIVMIG;
648 }
649
650 Branch_ptr pParent;
651 if (force == force_MIG)
652 {
653 pParent = Branch_ptr(new MBranch(pActive->m_rangePtr));
654 }
655 else
656 {
657 pParent = Branch_ptr(new DivMigBranch(pActive->m_rangePtr));
658 }
659
660 pParent->m_eventTime = eventT;
661 pParent->CopyPartitionsFrom(pActive);
662 pParent->SetPartition(force, pop);
663
664 pParent->SetChild(0, pActive);
665 pActive->SetParent(0, pParent);
666
667 m_timeList.Collate(pParent);
668 return pParent;
669 }
670
671 //------------------------------------------------------------------------------------
672
DiseaseMutate(double eventT,long endstatus,long maxevents,Branch_ptr pActive)673 Branch_ptr Tree::DiseaseMutate(double eventT, long endstatus, long maxevents, Branch_ptr pActive)
674 {
675 if (m_timeList.HowMany(btypeDisease) >= maxevents)
676 {
677 dis_overrun e;
678 throw e;
679 }
680
681 Branch_ptr pParent = Branch_ptr(new DBranch(pActive->m_rangePtr));
682 pParent->m_eventTime = eventT;
683 pParent->CopyPartitionsFrom(pActive);
684 pParent->SetPartition(force_DISEASE, endstatus);
685
686 pParent->SetChild(0, pActive);
687 pActive->SetParent(0, pParent);
688
689 m_timeList.Collate(pParent);
690
691 return pParent;
692
693 } // DiseaseMutate
694
695 //------------------------------------------------------------------------------------
696
TransitionEpoch(double eventT,long newpop,long maxevents,Branch_ptr pActive)697 Branch_ptr Tree::TransitionEpoch(double eventT, long newpop, long maxevents, Branch_ptr pActive)
698 {
699 if (m_timeList.HowMany(btypeEpoch) >= maxevents)
700 {
701 epoch_overrun e;
702 throw e;
703 }
704
705 Branch_ptr pParent = Branch_ptr(new EBranch(pActive->m_rangePtr));
706 pParent->m_eventTime = eventT;
707 pParent->CopyPartitionsFrom(pActive);
708 pParent->SetPartition(force_DIVMIG, newpop);
709
710 pParent->SetChild(0, pActive);
711 pActive->SetParent(0, pParent);
712
713 m_timeList.Collate(pParent);
714 return pParent;
715
716 } // TransitionEpoch
717
718 //------------------------------------------------------------------------------------
719
FindBranchesStartingOnOpenInterval(double starttime,double endtime)720 vector<Branch_ptr> Tree::FindBranchesStartingOnOpenInterval(
721 double starttime, double endtime)
722 {
723 vector<Branch_ptr> branches;
724 Branchiter br;
725 for(br = m_timeList.BeginBranch(); br != m_timeList.EndBranch(); ++br)
726 {
727 if ((*br)->m_eventTime <= starttime) continue;
728 if ((*br)->m_eventTime >= endtime) break;
729 branches.push_back(*br);
730 }
731
732 return branches;
733
734 } // FindBranchesStartingOnOpenInterval
735
736 //------------------------------------------------------------------------------------
737
FindEpochBranchesAt(double time)738 vector<Branch_ptr> Tree::FindEpochBranchesAt(double time)
739 {
740 vector<Branch_ptr> branches;
741 Branchiter br;
742 for(br = m_timeList.BeginBranch(); br != m_timeList.EndBranch(); ++br)
743 {
744 if ((*br)->m_eventTime == time && (*br)->Event() == btypeEpoch)
745 branches.push_back(*br);
746 }
747
748 return branches;
749
750 } // FindEpochBranchesAt
751
752 //------------------------------------------------------------------------------------
753
FindBranchesStartingRootwardOf(double time)754 vector<Branch_ptr> Tree::FindBranchesStartingRootwardOf(double time)
755 {
756 vector<Branch_ptr> branches;
757 Branchiter br;
758 for(br = m_timeList.BeginBranch(); br != m_timeList.EndBranch(); ++br)
759 {
760 if ((*br)->m_eventTime > time)
761 branches.push_back(*br);
762 }
763
764 return branches;
765
766 } // FindBranchesStartingRootwardOf
767
768 //------------------------------------------------------------------------------------
769
SwapSiteDLs()770 void Tree::SwapSiteDLs()
771 {
772 Individual ind;
773
774 // Pick an individual with phase-unknown sites at random.
775 do {
776 ind = m_individuals[m_randomSource->Long(m_individuals.size())];
777 } while (!ind.AnyPhaseUnknownSites());
778
779 // Pick a phase marker at random.
780 pair<long, long> locusmarker = ind.PickRandomPhaseMarker(*m_randomSource);
781 long locus = locusmarker.first;
782 long posn = locusmarker.second;
783
784 // Pick two tips in the individual at random.
785 vector<Branch_ptr> haps = ind.GetAllTips();
786 long size = haps.size();
787 long rnd1 = m_randomSource->Long(size);
788 long rnd2 = m_randomSource->Long(size - 1);
789 if (rnd1 <= rnd2)
790 rnd2++;
791
792 Branch_ptr pTip1 = haps[rnd1];
793 Branch_ptr pTip2 = haps[rnd2];
794
795 // NB: Swap first DLCell only!
796
797 assert(pTip1->GetNcells(0) == pTip2->GetNcells(0));
798
799 // Get the data likelihood arrays.
800 Cell_ptr dlcell1 = pTip1->GetDLCell(locus, markerCell, false);
801 Cell_ptr dlcell2 = pTip2->GetDLCell(locus, markerCell, false);
802
803 // Swap the cells at the marker position.
804 dlcell1->SwapDLs(dlcell2, posn);
805
806 // Take care of needed data-likelihood corrections.
807 const Locus & loc = (*m_pLocusVec)[locus];
808 m_aliases[locus] = loc.GetDLCalc()->RecalculateAliases(*this, loc);
809
810 MarkForDLRecalc(pTip1);
811 MarkForDLRecalc(pTip2);
812
813 } // Tree::SwapSiteDLs
814
815 //------------------------------------------------------------------------------------
816 // PickNewSiteDLs is used when haplotype arranging for non-50/50 haplotype probabilities.
817 // Otherwise, SwapSiteDLs is used.
818
PickNewSiteDLs()819 void Tree::PickNewSiteDLs()
820 {
821 long ind;
822
823 // Pick an individual with phase-unknown sites at random.
824 do {
825 ind = m_randomSource->Long(m_individuals.size());
826 } while (!m_individuals[ind].MultipleTraitHaplotypes());
827
828 // Pick a locus/marker at random.
829 pair<string, long> locusmarker = m_individuals[ind].PickRandomHaplotypeMarker();
830 string lname = locusmarker.first;
831 long marker = locusmarker.second;
832
833 // Have the individual in question pick up a new set of haplotypes.
834 m_individuals[ind].ChooseNewHaplotypesFor(lname, marker);
835
836 // Finish this routine in a subroutine, because chances are good that the locus in question
837 // is a member of m_pMovingLocusVec, which we only have if we're a RecTree.
838
839 ReassignDLsFor(lname, marker, ind);
840 return;
841
842 } // PickNewSiteDLs
843
844 //------------------------------------------------------------------------------------
845
ReassignDLsFor(string lname,long marker,long ind)846 void Tree::ReassignDLsFor(string lname, long marker, long ind)
847 {
848 // Find out which locus we're dealing with. We have its name. If the locus in question
849 // is in m_pMovingLocusVec, we've already dealt with it in RecTree::ReassignDLsFor.
850 long locus = FLAGLONG;
851 for (unsigned long lnum = 0; lnum < m_pLocusVec->size(); lnum++)
852 {
853 if ((*m_pLocusVec)[lnum].GetName() == lname)
854 {
855 locus = lnum;
856 }
857 }
858 if (locus == FLAGLONG)
859 {
860 throw implementation_error("We have the segment name " + lname
861 + ", but no segment matches it. Did something happen"
862 " to this tree? To the segment in question?");
863 }
864
865 vector<Branch_ptr> haps = m_individuals[ind].GetAllTips();
866 vector<LocusCell> cells = m_individuals[ind].GetLocusCellsFor(lname, marker);
867
868 for (unsigned long tip = 0; tip < haps.size(); tip++)
869 {
870 Cell_ptr origcell = haps[tip]->GetDLCell(locus, markerCell, false);
871 Cell_ptr newcell = cells[tip][0];
872 origcell->SetSiteDLs(marker, newcell->GetSiteDLs(marker));
873 MarkForDLRecalc(haps[tip]);
874 }
875
876 // take care of needed data-likelihood corrections
877 const Locus & loc = (*m_pLocusVec)[locus];
878 m_aliases[locus] = loc.GetDLCalc()->RecalculateAliases(*this, loc);
879
880 } // Tree::ReassignDLsFor
881
882 //------------------------------------------------------------------------------------
883
MarkForDLRecalc(Branch_ptr markbr)884 void Tree::MarkForDLRecalc(Branch_ptr markbr)
885 {
886 markbr->SetUpdateDL();
887 markbr->MarkParentsForDLCalc();
888
889 } // Tree::MarkForDLRecalc
890
891 //------------------------------------------------------------------------------------
892 // This version is used to limit a multi-locus case to just the locus under consideration.
893
GetLocusSubtrees(rangepair span) const894 rangevector Tree::GetLocusSubtrees(rangepair span) const
895 {
896 rangevector subtrees;
897 subtrees.push_back(span);
898 return subtrees;
899
900 } // GetLocusSubtrees
901
902 //------------------------------------------------------------------------------------
903
GetTips(StringVec1d & tipnames) const904 vector<Branch_ptr> Tree::GetTips(StringVec1d & tipnames) const
905 {
906 vector<Branch_ptr> tips;
907 StringVec1d::iterator tname;
908 for(tname = tipnames.begin(); tname != tipnames.end(); ++tname)
909 tips.push_back(GetTip(*tname));
910
911 return tips;
912
913 } // GetTips
914
915 //------------------------------------------------------------------------------------
916
GetTip(const string & name) const917 Branch_ptr Tree::GetTip(const string & name) const
918 {
919 return m_timeList.GetTip(name);
920 } // GetTip
921
922 //------------------------------------------------------------------------------------
923
NoPhaseUnknownSites() const924 bool Tree::NoPhaseUnknownSites() const
925 {
926 IndVec::const_iterator ind;
927 unsigned long locus;
928 for (locus = 0; locus < m_pLocusVec->size(); ++locus)
929 {
930 for(ind = m_individuals.begin(); ind != m_individuals.end(); ++ind)
931 if (ind->AnyPhaseUnknownSites()) return false;
932 }
933
934 return true;
935
936 } // Tree:NoPhaseUnknownSites
937
938 //------------------------------------------------------------------------------------
939
GetNsites() const940 long Tree::GetNsites() const
941 {
942 long nsites(0L);
943 vector<Locus>::iterator locus;
944 for(locus = m_pLocusVec->begin() ; locus != m_pLocusVec->end(); ++locus)
945 {
946 nsites += locus->GetNsites();
947 }
948
949 return nsites;
950
951 } // Tree::GetNsites
952
953 //------------------------------------------------------------------------------------
954 // Debugging function.
955
IsValidTree() const956 bool Tree::IsValidTree() const
957 {
958 // NB: This can only be called on a completed tree, not one in mid-rearrangement.
959 // Each Branch validates itself.
960 Branchconstiter brit = m_timeList.BeginBranch();
961 long ncuttable = 0;
962
963 for ( ; brit != m_timeList.EndBranch(); ++brit)
964 {
965 //LS TEST
966 //(*brit)->PrintInfo_Br();
967 if (!(*brit)->CheckInvariant())
968 {
969 return false;
970 }
971 ncuttable += (*brit)->Cuttable();
972 }
973
974 if (ncuttable != m_timeList.GetNCuttable())
975 {
976 return false;
977 }
978 return m_timeList.IsValidTimeList();
979
980 // if (FindContradictParts()) return false;
981 // return true;
982
983 } // IsValidTree
984
985 //------------------------------------------------------------------------------------
986
ConsistentWithParameters(const ForceParameters & fp) const987 bool Tree::ConsistentWithParameters(const ForceParameters& fp) const
988 {
989 Branchconstiter brit;
990 vector<double> epochtimes = fp.GetEpochTimes();
991 for (brit=m_timeList.FirstBody(); brit != m_timeList.EndBranch();
992 brit = m_timeList.NextBody(brit))
993 {
994 if ((*brit)->Event() == btypeEpoch) {
995 if (find(epochtimes.begin(), epochtimes.end(),
996 (*brit)->m_eventTime) == epochtimes.end()) {
997 return false; // tree inconsistent with parameters
998 }
999 }
1000 }
1001 return true;
1002 }
1003
1004 //------------------------------------------------------------------------------------
1005
operator ==(const Tree & src) const1006 bool Tree::operator==(const Tree & src) const
1007 {
1008 if (m_pLocusVec->size() != src.m_pLocusVec->size()) return false;
1009 unsigned long i;
1010 for (i = 0; i < m_pLocusVec->size(); ++i)
1011 {
1012 if ((*m_pLocusVec)[i].GetNsites() != (*src.m_pLocusVec)[i].GetNsites()) return false;
1013 }
1014 if (m_overallDL != src.m_overallDL) return false;
1015 if (m_timeList != src.m_timeList) return false;
1016
1017 // WARNING I don't guarantee these are the same tree, especially if
1018 // we are in mid-rearrangement, but they are at least very similar. --Mary
1019 return true;
1020
1021 } // operator==
1022
1023 //------------------------------------------------------------------------------------
1024
SimulateDataIfNeeded()1025 bool Tree::SimulateDataIfNeeded()
1026 {
1027 bool simulated = false;
1028 for (unsigned long loc = 0; loc < m_pLocusVec->size(); loc++)
1029 {
1030 Locus & locus = (*m_pLocusVec)[loc];
1031 if (locus.GetShouldSimulate())
1032 {
1033 simulated = true;
1034 locus.SimulateData(*this, m_totalSites);
1035 }
1036 }
1037
1038 // Redo the aliases.
1039 SetupAliases(*m_pLocusVec);
1040
1041 // Note that we have to calculate all DLCells.
1042 m_timeList.SetAllUpdateDLs();
1043
1044 return simulated;
1045 }
1046
1047 //------------------------------------------------------------------------------------
1048 // Debugging function.
1049
DLCheck(const Tree & other) const1050 void Tree::DLCheck(const Tree & other) const
1051 {
1052 cout << m_timeList.DLCheck(other.GetTimeList()) << endl;
1053 } // DLCheck
1054
1055 //------------------------------------------------------------------------------------
1056 // Debugging function.
1057
PrintStickThetasToFile(ofstream & of) const1058 void Tree::PrintStickThetasToFile(ofstream & of) const
1059 {
1060 m_timeManager->PrintStickThetasToFile(of);
1061 } // PrintStickThetasToFile
1062
1063 //------------------------------------------------------------------------------------
1064 // Debugging function.
1065
PrintStickFreqsToFile(ofstream & of) const1066 void Tree::PrintStickFreqsToFile(ofstream & of) const
1067 {
1068 m_timeManager->PrintStickFreqsToFile(of);
1069 } // PrintStickFreqsToFile
1070
1071 //------------------------------------------------------------------------------------
1072 // Debugging function.
1073
PrintStickFreqsToFileAtTime(ofstream & of,double time) const1074 void Tree::PrintStickFreqsToFileAtTime(ofstream & of, double time) const
1075 {
1076 m_timeManager->PrintStickFreqsToFileAtTime(of, time);
1077 } // PrintStickFreqsToFileAtTime
1078
1079 //------------------------------------------------------------------------------------
1080 // Debugging function.
1081
PrintStickThetasToFileForJoint300(ofstream & of) const1082 void Tree::PrintStickThetasToFileForJoint300(ofstream & of) const
1083 {
1084 m_timeManager->PrintStickThetasToFileForJoint300(of);
1085 } // PrintStickThetasToFileForJoint300
1086
1087 //------------------------------------------------------------------------------------
1088 // Debugging function.
1089
PrintStickToFile(ofstream & of) const1090 void Tree::PrintStickToFile(ofstream & of) const
1091 {
1092 m_timeManager->PrintStickToFile(of);
1093 } // PrintStickAndSummaryToFile
1094
1095 //------------------------------------------------------------------------------------
1096
PrintDirectionalMutationEventCountsToFile(ofstream & of) const1097 void Tree::PrintDirectionalMutationEventCountsToFile(ofstream & of) const
1098 {
1099 m_timeList.PrintDirectionalMutationEventCountsToFile(of);
1100 } // PrintDirectionalMutationEventCountsToFile
1101
1102 //------------------------------------------------------------------------------------
1103 // Debugging function.
1104
PrintTimeTilFirstEventToFile(ofstream & of) const1105 void Tree::PrintTimeTilFirstEventToFile(ofstream & of) const
1106 {
1107 m_timeList.PrintTimeTilFirstEventToFile(of);
1108 } // PrintTimeTilFirstEventToFile
1109
1110 //------------------------------------------------------------------------------------
1111 // Debugging function.
1112
PrintTraitPhenotypeAtLastCoalescence(ofstream & of) const1113 void Tree::PrintTraitPhenotypeAtLastCoalescence(ofstream & of) const
1114 {
1115 m_timeList.PrintTraitPhenotypeAtLastCoalescence(of);
1116 } // PrintTraitPhenotypeAtLastCoalescence
1117
1118 //------------------------------------------------------------------------------------
1119
DestroyStick()1120 void Tree::DestroyStick()
1121 {
1122 m_timeManager->ClearStick();
1123 } // DestroyStick
1124
1125 //------------------------------------------------------------------------------------
1126
SetStickParams(const ForceParameters & fp)1127 void Tree::SetStickParams(const ForceParameters & fp)
1128 {
1129 m_timeManager->SetStickParameters(fp);
1130 } // SetStickParams
1131
1132 //------------------------------------------------------------------------------------
1133
UsingStick() const1134 bool Tree::UsingStick() const
1135 {
1136 return m_timeManager->UsingStick();
1137 } // UsingStick
1138
1139 //------------------------------------------------------------------------------------
1140
ScoreStick(TreeSummary & treesum) const1141 void Tree::ScoreStick(TreeSummary & treesum) const
1142 {
1143 m_timeManager->ScoreStick(treesum);
1144 } // ScoreStick
1145
1146 //------------------------------------------------------------------------------------
1147
XpartThetasAtT(double time,const ForceParameters & fp) const1148 DoubleVec1d Tree::XpartThetasAtT(double time, const ForceParameters & fp) const
1149 {
1150 return m_timeManager->XpartThetasAtT(time, fp);
1151 } // XpartThetasAtT
1152
1153 //------------------------------------------------------------------------------------
1154
PartitionThetasAtT(double time,const force_type force,const ForceParameters & fp) const1155 DoubleVec1d Tree::PartitionThetasAtT(double time, const force_type force, const ForceParameters & fp) const
1156 {
1157 return m_timeManager->PartitionThetasAtT(time, force, fp);
1158 } // XpartThetasAtT
1159
1160 //------------------------------------------------------------------------------------
1161
SetNewTimesFrom(Branchiter startpoint,const DoubleVec1d & newtimes)1162 void Tree::SetNewTimesFrom(Branchiter startpoint, const DoubleVec1d & newtimes)
1163 {
1164 assert(IsValidTree());
1165 assert(!newtimes.empty());
1166
1167 bool updatefirstinvalid(startpoint == m_firstInvalid);
1168
1169 // First retime and remove all the branches from startpoint on...
1170 vector<Branch_ptr> newbranches;
1171 Branchiter branch(startpoint);
1172 DoubleVec1d::const_iterator newtime(newtimes.begin());
1173
1174 for( ; newtime != newtimes.end(); ++newtime)
1175 {
1176 assert(branch != m_timeList.EndBranch());
1177 // No special handling for recombinations needed; assumed handled in calling code that sets up "newtimes".
1178 (*branch)->m_eventTime = *newtime;
1179 Branch_ptr newbr(*branch);
1180 branch = m_timeList.NextBody(branch);
1181 m_timeList.Remove(newbr);
1182 newbranches.push_back(newbr);
1183 }
1184
1185 // Now put them back in, setting datalikelihood update flags as needed.
1186 vector<Branch_ptr>::iterator newbranch;
1187 for(newbranch = newbranches.begin(); newbranch != newbranches.end(); ++newbranch)
1188 {
1189 Branchiter brit(m_timeList.Collate(*newbranch));
1190 if (newbranch == newbranches.begin() && updatefirstinvalid)
1191 {
1192 m_firstInvalid = brit;
1193 }
1194 m_timeList.SetUpdateDLs(*newbranch);
1195 }
1196
1197 assert(IsValidTree());
1198
1199 } // SetNewTimesFrom
1200
1201 //------------------------------------------------------------------------------------
1202 // Method to, if necessary, shrink interval lengths, usually at the root, in an attempt to prevent explosions in the
1203 // next maximization. See the extensive comment in ChainManager::DoChain(), which is where this method is called.
1204 //
1205 // TEMPERATURE arg used only for debugging.
1206
GroomForGrowth(const DoubleVec1d & thetas,const DoubleVec1d & growths,double temperature)1207 bool Tree::GroomForGrowth(const DoubleVec1d & thetas, const DoubleVec1d & growths, double temperature)
1208 {
1209 if (thetas.empty() || growths.empty() || thetas.size() != growths.size())
1210 throw implementation_error("Tree::GroomForGrowth() received an invalid theta and/or growth vector");
1211
1212 bool treeWasModified = false;
1213 double starttime(0.0), cumulativeShift(0.0), prevIntervalLength(DBL_BIG);
1214 Branchiter brit;
1215 BranchBuffer lineages(registry.GetDataPack());
1216 for (brit = m_timeList.FirstTip(); brit != m_timeList.EndBranch(); brit = m_timeList.NextTip(brit))
1217 lineages.UpdateBranchCounts((*brit)->m_partitions, true);
1218
1219 // Iterate through the timelist, looking for fatal intervals.
1220 for (brit = m_timeList.FirstBody(); brit != m_timeList.EndBranch(); brit = m_timeList.NextBody(brit))
1221 {
1222 Branch_ptr pBranch = *brit;
1223 unsigned long nchildren(0UL), worst_pop(0UL), npops = thetas.size();
1224 double smallestExpectedLength(DBL_BIG); // Smallest E[dt] among the pops.
1225 double expectedLength;
1226 LongVec1d xpartitions = lineages.GetBranchXParts();
1227 string msg;
1228
1229 // Once we shrink an interval, we need to shift all intervals below it by this amount.
1230 // The amount accumulates.
1231 if (cumulativeShift > 0.0)
1232 {
1233 pBranch->m_eventTime -= cumulativeShift;
1234 pBranch->SetUpdateDL();
1235 }
1236
1237 // Determine how many child branches we will need to update.
1238 switch(pBranch->Event())
1239 {
1240 case btypeCoal:
1241 nchildren = static_cast<unsigned long>(NELEM);
1242 break;
1243 case btypeMig:
1244 case btypeDivMig:
1245 case btypeDisease:
1246 case btypeEpoch:
1247 nchildren = 1UL;
1248 break;
1249 case btypeRec:
1250 nchildren = 0UL; // Not literally true; this is a "magic value".
1251 break;
1252 case btypeBase:
1253 case btypeTip:
1254 {
1255 // MDEBUG We assume that we won't find a Tip; this will need to be fixed when serial sampling
1256 // is added. Unknown event type; we don't know how many children or parents this event type has.
1257 assert(false);
1258 string msg = " Unrecognized branch type (\""
1259 +ToString(pBranch->Event())
1260 + "\") received by Tree::Groom().";
1261 throw implementation_error(msg);
1262 break;
1263 }
1264 }
1265 double intervalLength = pBranch->m_eventTime - starttime;
1266 if (intervalLength <= 0.0)
1267 {
1268 if (0.0 == intervalLength && btypeRec == pBranch->Event() &&
1269 prevIntervalLength > 0.0)
1270 {
1271 // This is perfectly okay--rec branches come in pairs w/equal timestamps.
1272 lineages.UpdateBranchCounts(pBranch->m_partitions, true);
1273 // Recombination: The opposite of a coalescence--it has one child and two parents.
1274 // Each parent branch comes through separately and gets updated. We need to avoid
1275 // updating their shared child twice, so we choose to update the child when the
1276 // "left" parent comes through. (nchildren was set to the "magic value" of 0.)
1277 if (pBranch->Child(0)->Parent(0) == pBranch)
1278 {
1279 lineages.UpdateBranchCounts(pBranch->Child(0)->m_partitions, false);
1280 }
1281 starttime = pBranch->m_eventTime;
1282 prevIntervalLength = intervalLength;
1283 continue;
1284 }
1285 if (0.0 == intervalLength && btypeEpoch == pBranch->Event())
1286 {
1287 lineages.UpdateBranchCounts(pBranch->m_partitions, true);
1288 lineages.UpdateBranchCounts(pBranch->Child(0)->m_partitions, false);
1289 starttime = pBranch->m_eventTime;
1290 prevIntervalLength = intervalLength;
1291 continue;
1292 }
1293 string msg = "Tree::GroomForGrowth(), encountered a time interval length of ";
1294 msg += ToString(intervalLength) + ", followed by an event of type ";
1295 msg += ToString(pBranch->Event()) + ".";
1296 throw impossible_error(msg);
1297 }
1298
1299 // Calculate and store E[dt] for each population. Calculate E[dt] for coalescence + growth,
1300 // regardless of which event occurs at the bottom of the interval, because coalescence
1301 // + growth has the most extreme/sensitive values for lnWait() and DlnWait().
1302 for (unsigned long pop = 0; pop < npops; pop++)
1303 {
1304 const unsigned long k = xpartitions[pop];
1305 if (k > 1 && growths[pop] > 0.0)
1306 {
1307 expectedLength =
1308 (1.0 / growths[pop])
1309 *
1310 ExpE1(SafeProductWithExp(k * (k - 1) / (growths[pop] * thetas[pop]), growths[pop] * starttime));
1311 }
1312 else
1313 {
1314 continue;
1315 }
1316
1317 if (expectedLength < smallestExpectedLength)
1318 {
1319 smallestExpectedLength = expectedLength;
1320 worst_pop = pop; // For debugging message only.
1321 }
1322 }
1323
1324 // Note: Calling an interval fatal if it's 100 times its new expectation value is handwaving. So far,
1325 // 100 seems to be OK, but this value might need to be reduced to 20 or 10 if growth is very large.
1326 if (smallestExpectedLength < DBL_BIG && intervalLength >= 100.0 * smallestExpectedLength)
1327 {
1328 msg = "\nTemperature " + ToString(temperature) + ", population "
1329 + ToString(worst_pop) + ", k = " + ToString(xpartitions[worst_pop])
1330 + ", interval length is " + ToString(intervalLength)
1331 + ", expected value is " + ToString(smallestExpectedLength)
1332 + ". Interval ends with a " + ToString(pBranch->Event())
1333 + " event. Time stamps are " + ToString(starttime) + " and "
1334 + ToString(pBranch->m_eventTime) + "; theta = " + ToString(thetas[worst_pop])
1335 + " and g = " + ToString(growths[worst_pop]) + ". ";
1336
1337 // Shrink this overlong interval to a value that is between 1/2 and 3/2 times the
1338 // expectation value under the new parameter values for the "worst" population.
1339 double factor = 0.5 + m_randomSource->Float();
1340 double newIntervalLength = factor * smallestExpectedLength;
1341
1342 if (newIntervalLength < starttime * numeric_limits<double>::epsilon())
1343 {
1344 msg += "Tried to shrink this interval to a length of "
1345 + ToString(newIntervalLength) + ", but this is smaller than the "
1346 + "minimum length of "
1347 + ToString(starttime * numeric_limits<double>::epsilon())
1348 + " that can be added to " + ToString(starttime)
1349 + " without being lost to rounding error. Giving up and "
1350 + "copying the cold tree into the tree for temperature "
1351 + ToString(temperature) + "....\n";
1352 registry.GetRunReport().ReportDebug(msg);
1353 return false;
1354 }
1355
1356 pBranch->m_eventTime = starttime + newIntervalLength;
1357 pBranch->SetUpdateDL();
1358 treeWasModified = true;
1359
1360 msg += "Shrinking this interval to a length of "
1361 + ToString(newIntervalLength) + ", with a new ending time "
1362 + "stamp of " + ToString(pBranch->m_eventTime) + ".\n";
1363 registry.GetRunReport().ReportDebug(msg);
1364
1365 // Once we shrink an interval, we need to shift all intervals below it to reflect this.
1366 cumulativeShift += intervalLength - newIntervalLength;
1367
1368 }
1369
1370 // Prepare for next loop iteration (next branch in the timelist).
1371 for (unsigned long i = 0; i < nchildren; i++)
1372 lineages.UpdateBranchCounts(pBranch->Child(i)->m_partitions, false);
1373 lineages.UpdateBranchCounts(pBranch->m_partitions, true);
1374
1375 // Recombination: The opposite of a coalescence--it has one child and two parents.
1376 // Each parent branch comes through separately and gets updated. We need to avoid updating
1377 // their shared child twice, so we choose to update the child when the "left" parent comes
1378 // through. (nchildren was set to the "magic value" of 0.)
1379 if (btypeRec == pBranch->Event() && pBranch->Child(0)->Parent(0) == pBranch)
1380 {
1381 lineages.UpdateBranchCounts(pBranch->Child(0)->m_partitions, false);
1382 }
1383
1384 starttime = pBranch->m_eventTime;
1385 prevIntervalLength = intervalLength; // Used to detect successive zero-lengths, if any.
1386 }
1387
1388 if (treeWasModified)
1389 CalculateDataLikes();
1390
1391 return true;
1392 }
1393
1394 //------------------------------------------------------------------------------------
1395 // Method to, if necessary, shrink interval lengths, usually at the root, in an attempt to prevent
1396 // explosions in the next maximization. See the extensive comment in ChainManager::DoChain(),
1397 // which is where this method is called.
1398
GroomForLogisticSelection(const DoubleVec1d & thetas,double s,double temperature)1399 bool Tree::GroomForLogisticSelection(const DoubleVec1d & thetas,
1400 double s, // logistic selection coefficient
1401 double temperature) // temperature used only for debugging
1402
1403 {
1404 if (0.0 == s)
1405 return true;
1406 if (2 != thetas.size() || 0.0 == thetas[0] || 0.0 == thetas[1])
1407 {
1408 string msg = "Tree:GroomForLogisticSelection(), received invalid Theta vector.";
1409 throw implementation_error(msg);
1410 }
1411
1412 bool treeWasModified = false;
1413 double starttime(0.0), cumulativeShift(0.0), shrinkageFactor(1.0);
1414 double max_allowed_starttime(DBL_BIG);
1415 Branchiter brit, firstBadInterval, lastBadInterval;
1416 string msg;
1417
1418 if (s > 0.0)
1419 max_allowed_starttime = (EXPMAX - log(thetas[1])) / s;
1420 else
1421 max_allowed_starttime = (EXPMAX - log(thetas[0])) / (-s);
1422
1423 firstBadInterval = lastBadInterval = m_timeList.EndBranch();
1424
1425 // Iterate through the timelist, looking for fatal intervals.
1426 for (brit = m_timeList.FirstBody(); brit != m_timeList.EndBranch();
1427 brit = m_timeList.NextBody(brit))
1428 {
1429 if ((*brit)->m_eventTime >= max_allowed_starttime)
1430 {
1431 Branchiter brit2 = brit;
1432 firstBadInterval = brit;
1433 while (brit2 != m_timeList.EndBranch())
1434 {
1435 lastBadInterval = brit2;
1436 brit2 = m_timeList.NextBody(brit2);
1437 }
1438 break;
1439 }
1440 starttime = (*brit)->m_eventTime;
1441 }
1442
1443 if (firstBadInterval == m_timeList.EndBranch())
1444 return true; // nothing to shrink
1445
1446 if ((*lastBadInterval)->m_eventTime - starttime <= 0.0)
1447 {
1448 msg += "Tree::GroomForLogisticSelection(), encountered an interval ";
1449 msg += "(or sum of intervals) of length ";
1450 msg += ToString((*lastBadInterval)->m_eventTime - starttime) + ".";
1451 throw implementation_error(msg);
1452 }
1453
1454 shrinkageFactor = (max_allowed_starttime - starttime) /
1455 ((*lastBadInterval)->m_eventTime - starttime);
1456
1457 if (shrinkageFactor <= 0.0)
1458 {
1459 msg += "Tree::GroomForLogisticSelection(), computed an invalid shrinkage ";
1460 msg += "factor of " + ToString(shrinkageFactor) + ".";
1461 throw implementation_error(msg);
1462 }
1463
1464 // If there are any fatal intervals, iterate through these, and shrink them.
1465 for (brit = firstBadInterval; brit != m_timeList.EndBranch();
1466 brit = m_timeList.NextBody(brit))
1467 {
1468 Branch_ptr pBranch = *brit;
1469
1470 // Once we shrink an interval, we need to shift all intervals below it by this amount.
1471 // The amount accumulates.
1472 if (cumulativeShift > 0.0)
1473 {
1474 pBranch->m_eventTime -= cumulativeShift;
1475 pBranch->SetUpdateDL();
1476 }
1477
1478 double intervalLength = pBranch->m_eventTime - starttime;
1479 pBranch->m_eventTime = starttime + intervalLength * shrinkageFactor;
1480 pBranch->SetUpdateDL();
1481 treeWasModified = true;
1482
1483 // Once we shrink an interval, we need to shift all intervals below it to reflect this.
1484 cumulativeShift += intervalLength - intervalLength * shrinkageFactor;
1485
1486 starttime = pBranch->m_eventTime;
1487 }
1488
1489 if (treeWasModified)
1490 CalculateDataLikes();
1491
1492 return true;
1493 }
1494
1495 //------------------------------------------------------------------------------------
1496
Clone() const1497 Tree * PlainTree::Clone() const
1498 {
1499 return new PlainTree(*this, false);
1500 } // PlainTree::Clone
1501
1502 //------------------------------------------------------------------------------------
1503
MakeStump() const1504 Tree * PlainTree::MakeStump() const
1505 {
1506 return new PlainTree(*this, true);
1507 } // PlainTree stump maker
1508
1509 //------------------------------------------------------------------------------------
1510
CalculateDataLikes()1511 void PlainTree::CalculateDataLikes()
1512 {
1513 m_overallDL = CalculateDataLikesForFixedLoci();
1514 m_timeList.ClearUpdateDLs(); // Reset the updating flags.
1515 }
1516
1517 //------------------------------------------------------------------------------------
1518
Break(Branch_ptr pBranch)1519 void PlainTree::Break(Branch_ptr pBranch)
1520 {
1521 Branch_ptr pParent = pBranch->Parent(0);
1522
1523 if (pParent->CanRemove(pBranch))
1524 {
1525 Break(pParent); // Recursion
1526 m_timeList.Remove(pParent);
1527
1528 pParent = pBranch->Parent(1);
1529 if (pParent)
1530 {
1531 Break(pParent); // Recursion
1532 m_timeList.Remove(pParent);
1533 }
1534 }
1535 }
1536
1537 //____________________________________________________________________________________
1538