Context Navigation

tree.c

Visit:

Last change on this file was 10416, checked in by aboeckma, 11 years ago
Added mr bayes (no menu yet)
File size: 183.3 KB

Line
1	/*
2	* MrBayes 3
3	*
4	* (c) 2002-2010
5	*
6	* John P. Huelsenbeck
7	* Dept. Integrative Biology
8	* University of California, Berkeley
9	* Berkeley, CA 94720-3140
10	* johnh@berkeley.edu
11	*
12	* Fredrik Ronquist
13	* Swedish Museum of Natural History
14	* Box 50007
15	* SE-10405 Stockholm, SWEDEN
16	* fredrik.ronquist@nrm.se
17	*
18	* With important contributions by
19	*
20	* Paul van der Mark (paulvdm@sc.fsu.edu)
21	* Maxim Teslenko (maxim.teslenko@nrm.se)
22	*
23	* and by many users (run 'acknowledgements' to see more info)
24	*
25	* This program is free software; you can redistribute it and/or
26	* modify it under the terms of the GNU General Public License
27	* as published by the Free Software Foundation; either version 2
28	* of the License, or (at your option) any later version.
29	*
30	* This program is distributed in the hope that it will be useful,
31	* but WITHOUT ANY WARRANTY; without even the implied warranty of
32	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
33	* GNU General Public License for more details (www.gnu.org).
34	*
35	*/
36
37	#include <stdio.h>
38	#include <stdlib.h>
39	#include <time.h>
40	#include <math.h>
41	#include <string.h>
42	#include <stdarg.h>
43	#include <assert.h>
44	#include "bayes.h"
45	#include "best.h"
46	#include "mb.h"
47	#include "mcmc.h"
48	#include "mbmath.h"
49	#include "model.h"
50	#include "globals.h"
51	#include "tree.h"
52	#include "utils.h"
53	#include "command.h"
54
55	const char* const svnRevisionTreeC="$Rev: 498 $"; /* Revision keyword which is expended/updated by svn on each commit/update*/
56
57	/* local prototypes */
58	void DatedNodeDepths (TreeNode p, MrBFlt nodeDepths, int *index);
59	void DatedNodes (TreeNode p, TreeNode datedTips, int index);
60	void MarkDatedSubtreeNodes (TreeNode *p);
61	int NConstrainedTips (TreeNode *p);
62	int NDatedTips (TreeNode *p);
63	void PrintNode (char *s, int len, TreeNode *p, int isRooted);
64	void ResetPolyNode (PolyNode *p);
65	void ResetTreeNode (TreeNode *p);
66	void SetNodeDepths (Tree *t);
67
68	/* local global variable */
69	char noLabel[] = "";
70
71
72	/* AddToTreeList: Add tree at end of tree list */
73	int AddToTreeList (TreeList treeList, Tree tree)
74
75	{
76
77	TreeListElement listElement = (TreeListElement ) SafeCalloc (1, sizeof(TreeListElement));
78	if (!listElement)
79	return (ERROR);
80
81	listElement->order = (int *) SafeCalloc (tree->nIntNodes-1, sizeof(int));
82	if (!listElement->order)
83	return (ERROR);
84	listElement->next = NULL;
85
86	if (treeList->last == NULL)
87	treeList->last = treeList->first = listElement;
88	else
89	{
90	treeList->last->next = listElement;
91	treeList->last = listElement;
92	}
93
94	if (tree->isRooted)
95	StoreRTopology (tree, listElement->order);
96	else
97	StoreUTopology (tree, listElement->order);
98
99	return (NO_ERROR);
100	}
101
102
103
104
105
106	/* AllocatePolyTree: Allocate memory space for a polytomous tree */
107	PolyTree *AllocatePolyTree (int numTaxa)
108	{
109	int i;
110	PolyTree *pt;
111
112	pt = (PolyTree *) SafeCalloc (1, sizeof (PolyTree));
113	if (!pt)
114	return (NULL);
115
116	pt->memNodes = 2*numTaxa;
117	pt->nodes = (PolyNode ) SafeCalloc (2numTaxa, sizeof(PolyNode));
118	pt->allDownPass = (PolyNode *) SafeCalloc (3numTaxa, sizeof (PolyNode *));
119	pt->intDownPass = pt->allDownPass + 2*numTaxa;
120	if (pt->nodes == NULL \|\| pt->allDownPass == NULL)
121	{
122	free (pt->nodes);
123	free (pt->allDownPass);
124	free (pt);
125	return (NULL);
126	}
127
128	/* initialize nodes and set index and memoryIndex */
129	for (i=0; i<2*numTaxa; i++)
130	{
131	ResetPolyNode(&pt->nodes[i]);
132	pt->nodes[i].memoryIndex = i;
133	pt->nodes[i].index = i;
134	}
135
136	/* initialize tree properties */
137	pt->nNodes = pt->nIntNodes = 0;
138	pt->root = NULL;
139	pt->brlensDef = NO;
140	pt->isRooted = NO;
141	pt->isClock = NO;
142	pt->isCalibrated = NO;
143	pt->isRelaxed = NO;
144	pt->clockRate = 0.0;
145	strcpy(pt->name,"");
146
147	/* initialize bitsets */
148	pt->bitsets = NULL;
149
150	/* initialize relaxed clock parameters */
151	pt->nESets = 0;
152	pt->nEvents = NULL;
153	pt->position = NULL;
154	pt->rateMult = NULL;
155	pt->eSetName = NULL;
156
157	pt->nBSets = 0;
158	pt->effectiveBrLen = NULL;
159	pt->bSetName = NULL;
160
161	/* initialize population size set parameters */
162	pt->popSizeSet = NO;
163	pt->popSize = NULL;
164	pt->popSizeSetName = NULL;
165
166	return (pt);
167	}
168
169
170
171
172
173	/* AllocatePolyTreeRelClockParams: Allocate space for relaxed clock parameters */
174	int AllocatePolyTreeRelClockParams (PolyTree *pt, int nBSets, int nESets)
175	{
176	int i;
177
178	/* free previous clock params if any */
179	FreePolyTreeRelClockParams (pt);
180
181	/* set number of sets */
182	pt->nBSets = nBSets;
183	pt->nESets = nESets;
184
185	/* we do not allocate space for the actual names here; these will be NULL pointers */
186
187	/* take care of branch params */
188	if (pt->nBSets > 0)
189	{
190	pt->bSetName = (char *) SafeCalloc (pt->nBSets, sizeof (char ));
191	pt->effectiveBrLen = (MrBFlt *) SafeCalloc (pt->nBSets, sizeof (MrBFlt ));
192	for (i=0; i<pt->nBSets; i++)
193	pt->effectiveBrLen[i] = (MrBFlt *) SafeCalloc (pt->memNodes, sizeof(MrBFlt));
194	}
195
196	/* take care of breakpoint params */
197	if (pt->nESets > 0)
198	{
199	pt->eSetName = (char *) SafeCalloc (pt->nESets, sizeof(char ));
200	pt->nEvents = (int *) SafeCalloc (pt->nESets, sizeof(int ));
201	pt->position = (MrBFlt *) SafeCalloc (pt->nESets, sizeof(MrBFlt ));
202	pt->rateMult = (MrBFlt *) SafeCalloc (pt->nESets, sizeof(MrBFlt ));
203	for (i=0; i<pt->nESets; i++)
204	{
205	pt->nEvents[i] = (int *) SafeCalloc (pt->memNodes, sizeof(int));
206	pt->position[i] = (MrBFlt *) SafeCalloc (pt->memNodes, sizeof(MrBFlt ));
207	pt->rateMult[i] = (MrBFlt *) SafeCalloc (pt->memNodes, sizeof(MrBFlt ));
208	}
209	}
210
211	return (NO_ERROR);
212	}
213
214
215
216
217
218	/* AllocatePolyTreePartitions: Allocate space for and set partitions for polytomous tree */
219	int AllocatePolyTreePartitions (PolyTree *pt)
220	{
221	int i, nLongsNeeded, numTaxa;
222
223	/* get some handy numbers */
224	numTaxa = pt->memNodes/2;
225	nLongsNeeded = (numTaxa -1) / nBitsInALong + 1;
226
227	/* allocate space */
228	pt->bitsets = (SafeLong ) SafeRealloc ((void )pt->bitsets, pt->memNodesnLongsNeededsizeof(SafeLong));
229	if (pt->bitsets == NULL)
230	return (ERROR);
231	for (i=0; i<pt->memNodes*nLongsNeeded; i++)
232	pt->bitsets[i] = 0;
233
234	/* set node partition pointers */
235	for (i=0; i<pt->memNodes; i++)
236	pt->nodes[i].partition = pt->bitsets + i*nLongsNeeded;
237
238	/* clear and set partitions; if the tree is empty, nothing is set */
239	ResetPolyTreePartitions(pt);
240
241	return (NO_ERROR);
242	}
243
244
245
246
247
248	/* AllocateTree: Allocate memory space for a tree (unrooted or rooted) */
249	Tree *AllocateTree (int numTaxa)
250	{
251	int i;
252	Tree *t;
253
254	t = (Tree *) SafeCalloc (1, sizeof (Tree));
255	if (t == NULL)
256	return NULL;
257
258	/* initialize basic tree properties */
259	t->memNodes = 2*numTaxa;
260	strcpy (t->name, "");
261
262	t->isRooted = NO;
263	t->isClock = NO;
264
265	t->checkConstraints = NO;
266	t->nConstraints = 0;
267	t->nLocks = 0;
268	t->isCalibrated = NO;
269	t->nNodes = t->nIntNodes = 0;
270	t->nRelParts = 0;
271	t->relParts = NULL;
272
273	/* initialize pointers */
274	t->bitsets = NULL;
275	t->flags = NULL;
276	t->constraints = NULL;
277
278	/* allocate and initialize nodes and node arrays (enough for both rooted and unrooted trees) */
279	t->nNodes = 0;
280	t->nIntNodes = 0;
281	if ((t->nodes = (TreeNode ) SafeCalloc (2numTaxa, sizeof (TreeNode))) == NULL)
282	{
283	free (t);
284	return NULL;
285	}
286	if ((t->allDownPass = (TreeNode *) SafeCalloc (3numTaxa, sizeof (TreeNode *))) == NULL)
287	{
288	free (t->nodes);
289	free (t);
290	return NULL;
291	}
292	t->intDownPass = t->allDownPass + t->memNodes;
293
294	/* initialize nodes and set index and memoryIndex */
295	for (i=0; i<t->memNodes; i++)
296	{
297	ResetTreeNode(&t->nodes[i]);
298	t->nodes[i].memoryIndex = i;
299	t->nodes[i].index = i;
300	}
301
302	return t;
303	}
304
305
306
307
308
309	/* AllocateFixedTree: Allocate memory space for a fixed unrooted or rooted tree */
310	Tree *AllocateFixedTree (int numTaxa, int isRooted)
311	{
312	int i;
313	Tree *t;
314
315	t = (Tree *) SafeCalloc (1, sizeof (Tree));
316	if (t == NULL)
317	return NULL;
318
319	/* initialize basic tree properties */
320	if (isRooted == YES)
321	t->memNodes = 2*numTaxa;
322	else
323	t->memNodes = 2*numTaxa - 2;
324	strcpy (t->name, "");
325
326	t->isRooted = isRooted;
327	t->isClock = NO;
328
329	t->checkConstraints = NO;
330	t->nConstraints = 0;
331	t->nLocks = 0;
332	t->isCalibrated = NO;
333	t->nNodes = t->nIntNodes = 0;
334	t->nRelParts = 0;
335	t->relParts = NULL;
336
337	/* initialize pointers */
338	t->bitsets = NULL;
339	t->flags = NULL;
340	t->constraints = NULL;
341
342	/* allocate and initialize nodes and node arrays (enough for both rooted and unrooted trees) */
343	if (t->isRooted)
344	{
345	t->nNodes = 2*numTaxa;
346	t->nIntNodes = numTaxa - 1;
347	}
348	else
349	{
350	t->nNodes = 2*numTaxa - 2;
351	t->nIntNodes = numTaxa - 2;
352	}
353	if ((t->nodes = (TreeNode *) SafeCalloc (t->nNodes, sizeof (TreeNode))) == NULL)
354	{
355	free (t);
356	return NULL;
357	}
358	if ((t->allDownPass = (TreeNode *) SafeCalloc (t->nNodes + t->nIntNodes, sizeof (TreeNode ))) == NULL)
359	{
360	free (t->nodes);
361	free (t);
362	return NULL;
363	}
364	t->intDownPass = t->allDownPass + t->nNodes;
365
366	/* initialize nodes and set index and memoryIndex */
367	for (i=0; i<t->memNodes; i++)
368	{
369	ResetTreeNode(&t->nodes[i]);
370	t->nodes[i].memoryIndex = i;
371	t->nodes[i].index = i;
372	}
373
374	return t;
375	}
376
377
378
379
380
381	/* AllocateTreePartitions: Allocate space for and set partitions for tree */
382	int AllocateTreePartitions (Tree *t)
383	{
384	int i, nLongsNeeded, numTaxa;
385	TreeNode *p;
386
387	/* get some handy numbers */
388	if (t->isRooted == YES)
389	numTaxa = t->nNodes - t->nIntNodes - 1;
390	else
391	numTaxa = t->nNodes - t->nIntNodes;
392	nLongsNeeded = (numTaxa - 1) / nBitsInALong + 1;
393
394	/* reallocate space */
395	t->bitsets = (SafeLong ) SafeRealloc ((void ) t->bitsets, (size_t)(t->nNodesnLongsNeededsizeof(SafeLong)));
396	if (!t->bitsets)
397	return (ERROR);
398
399	/* clear bit fields */
400	for (i=0; i<t->nNodes*nLongsNeeded; i++)
401	t->bitsets[0] = 0;
402
403	/* set node pointers to bit fields */
404	for (i=0; i<t->nNodes; i++)
405	{
406	p = t->allDownPass[i];
407	p->partition = t->bitsets + i*nLongsNeeded;
408	}
409
410	/* set partition specifiers for terminals */
411	ResetTreePartitions(t);
412
413	return (NO_ERROR);
414	}
415
416
417
418
419
420	int AreTopologiesSame (Tree t1, Tree t2)
421
422	{
423	int i, j, k, nLongsNeeded, nTaxa, revertBits;
424	SafeLong bitsets, mask;
425	TreeNode p, q;
426
427	if (t1->nNodes != t2->nNodes)
428	return (NO);
429	if (t1->nIntNodes != t2->nIntNodes)
430	return (NO);
431
432	if (t1->isRooted == NO && (t1->root->index != t2->root->index))
433	revertBits = YES;
434	else
435	revertBits = NO;
436
437	if (t1->isRooted == YES)
438	nTaxa = t1->nNodes - t1->nIntNodes - 1;
439	else
440	nTaxa = t1->nNodes - t1->nIntNodes;
441
442	/* allocate space */
443	nLongsNeeded = (nTaxa - 1) / nBitsInALong + 1;
444	bitsets = (SafeLong ) SafeCalloc (4nLongsNeeded*nTaxa+nLongsNeeded, sizeof(SafeLong));
445	mask = bitsets + 4nLongsNeedednTaxa;
446
447	/* set mask */
448	for (i=0; i<nTaxa; i++)
449	SetBit(i, mask);
450
451	/* set partition pointers */
452	for (i=0; i<t1->nNodes; i++)
453	{
454	p = t1->allDownPass[i];
455	q = t2->allDownPass[i];
456	p->partition = bitsets + i2nLongsNeeded;
457	q->partition = p->partition + nLongsNeeded;
458	}
459
460	/* set terminal partitions */
461	for (i=0; i<t1->nNodes; i++)
462	{
463	p = t1->allDownPass[i];
464	q = t2->allDownPass[i];
465	if (p->right == NULL)
466	SetBit (p->index, p->partition);
467	if (q->right == NULL)
468	SetBit (q->index, q->partition);
469	}
470
471	/* set internal partitions */
472	for (i=0; i<t1->nIntNodes; i++)
473	{
474	p = t1->intDownPass[i];
475	q = t2->intDownPass[i];
476	for (j=0; j<nLongsNeeded; j++)
477	{
478	p->partition[j] = p->left->partition[j] \| p->right->partition[j];
479	q->partition[j] = q->left->partition[j] \| q->right->partition[j];
480	}
481	}
482
483	/* check for congruence */
484	for (i=0; i<t1->nIntNodes; i++)
485	{
486	p = t1->intDownPass[i];
487	if (t1->isRooted == NO && IsBitSet(t2->root->index,p->partition))
488	FlipBits(p->partition,nLongsNeeded, mask);
489	for (j=0; j<t2->nIntNodes; j++)
490	{
491	q = t2->intDownPass[j];
492	for (k=0; k<nLongsNeeded; k++)
493	{
494	if (p->partition[k] != q->partition[k])
495	break;
496	}
497	if (k == nLongsNeeded)
498	break;
499	}
500	if (j == t2->nIntNodes)
501	{
502	free (bitsets);
503	return (NO);
504	}
505	}
506
507	free (bitsets);
508	return (YES);
509	}
510
511
512
513
514
515	int AreTreesSame (Tree t1, Tree t2)
516
517	{
518	int i, j, k, nLongsNeeded, nTaxa, revertBits;
519	SafeLong bitsets, mask;
520	TreeNode p, q;
521
522	extern void ShowNodes(TreeNode*, int, int);
523
524	if (t1->nNodes != t2->nNodes)
525	return (NO);
526	if (t1->nIntNodes != t2->nIntNodes)
527	return (NO);
528
529	if (t1->isRooted == NO && (t1->root->index != t2->root->index))
530	revertBits = YES;
531	else
532	revertBits = NO;
533
534	if (t1->isRooted == YES)
535	nTaxa = t1->nNodes - t1->nIntNodes - 1;
536	else
537	nTaxa = t1->nNodes - t1->nIntNodes;
538
539	/* allocate space */
540	nLongsNeeded = (nTaxa - 1) / nBitsInALong + 1;
541	bitsets = (SafeLong ) SafeCalloc (4nLongsNeeded*nTaxa+nLongsNeeded, sizeof(SafeLong));
542	mask = bitsets + 4nLongsNeedednTaxa;
543
544	/* set mask */
545	for (i=0; i<nTaxa; i++)
546	SetBit(i, mask);
547
548	/* set partition pointers */
549	for (i=0; i<t1->nNodes; i++)
550	{
551	p = t1->allDownPass[i];
552	q = t2->allDownPass[i];
553	p->partition = bitsets + i2nLongsNeeded;
554	q->partition = p->partition + nLongsNeeded;
555	}
556
557	/* set terminal partitions */
558	for (i=0; i<t1->nNodes; i++)
559	{
560	p = t1->allDownPass[i];
561	q = t2->allDownPass[i];
562	if (p->right == NULL)
563	SetBit (p->index, p->partition);
564	if (q->right == NULL)
565	SetBit (q->index, q->partition);
566	}
567
568	/* set internal partitions */
569	for (i=0; i<t1->nIntNodes; i++)
570	{
571	p = t1->intDownPass[i];
572	q = t2->intDownPass[i];
573	for (j=0; j<nLongsNeeded; j++)
574	{
575	p->partition[j] = p->left->partition[j] \| p->right->partition[j];
576	q->partition[j] = q->left->partition[j] \| q->right->partition[j];
577	}
578	}
579
580	/* check for congruence */
581	for (i=0; i<t1->nNodes; i++)
582	{
583	p = t1->allDownPass[i];
584	if (p->anc == NULL && t1->isRooted == YES)
585	continue;
586	if (t1->isRooted == NO && IsBitSet(t2->root->index,p->partition))
587	FlipBits(p->partition,nLongsNeeded, mask);
588	for (j=0; j<t2->nNodes; j++)
589	{
590	q = t2->allDownPass[j];
591	for (k=0; k<nLongsNeeded; k++)
592	{
593	if (p->partition[k] != q->partition[k])
594	break;
595	}
596	if (k == nLongsNeeded && AreDoublesEqual (p->length, q->length, 0.000001) == YES)
597	break;
598	else if (k == nLongsNeeded)
599	{
600	free (bitsets);
601	return (NO);
602	}
603	}
604	if (j == t2->nNodes)
605	{
606	free (bitsets);
607	return (NO);
608	}
609	}
610	free (bitsets);
611	return (YES);
612	}
613
614
615
616
617
618	/*----------------------------------------------------------------
619	\|
620	\| BuildConstraintTree: Build constraint tree. The tree t is
621	\| needed only to hold information about constraints and
622	\| included taxa.
623	\|
624	----------------------------------------------------------------*/
625	int BuildConstraintTree (Tree t, PolyTree pt, char **localTaxonNames)
626
627	{
628
629	int i, j, k, constraintId, nLongsNeeded, nextNode;
630	SafeLong constraintPartition, mask;
631	PolyNode pp, qq, rr, ss, *tt;
632
633	pt->isRooted = t->isRooted;
634
635	nLongsNeeded = (numLocalTaxa - 1) / nBitsInALong + 1;
636	constraintPartition = (SafeLong ) SafeCalloc (2nLongsNeeded, sizeof(SafeLong));
637	if (!constraintPartition)
638	{
639	MrBayesPrint ("%s Problems allocating constraintPartition in BuildConstraintTree", spacer);
640	return (ERROR);
641	}
642	mask = constraintPartition + nLongsNeeded;
643
644	/* calculate mask (needed to take care of unused bits when flipping partitions) */
645	for (i=0; i<numLocalTaxa; i++)
646	SetBit (i, mask);
647
648	/* reset all nodes */
649	for (i=0; i<2*numLocalTaxa; i++)
650	{
651	pp = &pt->nodes[i];
652	pp->isDated = NO;
653	pp->calibration = NULL;
654	pp->age = -1.0;
655	pp->isLocked = NO;
656	pp->lockID = -1;
657	pp->index = i;
658	}
659
660	/* build a bush */
661	pt->root = &pt->nodes[numLocalTaxa];
662	for (i=0; i<numLocalTaxa; i++)
663	{
664	pp = &pt->nodes[i];
665	pp->index = i;
666	pp->left = NULL;
667	if (i == numLocalTaxa - 1)
668	pp->sib = NULL;
669	else
670	pp->sib = &pt->nodes[i+1];
671	pp->anc = pt->root;
672	}
673	pp = pt->root;
674	pp->left = &pt->nodes[0];
675	pp->anc = pp->sib = NULL;
676	pt->nNodes = numLocalTaxa + 1;
677	pt->nIntNodes = 1;
678
679	/* make sure the outgroup is the right-most node */
680	pt->nodes[localOutGroup].index = numLocalTaxa - 1;
681	pt->nodes[numLocalTaxa - 1].index = localOutGroup;
682
683	/* allocate and set partition specifiers in bush */
684	GetPolyDownPass(pt);
685	AllocatePolyTreePartitions(pt);
686
687	/* set terminal taxon labels */
688	for (i=0; i<pt->nNodes; i++)
689	{
690	pp = pt->allDownPass[i];
691	if (pp->index < numLocalTaxa)
692	strcpy (pp->label, localTaxonNames[pp->index]);
693	}
694
695	/* resolve the bush according to constraints */
696	/* for now, satisfy all constraints */
697	/* for now, bail out if constraints are not compatible */
698	/* Eventually, we might want to be build a parsimony (WAB) or compatibility (WIB) matrix and
699	draw a starting tree from the universe according to the score of the tree. A simple way of accomplishing
700	approximately this is to use sequential addition, with probabilities in each step determined
701	by the parsimony or compatibility score of the different possibilities. */
702	nextNode = numLocalTaxa + 1;
703	for (constraintId=0; constraintId<numDefinedConstraints; constraintId++)
704	{
705	if (t->constraints[constraintId] == NO \|\| definedConstraintsType[constraintId] != HARD )
706	continue;
707
708	/* initialize bits in partition to add; get rid of deleted taxa in the process */
709	ClearBits(constraintPartition, nLongsNeeded);
710	for (i=j=0; i<numTaxa; i++)
711	{
712	if (taxaInfo[i].isDeleted == YES)
713	continue;
714	if (IsBitSet(i, definedConstraint[constraintId]) == YES)
715	SetBit(j, constraintPartition);
716	j++;
717	}
718	assert (j == numLocalTaxa);
719
720	/* make sure outgroup is outside constrained partition if the tree is unrooted */
721	if (t->isRooted == NO && IsBitSet(localOutGroup, constraintPartition))
722	FlipBits(constraintPartition, nLongsNeeded, mask);
723
724	/* check that partition should be included */
725	k = NumBits(constraintPartition, nLongsNeeded);
726	if (k == 0)
727	{
728	MrBayesPrint ("%s WARNING: Constraint '%s' refers only to deleted taxa\n", spacer, constraintNames[constraintId]);
729	MrBayesPrint ("%s and will be disregarded\n", spacer);
730	t->constraints[constraintId] = NO;
731	t->nLocks--;
732	continue;
733	}
734	if (k == 1)
735	{
736	MrBayesPrint ("%s WARNING: Constraint '%s' refers to a single tip and\n", spacer, constraintNames[constraintId]);
737	MrBayesPrint ("%s will be disregarded\n", spacer);
738	t->constraints[constraintId] = NO;
739	t->nLocks--;
740	continue;
741	}
742
743	/* check if root in rooted tree (we allow this to enable inference of ancestral states) */
744	if (k == numLocalTaxa && t->isRooted == YES)
745	{
746	pt->root->isLocked = YES;
747	pt->root->lockID = constraintId;
748	continue;
749	}
750
751	/* check if interior root in unrooted tree (we allow this to enable inference of ancestral states) */
752	if ((k == numLocalTaxa - 1 \|\| k == numLocalTaxa) && t->isRooted == NO)
753	{
754	pt->root->isLocked = YES;
755	pt->root->lockID = constraintId;
756	continue;
757	}
758
759	/* find first included terminal */
760	k = FirstTaxonInPartition (constraintPartition, nLongsNeeded);
761	for (i=0; pt->nodes[i].index != k; i++)
762	;
763	pp = &pt->nodes[i];
764
765	/* go down until node is not included in constraint */
766	do {
767	qq = pp;
768	pp = pp->anc;
769	} while (IsPartNested(pp->partition, constraintPartition, nLongsNeeded));
770
771	/* check that the node has not yet been included */
772	for (i=0; i<nLongsNeeded; i++)
773	{
774	if (qq->partition[i] != constraintPartition[i])
775	break;
776	}
777	if (i==nLongsNeeded)
778	{
779	MrBayesPrint ("%s WARNING: Constraint '%s' is a duplicate of another constraint\n", spacer, constraintNames[constraintId]);
780	MrBayesPrint ("%s and will be ignored\n", spacer);
781	t->nLocks--;
782	t->constraints[i] = NO;
783	continue;
784	}
785
786	/* create a new node */
787	tt = &pt->nodes[nextNode++];
788	tt->anc = pp;
789	tt->isLocked = YES;
790	tt->lockID = constraintId;
791	for (i=0; i<nLongsNeeded; i++)
792	tt->partition[i] = constraintPartition[i];
793	pt->nIntNodes++;
794	pt->nNodes++;
795
796	/* sort descendant nodes in two connected groups: included and excluded */
797	/* if there is a descendant that overlaps (incompatible) then return error */
798	rr = ss = NULL;
799	qq = pp->left;
800	do {
801	if (IsPartNested(qq->partition, constraintPartition, nLongsNeeded))
802	{
803	if (ss != NULL)
804	ss->sib = qq;
805	else
806	tt->left = qq;
807	ss = qq;
808	qq->anc = tt;
809	}
810	else if (IsPartCompatible(qq->partition, constraintPartition, nLongsNeeded))
811	{
812	if (rr != NULL)
813	rr->sib = qq;
814	else
815	tt->sib = qq;
816	rr = qq;
817	}
818	else
819	{
820	free (constraintPartition);
821	return (ERROR);
822	}
823	qq = qq->sib;
824	} while (qq != NULL);
825	pp->left = tt;
826	rr->sib = ss->sib = NULL;
827	}
828
829	/* relabel interior nodes */
830	GetPolyDownPass(pt);
831	for (i=0; i<pt->nIntNodes; i++)
832	pt->intDownPass[i]->index = i + numLocalTaxa;
833
834	/* exit */
835	free (constraintPartition);
836	FreePolyTreePartitions(pt);
837	return (NO_ERROR);
838	}
839
840
841
842
843
844	/*----------------------------------------------
845	\|
846	\| BuildRandomRTopology: Builds a random rooted
847	\| topology. Will set indices in t->nodes
848	\| such that they are from 0 to n-1 for the n tips
849	\| and from n to 2n-2 for the n-1 interior
850	\| nodes. Last is root. Does not touch labels
851	\| of tips.
852	\|
853	----------------------------------------------*/
854	int BuildRandomRTopology (Tree t, SafeLong seed)
855	{
856	int i, j, nTips;
857	TreeNode p, q, *r;
858
859	nTips = t->nNodes - t->nIntNodes - 1;
860
861	for (i=0; i<t->nNodes; i++)
862	{
863	p = &t->nodes[i];
864	p->index = i;
865	p->left = p->right = p->anc = NULL;
866	}
867
868	/* connect the first two tip nodes */
869	q = &t->nodes[0];
870	r = &t->nodes[1];
871	p = &t->nodes[nTips];
872	q->anc = r->anc = p;
873	p->left = q;
874	p->right = r;
875	q = &t->nodes[2*nTips-1];
876	p->anc = q;
877	q->left = p;
878
879	for (i=2; i<nTips; i++)
880	{
881	q = &t->nodes[i];
882	r = &t->nodes[i-2+nTips+1];
883	q->anc = r;
884	r->left = q;
885	j = (int) (RandomNumber(seed) * (2 * i - 1));
886	if (j < i)
887	p = &t->nodes[j];
888	else
889	p = &t->nodes[j-i + nTips];
890	r->right = p;
891	r->anc = p->anc;
892	if (p->anc != NULL)
893	{
894	if (p->anc->left == p)
895	p->anc->left = r;
896	else
897	p->anc->right = r;
898	}
899	p->anc = r;
900	}
901
902	/* set root and get downpass */
903	t->root = &t->nodes[2*nTips-1];
904	GetDownPass(t);
905
906	/* relabel interior nodes */
907	for (i=0; i<t->nIntNodes; i++)
908	t->intDownPass[i]->index = i+nTips;
909
910	return (NO_ERROR);
911	}
912
913
914
915
916
917	/*----------------------------------------------
918	\|
919	\| BuildRandomUTopology: Builds a random unrooted
920	\| topology. Assumes that indices are set
921	\| in t->nodes from 0 to n-1 for the n tips
922	\| and from n to 2n-3 for the n-2 interior
923	\| nodes. Move the calculation root after
924	\| this routine to get the right root.
925	\|
926	----------------------------------------------*/
927	int BuildRandomUTopology (Tree t, SafeLong seed)
928	{
929	int i, j, nTips;
930	TreeNode p, q, *r;
931
932	nTips = t->nNodes - t->nIntNodes;
933
934	for (i=0; i<t->nNodes; i++)
935	{
936	p = &t->nodes[i];
937	p->index = i;
938	p->left = p->right = p->anc = NULL;
939	}
940
941	/* connect the first three nodes, assuming 0 is calc root */
942	q = &t->nodes[1];
943	r = &t->nodes[2];
944	p = &t->nodes[nTips];
945	q->anc = r->anc = p;
946	p->left = q;
947	p->right = r;
948	q = &t->nodes[0];
949	p->anc = q;
950	q->left = p;
951
952	for (i=3; i<nTips; i++)
953	{
954	q = &t->nodes[i];
955	r = &t->nodes[i - 3 + nTips + 1];
956	q->anc = r;
957	r->left = q;
958	j = (int) (RandomNumber(seed) * (2 * i - 3));
959	if (j < i - 1)
960	p = &t->nodes[j+1];
961	else
962	p = &t->nodes[j+1-i + nTips];
963	r->right = p;
964	r->anc = p->anc;
965	if (p->anc->left == p)
966	p->anc->left = r;
967	else
968	p->anc->right = r;
969	p->anc = r;
970	}
971
972	t->root = &t->nodes[0];
973
974	/* get downpass */
975	GetDownPass(t);
976
977	/* relabel interior nodes */
978	for (i=0; i<t->nIntNodes; i++)
979	t->intDownPass[i]->index = i+nTips;
980
981	return (NO_ERROR);
982	}
983
984
985
986
987
988	/*----------------------------------------------------------------
989	\|
990	\| CheckConstraints: Check that tree complies with constraints
991	\|
992	----------------------------------------------------------------*/
993	int CheckConstraints (Tree *t)
994
995	{
996
997	int a, i, j, k, nLongsNeeded;
998	SafeLong constraintPartition, mask;
999	TreeNode *p=NULL;
1000
1001	if (t->checkConstraints == NO)
1002	return (NO_ERROR);
1003
1004	/* allocate space */
1005	nLongsNeeded = (numLocalTaxa - 1) / nBitsInALong + 1;
1006	constraintPartition = (SafeLong ) SafeCalloc (2nLongsNeeded, sizeof(SafeLong));
1007	if (!constraintPartition)
1008	{
1009	MrBayesPrint ("%s Problems allocating constraintPartition in CheckConstraints", spacer);
1010	return (ERROR);
1011	}
1012	mask = constraintPartition + nLongsNeeded;
1013
1014	/* set mask (needed to reset unused bits when flipping partitions) */
1015	for (i=0; i<numLocalTaxa; i++)
1016	SetBit (i, mask);
1017
1018	if (AllocateTreePartitions(t) == ERROR)
1019	{
1020	MrBayesPrint ("%s Problems allocating tree partitions in CheckConstraints", spacer);
1021	return (ERROR);
1022	}
1023
1024	for (a=0; a<numDefinedConstraints; a++)
1025	{
1026	if (t->constraints[a] == NO \|\| definedConstraintsType[a] != HARD)
1027	continue;
1028
1029	/* set bits in partition to check */
1030	ClearBits(constraintPartition, nLongsNeeded);
1031	for (j=k=0; j<numTaxa; j++)
1032	{
1033	if (taxaInfo[j].isDeleted == YES)
1034	continue;
1035	if (IsBitSet(j, definedConstraint[a]) == YES)
1036	SetBit(k, constraintPartition);
1037	k++;
1038	}
1039
1040	/* make sure outgroup is outside constrained partition if unrooted tree */
1041	if (t->isRooted == NO && IsBitSet(localOutGroup, constraintPartition))
1042	FlipBits(constraintPartition, nLongsNeeded, mask);
1043
1044	/* find the locked node */
1045	for (i=j=0; i<t->nNodes; i++)
1046	{
1047	if (t->allDownPass[i]->isLocked == YES && t->allDownPass[i]->lockID == a)
1048	{
1049	p = t->allDownPass[i];
1050	j++;
1051	}
1052	}
1053
1054	if (j != 1)
1055	{
1056	MrBayesPrint ("%s Tree has %d locks with id %d identifying constraint '%s'\n", spacer, j, a, constraintNames[a]);
1057	free (constraintPartition);
1058	FreeTreePartitions(t);
1059	return (ERROR);
1060	}
1061
1062	/* check that locked node is correct */
1063	for (i=0; i<nLongsNeeded; i++)
1064	{
1065	if (p->partition[i] != constraintPartition[i])
1066	{
1067	MrBayesPrint ("%s Lock %d is set for the wrong node [this is a bug]\n", spacer, a);
1068	free (constraintPartition);
1069	FreeTreePartitions(t);
1070	return (ERROR);
1071	}
1072	}
1073	}
1074
1075	FreeTreePartitions (t);
1076	free (constraintPartition);
1077	return (NO_ERROR);
1078	}
1079
1080
1081
1082
1083	/*----------------------------------------------------------------
1084	\|
1085	\| CheckSetConstraints: Check and set tree constraints
1086	\|
1087	----------------------------------------------------------------*/
1088	int CheckSetConstraints (Tree *t)
1089
1090	{
1091
1092	int a, i, j, k, nLongsNeeded, foundIt;
1093	SafeLong constraintPartition, mask;
1094	TreeNode *p;
1095
1096	if (t->checkConstraints == NO)
1097	return (NO_ERROR);
1098
1099	/* reset all existing locks, if any */
1100	t->nLocks=0;
1101	for (i=0; i<t->nNodes; i++)
1102	{
1103	p = t->allDownPass[i];
1104	p->isLocked = NO;
1105	p->lockID = -1;
1106	if (p->left != NULL)
1107	{
1108	p->calibration = NULL;
1109	p->isDated = NO;
1110	p->age = -1.0;
1111	}
1112	}
1113
1114	/* allocate space */
1115	if (AllocateTreePartitions (t) == ERROR)
1116	{
1117	MrBayesPrint ("%s Problems allocating tree bitsets", spacer);
1118	return ERROR;
1119	}
1120
1121	nLongsNeeded = ((numLocalTaxa - 1) / nBitsInALong) + 1;
1122	constraintPartition = (SafeLong ) SafeCalloc (2nLongsNeeded, sizeof(SafeLong));
1123	if (!constraintPartition)
1124	{
1125	MrBayesPrint ("%s Problems allocating constraintPartition", spacer);
1126	FreeTreePartitions(t);
1127	return ERROR;
1128	}
1129	mask = constraintPartition + nLongsNeeded;
1130
1131	/* set mask (needed to take care of unused bits when flipping partitions) */
1132	for (i=0; i<numLocalTaxa; i++)
1133	SetBit (i, mask);
1134
1135	for (a=0; a<numDefinedConstraints; a++)
1136	{
1137	if (modelParams[t->relParts[0]].activeConstraints[a] == NO \|\| definedConstraintsType[a] != HARD )
1138	continue;
1139
1140	/* set bits in partition to add */
1141	ClearBits(constraintPartition, nLongsNeeded);
1142	for (i=j=0; i<numTaxa; i++)
1143	{
1144	if (taxaInfo[i].isDeleted == YES)
1145	continue;
1146	if (IsBitSet(i, definedConstraint[a]) == YES)
1147	SetBit(j, constraintPartition);
1148	j++;
1149	}
1150
1151	k = NumBits(constraintPartition, nLongsNeeded);
1152	if (k == 0 \|\| k == 1)
1153	continue;
1154
1155	/* make sure outgroup is outside constrained partition (marked 0) */
1156	if (t->isRooted == NO && IsBitSet(localOutGroup, constraintPartition) == YES)
1157	FlipBits(constraintPartition, nLongsNeeded, mask);
1158
1159	/* find the node that should be locked */
1160	foundIt = NO;
1161	for (i=0; i<t->nIntNodes; i++)
1162	{
1163	p = t->intDownPass[i];
1164	for (j=0; j<nLongsNeeded; j++)
1165	{
1166	if (p->partition[j] != constraintPartition[j])
1167	break;
1168	}
1169
1170	if (j == nLongsNeeded)
1171	{
1172	foundIt = YES;
1173	p->isLocked = YES;
1174	p->lockID = a;
1175	if (nodeCalibration[a].prior != unconstrained)
1176	{
1177	p->isDated = YES;
1178	p->calibration = &nodeCalibration[a];
1179	}
1180	t->nLocks++;
1181	break;
1182	}
1183	}
1184
1185	if (foundIt == NO)
1186	{
1187	MrBayesPrint ("%s Tree breaks constraint '%s'\n", spacer, constraintNames[a]);
1188	FreeTreePartitions (t);
1189	free (constraintPartition);
1190	return (ERROR);
1191	}
1192	}
1193
1194	/* exit */
1195	FreeTreePartitions(t);
1196	free (constraintPartition);
1197	return (NO_ERROR);
1198	}
1199
1200
1201
1202
1203
1204	/*-----------------------------------------------------------------------
1205	\|
1206	\| ColorClusters: Recursive function to color the clusters in a tree by
1207	\| assigning numbers to them in their variable x
1208	\|
1209	------------------------------------------------------------------------*/
1210	void ColorClusters (TreeNode p, int index)
1211	{
1212	if (p!=NULL)
1213	{
1214	if (p->isLocked == YES \|\| p->anc == NULL \|\| p->anc->anc == NULL)
1215	p->x = (++(*index));
1216	else
1217	p->x = p->anc->x;
1218	ColorClusters(p->left, index);
1219	ColorClusters(p->right, index);
1220	}
1221	}
1222
1223
1224
1225
1226
1227	/* CopyPolyNodes: Copies everything except pointers and memoryIndex */
1228	void CopyPolyNodes (PolyNode p, PolyNode q, int nLongsNeeded)
1229	{
1230	p->index = q->index;
1231	p->mark = q->mark;
1232	p->length = q->length;
1233	p->x = q->x;
1234	p->y = q->y;
1235	p->isDated = q->isDated;
1236	p->calibration = q->calibration;
1237	p->age = q->age;
1238	p->isLocked = q->isLocked;
1239	p->lockID = q->lockID;
1240	strcpy (p->label, q->label);
1241	if( nLongsNeeded!=0 )
1242	{
1243	assert(p->partition);
1244	assert(q->partition);
1245	memcpy(p->partition,q->partition, nLongsNeeded*sizeof(SafeLong));
1246	}
1247	p->support = q->support;
1248	p->f = q->f;
1249	}
1250
1251
1252
1253
1254
1255	void CopySubtreeToTree (Tree subtree, Tree t)
1256
1257	{
1258
1259	int i, /j,/ k;
1260	TreeNode p, q=NULL, *r;
1261
1262	for (i=/j=/0; i<subtree->nNodes - 1; i++)
1263	{
1264	p = subtree->allDownPass[i];
1265
1266	for (k=0; k<t->nNodes; k++)
1267	{
1268	q = t->allDownPass[k];
1269	if (q->index == p->index)
1270	break;
1271	}
1272	q->length = p->length;
1273	q->marked = YES;
1274	if (p->left != NULL && p->right != NULL)
1275	{
1276	for (k=0; k<t->nNodes; k++)
1277	{
1278	r = t->allDownPass[k];
1279	if (r->index == p->left->index)
1280	{
1281	q->left = r;
1282	r->anc = q;
1283	}
1284	else if (r->index == p->right->index)
1285	{
1286	q->right = r;
1287	r->anc = q;
1288	}
1289	}
1290	}
1291	}
1292
1293	p = subtree->root;
1294
1295	for (k=0; k<t->nNodes; k++)
1296	{
1297	q = t->allDownPass[k];
1298	if (q->index == p->index)
1299	break;
1300	}
1301
1302	if (q->left->marked == YES)
1303	{
1304	for (k=0; k<t->nIntNodes; k++)
1305	{
1306	r = t->intDownPass[k];
1307	if (r->index == p->left->index)
1308	{
1309	q->left = r;
1310	r->anc = q;
1311	}
1312	}
1313	}
1314	else if (q->right->marked == YES)
1315	{
1316	for (k=0; k<t->nIntNodes; k++)
1317	{
1318	r = t->intDownPass[k];
1319	if (r->index == p->left->index)
1320	{
1321	q->right = r;
1322	r->anc = q;
1323	}
1324	}
1325	}
1326	}
1327
1328
1329
1330
1331
1332	/*-----------------------------------------------------------------
1333	\|
1334	\| CopyToPolyTreeFromPolyTree: copies second tree to first tree
1335	\|
1336	-----------------------------------------------------------------*/
1337	int CopyToPolyTreeFromPolyTree (PolyTree to, PolyTree from)
1338	{
1339	int i, j, k, nLongsNeeded;
1340	PolyNode p, q;
1341
1342	/* check we have enough memory */
1343	assert (to->memNodes >= from->nNodes);
1344	if( from->bitsets==NULL \|\| to->bitsets==NULL )
1345	{
1346	nLongsNeeded=0;
1347	}
1348	else
1349	{
1350	assert (to->memNodes >= from->memNodes);/Otherwise partition length woould not be long enough for nodes in "to" /
1351	nLongsNeeded = (from->memNodes/2 - 1) / nBitsInALong + 1;
1352	}
1353
1354
1355	/* copy nodes */
1356	for (i=0; i<from->nNodes; i++)
1357	{
1358	/* copy pointers */
1359	p = from->nodes + i;
1360	q = to->nodes + i;
1361
1362	if (p->anc != NULL)
1363	q->anc = to->nodes + p->anc->memoryIndex;
1364	else
1365	{
1366	q->anc = NULL;
1367	to->root = q;
1368	}
1369
1370	if (p->left != NULL)
1371	q->left = to->nodes + p->left->memoryIndex;
1372	else
1373	q->left = NULL;
1374
1375	if (p->sib != NULL)
1376	q->sib = to->nodes + p->sib->memoryIndex;
1377	else
1378	q->sib = NULL;
1379
1380	/* Copy everything else except memoryIndex */
1381	CopyPolyNodes (q, p, nLongsNeeded);
1382	}
1383
1384	/* fill node arrays */
1385	/* copy tree properties */
1386	to->nNodes = from->nNodes;
1387	to->nIntNodes = from->nIntNodes;
1388	to->isRooted = from->isRooted;
1389	to->isClock = from->isClock;
1390	to->isRelaxed = from->isRelaxed;
1391	to->clockRate = from->clockRate;
1392	strcpy (to->name, from->name);
1393
1394	GetPolyDownPass (to);
1395
1396	/* copy partitions */
1397	if (from->bitsets)
1398	{
1399	if (!to->bitsets)
1400	AllocatePolyTreePartitions(to);
1401	else
1402	ResetPolyTreePartitions(to);
1403	}
1404
1405	/* copy relaxed clock parameters */
1406	FreePolyTreeRelClockParams (to);
1407
1408	if (from->nBSets + from->nESets > 0)
1409	AllocatePolyTreeRelClockParams (to, from->nBSets, from->nESets);
1410
1411	for (i=0; i<to->nBSets; i++)
1412	{
1413	to->bSetName[i] = (char *) SafeCalloc (strlen(from->bSetName[i])+2, sizeof(char));
1414	strcpy (to->bSetName[i], from->bSetName[i]);
1415	for (j=0; j<from->nNodes; j++)
1416	to->effectiveBrLen[i][j] = from->effectiveBrLen[i][j];
1417	}
1418
1419	for (i=0; i<to->nESets; i++)
1420	{
1421	to->eSetName[i] = (char *) SafeCalloc (strlen(from->eSetName[i])+2, sizeof(char));
1422	strcpy (to->eSetName[i], from->eSetName[i]);
1423	for (j=0; j<from->nNodes; j++)
1424	{
1425	to->nEvents[i][j] = from->nEvents[i][j];
1426	if (to->nEvents[i][j] > 0)
1427	{
1428	to->position[i][j] = (MrBFlt *) SafeCalloc (to->nEvents[i][j], sizeof (MrBFlt));
1429	to->rateMult[i][j] = (MrBFlt *) SafeCalloc (to->nEvents[i][j], sizeof (MrBFlt));
1430	for (k=0; k<to->nEvents[i][j]; k++)
1431	{
1432	to->position[i][j][k] = from->position[i][j][k];
1433	to->rateMult[i][j][k] = from->rateMult[i][j][k];
1434	}
1435	}
1436	}
1437	}
1438
1439	/* copy population size parameters */
1440	FreePolyTreePopSizeParams(to);
1441	to->popSizeSet = from->popSizeSet;
1442	if (to->popSizeSet == YES)
1443	{
1444	to->popSize = (MrBFlt *) SafeCalloc (to->nNodes-1, sizeof(MrBFlt));
1445	for (i=0; i<to->nNodes-1; i++)
1446	to->popSize[i] = from->popSize[i];
1447	to->popSizeSetName = (char *) SafeCalloc (strlen(from->popSizeSetName) + 1, sizeof(char));
1448	strcpy (to->popSizeSetName, from->popSizeSetName);
1449	}
1450
1451	return (NO_ERROR);
1452
1453	}
1454
1455
1456
1457
1458
1459	/*-----------------------------------------------------------------
1460	\|
1461	\| CopyToSpeciesTreeFromPolyTree: copies second tree (polytomous) to
1462	\| first tree, which is a species tree. The species tree needs to
1463	\| be allocated enough space first to hold the resulting tree.
1464	\|
1465	-----------------------------------------------------------------*/
1466	int CopyToSpeciesTreeFromPolyTree (Tree to, PolyTree from)
1467
1468	{
1469
1470	int i;
1471	PolyNode *p;
1472	TreeNode q, q1;
1473	#if !defined (NDEBUG)
1474	int j;
1475	#endif
1476
1477	/* make sure assumptions are correct */
1478	assert (from->isRooted == YES);
1479	assert (from->isClock == YES);
1480	assert (from->nNodes - from->nIntNodes == numSpecies);
1481	assert (to->memNodes == 2*numSpecies);
1482	assert (to->nIntNodes == from->nIntNodes);
1483	assert (to->nNodes == from->nNodes + 1);
1484
1485	/* make sure indices are set correctly for from nodes */
1486	#if !defined (NDEBUG)
1487	for (i=0; i<from->nNodes; i++)
1488	{
1489	for (j=0; j<from->nNodes; j++)
1490	{
1491	p = from->allDownPass[j];
1492	if (p->index == i)
1493	break;
1494	}
1495	assert (j != from->nNodes);
1496	assert (!(p->left == NULL && p->index >= numSpecies));
1497	}
1498	#endif
1499
1500	/* copy nodes */
1501	for (i=0; i<from->nNodes; i++)
1502	{
1503	/* copy pointers */
1504	p = from->allDownPass[i];
1505	q = to->nodes + p->index;
1506
1507	if (p->anc != NULL)
1508	q->anc = to->nodes + p->anc->index;
1509	else
1510	q->anc = NULL;
1511
1512	if (p->left != NULL)
1513	q->left = to->nodes + p->left->index;
1514	else
1515	q->left = NULL;
1516
1517	if (p->left != NULL)
1518	q->right = to->nodes + p->left->sib->index;
1519	else
1520	q->right = NULL;
1521
1522	q->nodeDepth = p->depth;
1523	q->age = p->age;
1524	q->length = p->length;
1525	q->index = p->index;
1526	if (q->index < numSpecies)
1527	q->label = speciesNameSets[speciespartitionNum].names[q->index];
1528	else
1529	q->label = noLabel;
1530	}
1531
1532	/* fix root */
1533	p = from->root;
1534	q = to->nodes + p->index;
1535	q1 = to->nodes + from->nNodes; /* get the 'extra' root node that polytomous trees do not use */
1536	q->anc = q1;
1537	q1->index = from->nNodes;
1538	q1->left = q;
1539	q1->right = q1->anc = NULL;
1540	q1->isLocked = NO;
1541	q1->lockID = -1;
1542	q1->isDated = NO;
1543	q1->calibration = NULL;
1544	q1->age = -1.0;
1545	to->root = q1;
1546
1547	/* get downpass */
1548	GetDownPass (to);
1549
1550	/* a user tree might not come with node depths set */
1551	if (to->root->left->nodeDepth == 0.0)
1552	SetNodeDepths(to);
1553
1554	/* set partitions */
1555	if (to->bitsets)
1556	ResetTreePartitions(to);
1557
1558	return (NO_ERROR);
1559	}
1560
1561
1562
1563
1564
1565	/*-----------------------------------------------------------------
1566	\|
1567	\| CopyToTreeFromPolyTree: copies second tree (polytomous) to first
1568	\| tree (used to initialize constrained starting trees, e.g.).
1569	\| An unrooted source tree will be rooted on outgroup
1570	\| An unrooted source tree that needs to be copied to
1571	\| a rooted target tree will be randomly rooted on a node below
1572	\| all defined constraints. The to tree needs to be allocated
1573	\| enough space first to hold the resulting tree.
1574	\|
1575	-----------------------------------------------------------------*/
1576	int CopyToTreeFromPolyTree (Tree to, PolyTree from)
1577
1578	{
1579
1580	int i, j, nNodesNeeded;
1581	PolyNode *p;
1582	TreeNode q, q1;
1583
1584	/* refuse to arbitrarily root an input tree */
1585	assert (!(from->isRooted == NO && to->isRooted == YES));
1586	if ( (from->isRooted == NO) && (to->isRooted == YES) )
1587	{
1588	MrBayesPrint ("%s Fail to copy trees due to difference in rootedness of source and destination. \n", spacer);
1589	return (ERROR);
1590	}
1591	/* calculate space needed */
1592	if (from->isRooted == YES && to->isRooted == YES)
1593	nNodesNeeded = from->nNodes + 1;
1594	else if (from->isRooted == YES && to->isRooted == NO)
1595	nNodesNeeded = from->nNodes;
1596	else /* if (from->isRooted = NO && to->isRooted == NO) */
1597	nNodesNeeded = from->nNodes;
1598
1599	/* make sure assumptions are in order */
1600	assert (to->memNodes >= nNodesNeeded);
1601	assert (localOutGroup >= 0 && localOutGroup < numLocalTaxa);
1602	assert (numLocalTaxa == from->nNodes - from->nIntNodes);
1603	assert (!(from->isRooted == YES && from->nNodes != 2*from->nIntNodes + 1));
1604	assert (!(from->isRooted == NO && from->nNodes != 2*from->nIntNodes + 2));
1605
1606	/* make sure indices are set correctly for from nodes */
1607	for (i=0; i<from->nNodes; i++)
1608	{
1609	for (j=0; j<from->nNodes; j++)
1610	{
1611	p = from->allDownPass[j];
1612	if (p->index == i)
1613	break;
1614	}
1615	assert (j != from->nNodes);
1616	assert (!(p->left == NULL && p->index >= numLocalTaxa));
1617	}
1618
1619	/* deal with root */
1620	if (to->isRooted == NO && from->isRooted == YES)
1621	Deroot(from);
1622
1623	/* make sure calculation root is set correctly */
1624	if (to->isRooted == NO && MovePolyCalculationRoot (from, localOutGroup) == ERROR)
1625	return ERROR;
1626
1627	/* copy nodes */
1628	for (i=0; i<from->nNodes; i++)
1629	{
1630	/* copy pointers */
1631	p = from->allDownPass[i];
1632	q = to->nodes + p->index;
1633
1634	if (p->anc != NULL)
1635	q->anc = to->nodes + p->anc->index;
1636	else
1637	q->anc = NULL;
1638
1639	if (p->left != NULL)
1640	q->left = to->nodes + p->left->index;
1641	else
1642	q->left = NULL;
1643
1644	if (p->left != NULL)
1645	q->right = to->nodes + p->left->sib->index;
1646	else
1647	q->right = NULL;
1648
1649	q->isLocked = p->isLocked;
1650	q->lockID = p->lockID;
1651	q->isDated = p->isDated;
1652	q->calibration = p->calibration;
1653	q->age = p->age;
1654	q->nodeDepth = p->depth;
1655	q->length = p->length;
1656	q->index = p->index;
1657	if (q->index < numLocalTaxa)
1658	q->label = localTaxonNames[q->index];
1659	else
1660	q->label = noLabel;
1661	}
1662
1663	/* fix root */
1664	if (to->isRooted == NO)
1665	{
1666	p = from->root;
1667	q = to->nodes + p->index;
1668	q->anc = to->root = to->nodes + p->left->sib->sib->index;
1669	q->length = to->root->length;
1670	to->root->length = 0.0;
1671	to->root->left = q;
1672	to->root->right = to->root->anc = NULL;
1673	}
1674	else
1675	{
1676	p = from->root;
1677	q = to->nodes + p->index;
1678	q1 = to->nodes + from->nNodes; /* get the 'extra' root node that polytomous trees do not use */
1679	q->anc = q1;
1680	q1->index = from->nNodes;
1681	q1->left = q;
1682	q1->right = q1->anc = NULL;
1683	q1->isLocked = NO;
1684	q1->lockID = -1;
1685	q1->isDated = NO;
1686	q1->calibration = NULL;
1687	q1->age = -1.0;
1688	to->root = q1;
1689	}
1690
1691	/* get downpass */
1692	GetDownPass (to);
1693
1694	/* set node depths */
1695	if (to->isRooted == YES && to->root->left->nodeDepth == 0.0)
1696	SetNodeDepths(to);
1697
1698	/* set partitions */
1699	if (to->bitsets)
1700	ResetTreePartitions(to);
1701
1702	/* relaxed clock parameters are not stored in binary trees but in separate parameters */
1703
1704	return (NO_ERROR);
1705	}
1706
1707
1708
1709
1710
1711	/*-----------------------------------------------------------------
1712	\|
1713	\| CopyToTreeFromTree: copies second tree to first tree
1714	\| (used to initialize brlen sets for same topology)
1715	\| Note: partition information of nodes are not copied if
1716	\| either "from" or "to" tree does not have bitsets allocated
1717	\|
1718	-----------------------------------------------------------------*/
1719	int CopyToTreeFromTree (Tree to, Tree from)
1720
1721	{
1722
1723	int i, numTaxa, nLongsNeeded;
1724	TreeNode p, q;
1725
1726	numTaxa = from->nNodes - from->nIntNodes - (from->isRooted == YES ? 1 : 0);
1727	nLongsNeeded = (numTaxa - 1) / nBitsInALong + 1;
1728	if( from->bitsets==NULL \|\| to->bitsets==NULL )
1729	nLongsNeeded=0;
1730
1731	/* check that there is enough memory */
1732	assert (to->memNodes >= from->nNodes);
1733
1734	/* copy nodes (use index of p as memoryIndex for q) */
1735	for (i=0; i<from->nNodes; i++)
1736	{
1737	/* copy pointers */
1738	p = from->nodes + i;
1739	q = to->nodes + p->index;
1740
1741	if (p->anc != NULL)
1742	q->anc = to->nodes + p->anc->index;
1743	else
1744	{
1745	q->anc = NULL;
1746	to->root = q;
1747	}
1748
1749	if (p->left != NULL)
1750	q->left = to->nodes + p->left->index;
1751	else
1752	q->left = NULL;
1753
1754	if (p->right != NULL)
1755	q->right = to->nodes + p->right->index;
1756	else
1757	q->right = NULL;
1758
1759	CopyTreeNodes (q, p, nLongsNeeded);
1760	}
1761
1762	/* create new node arrays */
1763	to->nNodes = from->nNodes;
1764	to->nIntNodes = from->nIntNodes;
1765	GetDownPass(to);
1766
1767	/* copy tree properties (these should be constant most of them) */
1768	strcpy (to->name, from->name);
1769	to->isRooted = from->isRooted;
1770	to->isClock = from->isClock;
1771	to->isCalibrated = from->isCalibrated;
1772	to->checkConstraints = from->checkConstraints;
1773	to->nConstraints = from->nConstraints;
1774	to->constraints = from->constraints;
1775	to->nLocks = from->nLocks;
1776	to->nRelParts = from->nRelParts;
1777	to->relParts = from->relParts;
1778
1779	/* copy partitions */
1780	if (from->bitsets)
1781	{
1782	if (!to->bitsets)
1783	AllocateTreePartitions(to);
1784	else
1785	ResetTreePartitions(to);
1786	}
1787
1788	return (NO_ERROR);
1789
1790	}
1791
1792
1793
1794
1795
1796	/* Copeing node q to node p. nLongsNeeded is the length of partition that the node represents. */
1797	void CopyTreeNodes (TreeNode p, TreeNode q, int nLongsNeeded)
1798	{
1799
1800	/* copies everything except pointers and memoryIndex */
1801	p->index = q->index;
1802	p->scalerNode = q->scalerNode;
1803	p->upDateCl = q->upDateCl;
1804	p->upDateTi = q->upDateTi;
1805	p->marked = q->marked;
1806	p->length = q->length;
1807	p->nodeDepth = q->nodeDepth;
1808	p->x = q->x;
1809	p->y = q->y;
1810	p->isDated = q->isDated;
1811	p->calibration = q->calibration;
1812	p->age = q->age;
1813	p->isLocked = q->isLocked;
1814	p->lockID = q->lockID;
1815	p->d = q->d;
1816	//p->partition = q->partition;//the content should be copied not the pointers. Otherwise we are risking to have segmentation faults
1817	if( nLongsNeeded!=0 )
1818	{
1819	assert(p->partition);
1820	assert(q->partition);
1821	memcpy(p->partition,q->partition, nLongsNeeded*sizeof(SafeLong));
1822	}
1823	p->label = q->label;
1824	}
1825
1826
1827
1828
1829
1830	void CopyTreeToSubtree (Tree t, Tree subtree)
1831
1832	{
1833
1834	int i, j, k;
1835	TreeNode p, q, *r;
1836
1837	for (i=j=0; i<t->nNodes; i++)
1838	{
1839	p = t->allDownPass[i];
1840	if (p->marked == NO)
1841	continue;
1842
1843	q = &subtree->nodes[j++];
1844	q->index = p->index;
1845	q->length = p->length;
1846	if (p->left == NULL \|\| p->left->marked == NO)
1847	q->left = q->right = NULL;
1848	else
1849	{
1850	for (k=0; k<j-1; k++)
1851	{
1852	r = &subtree->nodes[k];
1853	if (r->index == p->left->index)
1854	{
1855	q->left = r;
1856	r->anc = q;
1857	}
1858	else if (r->index == p->right->index)
1859	{
1860	q->right = r;
1861	r->anc = q;
1862	}
1863	}
1864	}
1865
1866	if (p->anc->marked == NO)
1867	{
1868	r = &subtree->nodes[j++];
1869	subtree->root = r;
1870	r->anc = r->right = NULL;
1871	r->left = q;
1872	q->anc = r;
1873	r->length = 0.0;
1874	r->index = p->anc->index;
1875	}
1876
1877	}
1878
1879	GetDownPass (subtree);
1880
1881	subtree->isRooted = t->isRooted;
1882	subtree->nRelParts = t->nRelParts;
1883	subtree->relParts = t->relParts;
1884	}
1885
1886
1887
1888
1889
1890	/* DatedNodeDepths: Recursive function to get node depths */
1891	void DatedNodeDepths (TreeNode p, MrBFlt nodeDepths, int *index)
1892	{
1893	if (p != NULL)
1894	{
1895	if (p->left == NULL \|\| p->isDated == YES)
1896	nodeDepths[(*index)++] = p->nodeDepth;
1897	else
1898	{
1899	DatedNodeDepths (p->left, nodeDepths, index);
1900	DatedNodeDepths (p->right, nodeDepths, index);
1901	}
1902	}
1903	}
1904
1905
1906
1907
1908
1909	/* DatedNodes: Recursive function to get dated tips or interior nodes */
1910	void DatedNodes (TreeNode p, TreeNode datedNodes, int index)
1911	{
1912	if (p != NULL)
1913	{
1914	if (p->left != NULL && p->isDated == NO)
1915	{
1916	DatedNodes (p->left, datedNodes, index);
1917	DatedNodes (p->right, datedNodes, index);
1918	}
1919	datedNodes[(*index)++] = p;
1920	}
1921	}
1922
1923
1924
1925
1926
1927	/* Deroot: Deroot a rooted polytomous tree with branch lengths */
1928	int Deroot (PolyTree *pt)
1929	{
1930	PolyNode p, q, *r, tempNode;
1931	int i;
1932
1933	p = pt->root;
1934
1935	if (p->left->sib->sib != NULL)
1936	return (ERROR); /* tree is not rooted or it is polytomous */
1937
1938	if (p != &pt->nodes[pt->nNodes-1])
1939	{
1940	q = &pt->nodes[pt->nNodes-1];
1941	/* now swap content of p and q including pointers */
1942	tempNode = *q;
1943	q = p;
1944	*p = tempNode;
1945	/* swap back memoryindex */
1946	p->memoryIndex = q->memoryIndex;
1947	q->memoryIndex = tempNode.memoryIndex;
1948	/* all pointers to q should be pointers to p */
1949	for (i=0; i<pt->nNodes; i++)
1950	{
1951	r = &pt->nodes[i];
1952	if (r->left == q)
1953	r->left = p;
1954	if (r->anc == q)
1955	r->anc = p;
1956	if (r->sib == q)
1957	r->sib = p;
1958	}
1959	/* all pointers to p should be pointers to q; all these are anc pointers from the descendants of the root */
1960	pt->root = q;
1961	for (r=q->left; r!=NULL; r=r->sib)
1962	r->anc = q;
1963	/* finally set p to the new root */
1964	p = pt->root;
1965	}
1966
1967	/* make sure the left of the old root is interior and can be used as new root */
1968	if (p->left->left == NULL)
1969	{
1970	q = p->left;
1971	r = q->sib;
1972	p->left = r;
1973	r->sib = q;
1974	q->sib = NULL;
1975	}
1976
1977	pt->root = p->left;
1978	pt->root->left->sib->sib = p->left->sib;
1979	p->left->sib->length += pt->root->length;
1980	pt->root->length = 0.0;
1981	pt->root->sib = NULL;
1982	pt->root->anc = NULL;
1983
1984	pt->nNodes--;
1985	pt->nIntNodes--;
1986
1987	GetPolyDownPass(pt);
1988
1989	return (NO_ERROR);
1990	}
1991
1992
1993
1994
1995
1996	/* EraseTreeList: Erase all trees in treeList */
1997	void EraseTreeList (TreeList *treeList)
1998	{
1999	TreeListElement *listElement;
2000	TreeListElement *previous;
2001
2002	listElement = treeList->first;
2003	if (listElement != NULL)
2004	do
2005	{
2006	free (listElement->order);
2007	previous = listElement;
2008	listElement = listElement->next;
2009	free (previous);
2010	}
2011	while (listElement != NULL);
2012
2013	treeList->first = treeList->last = NULL;
2014	}
2015
2016
2017
2018
2019	void UpdateTreeWithClockrate (Tree *t, MrBFlt clockRate)
2020	{
2021	int i;
2022	TreeNode *p;
2023
2024	if( t->fromUserTree == NO )
2025	{
2026	/Set nodeDepth/
2027	for (i=0; i<t->nNodes; i++)
2028	{
2029	p = t->allDownPass[i];
2030	p->nodeDepth = p->age * clockRate;
2031	}
2032
2033	/* calculate branch lengths */
2034	for (i=0; i<t->nNodes; i++)
2035	{
2036	p = t->allDownPass[i];
2037	if (p->anc != NULL)
2038	{
2039	if (p->anc->anc != NULL)
2040	{
2041	p->length = p->anc->nodeDepth - p->nodeDepth;
2042	}
2043	else
2044	p->length = 0.0; //not a problem for root node.
2045	}
2046	}
2047	}
2048	else
2049	{
2050	for (i=0; i<t->nNodes-1; i++)
2051	{
2052	p = t->allDownPass[i];
2053	p->age = p->nodeDepth / clockRate;
2054	}
2055	}
2056
2057	}
2058
2059
2060
2061
2062	/*----------------------------------------------------------------
2063	\|
2064	\| findAllowedClockrate: Finds the range of clock rates allowed for the tree.
2065	\|
2066	\| @param t - tree to check (IN)
2067	\| @minClockRate - adress where minimum allowed clock rate is stored (OUT)
2068	\| @maxClockRate - adress where maximum allowed clock rate is stored (OUT)
2069	\|
2070	----------------------------------------------------------------*/
2071	void findAllowedClockrate (Tree t, MrBFlt minClockRate, MrBFlt *maxClockRate )
2072	{
2073	int i;
2074	TreeNode *p;
2075	MrBFlt min, max, tmp;
2076
2077	min=0.0;
2078	max=MRBFLT_MAX;
2079
2080	*minClockRate = 2.0;
2081	*maxClockRate = 1.0;
2082
2083	if( t->fromUserTree == NO )
2084	{
2085	for (i=0; i<t->nNodes-1; i++)
2086	{
2087	p = t->allDownPass[i];
2088	if (p->anc->anc != NULL)
2089	{
2090	tmp = BRLENS_MIN/(p->anc->age - p->age);
2091	assert(tmp > 0);
2092	if( tmp > min)
2093	min = tmp;
2094
2095	tmp = BRLENS_MAX/(p->anc->age - p->age);
2096	assert(tmp > 0);
2097	if( tmp > max)
2098	max = tmp;
2099	}
2100	}
2101	*minClockRate= min;
2102	*maxClockRate= max;
2103	}
2104	else
2105	{
2106	IsCalibratedClockSatisfied (t,minClockRate,maxClockRate, 0.001);
2107	}
2108
2109	}
2110
2111
2112
2113
2114
2115	/* FreePolyTree: Free memory space for a polytomous tree (unrooted or rooted) */
2116	void FreePolyTree (PolyTree *pt)
2117	{
2118	if (pt != NULL)
2119	{
2120	FreePolyTreePartitions(pt);
2121	FreePolyTreeRelClockParams(pt);
2122	FreePolyTreePopSizeParams(pt);
2123	free (pt->allDownPass);
2124	free (pt->nodes);
2125	free (pt);
2126	}
2127	}
2128
2129
2130
2131
2132
2133	/* FreePolyTreePartitions: Free memory space for polytomous tree partitions */
2134	void FreePolyTreePartitions (PolyTree *pt)
2135	{
2136	int i;
2137	if (pt != NULL && pt->bitsets != NULL)
2138	{
2139	for (i=0; i<pt->memNodes; i++)
2140	pt->nodes[i].partition = NULL;
2141	free (pt->bitsets);
2142	pt->bitsets = NULL;
2143	}
2144	}
2145
2146
2147
2148
2149
2150	/* FreePolyTreePopSizeParams: Free population size set parameters of polytree */
2151	void FreePolyTreePopSizeParams (PolyTree *pt)
2152	{
2153	if (pt->popSizeSet == YES)
2154	{
2155	free (pt->popSize);
2156	free (pt->popSizeSetName);
2157	}
2158	pt->popSizeSet = NO;
2159	pt->popSize = NULL;
2160	pt->popSizeSetName = NULL;
2161	}
2162
2163
2164
2165
2166
2167	/* FreePolyTreeRelClockParams: Free relaxed clock parameters of polytree */
2168	void FreePolyTreeRelClockParams (PolyTree *pt)
2169	{
2170	int i, j;
2171
2172	/* free breakpoint clock parameters */
2173	for (i=0; i<pt->nESets; i++)
2174	{
2175	for (j=0; j<pt->memNodes; j++)
2176	{
2177	if (pt->nEvents[i][j] > 0)
2178	{
2179	free (pt->position[i][j]);
2180	free (pt->rateMult[i][j]);
2181	}
2182	}
2183	free (pt->eSetName[i]);
2184	free (pt->nEvents[i]);
2185	free (pt->position[i]);
2186	free (pt->rateMult[i]);
2187	}
2188	free (pt->nEvents);
2189	free (pt->position);
2190	free (pt->rateMult);
2191	free (pt->eSetName);
2192	pt->nESets = 0;
2193	pt->nEvents = NULL;
2194	pt->position = NULL;
2195	pt->rateMult = NULL;
2196	pt->eSetName = NULL;
2197
2198	/* free branch clock parameters */
2199	for (i=0; i<pt->nBSets; i++)
2200	{
2201	free (pt->bSetName[i]);
2202	free (pt->effectiveBrLen[i]);
2203	}
2204	free (pt->effectiveBrLen);
2205	free (pt->bSetName);
2206	pt->nBSets = 0;
2207	pt->effectiveBrLen = NULL;
2208	pt->bSetName = NULL;
2209	}
2210
2211
2212
2213
2214
2215	/* FreeTree: Free memory space for a tree (unrooted or rooted) */
2216	void FreeTree (Tree *t)
2217	{
2218	if (t != NULL)
2219	{
2220	free (t->bitsets);
2221	free (t->flags);
2222	free (t->allDownPass);
2223	free (t->nodes);
2224	free (t);
2225	}
2226	}
2227
2228
2229
2230
2231
2232	/* FreeTreePartitions: Free memory space for tree partitions */
2233	void FreeTreePartitions (Tree *t)
2234	{
2235	int i;
2236
2237	if (t != NULL && t->bitsets != NULL)
2238	{
2239	free (t->bitsets);
2240	t->bitsets = NULL;
2241	for (i=0; i<t->memNodes; i++)
2242	t->nodes[i].partition = NULL;
2243	}
2244	}
2245
2246
2247
2248
2249
2250	/*-------------------------------------------------------------------------------------------
2251	\|
2252	\| GetDatedNodeDepths: Get an array containing the node depths of the dated tips,
2253	\| internal or external, plus dated root
2254	\|
2255	---------------------------------------------------------------------------------------------*/
2256	void GetDatedNodeDepths (TreeNode p, MrBFlt nodeDepths)
2257	{
2258	int index = 0;
2259
2260	assert (p != NULL);
2261
2262	nodeDepths[index++] = p->nodeDepth; /* include root node depth */
2263	if (p->left != NULL)
2264	{
2265	DatedNodeDepths (p->left, nodeDepths, &index);
2266	DatedNodeDepths (p->right, nodeDepths, &index);
2267	}
2268	}
2269
2270
2271
2272
2273
2274	/*-------------------------------------------------------------------------------------------
2275	\|
2276	\| GetDatedNodes: Get an array containing the dated tips,
2277	\| internal or external, and all interior nodes in the same subtree
2278	\|
2279	---------------------------------------------------------------------------------------------*/
2280	void GetDatedNodes (TreeNode p, TreeNode *datedNodes)
2281	{
2282	int index = 0;
2283
2284	assert (p != NULL);
2285
2286	if (p->left!= NULL)
2287	{
2288	DatedNodes (p->left, datedNodes, &index);
2289	DatedNodes (p->right, datedNodes, &index);
2290	}
2291	}
2292
2293
2294
2295
2296
2297	/* get down pass for tree t (wrapper function) */
2298	void GetDownPass (Tree *t)
2299
2300	{
2301
2302	int i, j;
2303
2304	i = j = 0;
2305	GetNodeDownPass (t, t->root, &i, &j);
2306
2307	}
2308
2309
2310
2311
2312
2313	/* get the actual down pass sequences */
2314	void GetNodeDownPass (Tree t, TreeNode p, int i, int j)
2315
2316	{
2317
2318	if (p != NULL )
2319	{
2320	GetNodeDownPass (t, p->left, i, j);
2321	GetNodeDownPass (t, p->right, i, j);
2322	if (p->left != NULL && p->right != NULL && p->anc != NULL)
2323	{
2324	t->intDownPass[(*i)++] = p;
2325	t->allDownPass[(*j)++] = p;
2326	}
2327	else if (p->left == NULL && p->right == NULL && p->anc != NULL)
2328	{
2329	t->allDownPass[(*j)++] = p;
2330	}
2331	else if (p->left != NULL && p->right == NULL && p->anc == NULL)
2332	{
2333	t->allDownPass[(*j)++] = p;
2334	}
2335	}
2336
2337	}
2338
2339
2340
2341
2342
2343	/* GetPolyAges: Get PolyTree node ages */
2344	void GetPolyAges (PolyTree *t)
2345	{
2346	int i;
2347	PolyNode *p;
2348
2349	GetPolyDepths (t); /* just to make sure... */
2350
2351	for (i=0; i<t->nNodes; i++)
2352	{
2353	p = t->allDownPass[i];
2354	p->age = p->depth / t->clockRate;
2355	}
2356	}
2357
2358
2359
2360
2361
2362	/* GetPolyDepths: Get PolyTree node depths */
2363	void GetPolyDepths (PolyTree *t)
2364	{
2365	int i;
2366	MrBFlt maxDepth;
2367	PolyNode *p;
2368
2369	maxDepth = t->root->depth = 0.0;
2370
2371	for (i=t->nNodes-2; i>=0; i--)
2372	{
2373	p = t->allDownPass[i];
2374	p->depth = p->anc->depth + p->length;
2375	if (p->depth > maxDepth)
2376	maxDepth = p->depth;
2377	}
2378
2379	for (i=0; i<t->nNodes; i++)
2380	{
2381	p = t->allDownPass[i];
2382	p->depth = maxDepth - p->depth;
2383	}
2384	}
2385
2386
2387
2388
2389
2390	/* get down pass for polytomous tree t (wrapper function) */
2391	void GetPolyDownPass (PolyTree *t)
2392
2393	{
2394
2395	int i, j;
2396
2397	i = j = 0;
2398	GetPolyNodeDownPass (t, t->root, &i, &j);
2399	assert( t->nIntNodes==j );
2400
2401	}
2402
2403
2404
2405
2406
2407	/* get the actual down pass sequences for a polytomous tree */
2408	void GetPolyNodeDownPass (PolyTree t, PolyNode p, int i, int j)
2409
2410	{
2411
2412	PolyNode *q;
2413
2414	if (p->left != NULL)
2415	{
2416	for (q=p->left; q!=NULL; q=q->sib)
2417	GetPolyNodeDownPass(t, q, i, j);
2418	}
2419
2420	t->allDownPass[(*i)++] = p;
2421	if (p->left != NULL )
2422	t->intDownPass[(*j)++] = p;
2423
2424	}
2425
2426
2427
2428
2429
2430	/* GetFromTreeList: Get first tree from a tree list and remove it from the list*/
2431	int GetFromTreeList (TreeList treeList, Tree tree)
2432
2433	{
2434	TreeListElement *listElement;
2435
2436	if (treeList->first == NULL)
2437	{
2438	MrBayesPrint ("%s Tree list empty\n", spacer);
2439	return (ERROR);
2440	}
2441	if (tree->isRooted == YES)
2442	RetrieveRTopology (tree, treeList->first->order);
2443	else
2444	{
2445	RetrieveUTopology (tree, treeList->first->order);
2446	if (localOutGroup != 0)
2447	MoveCalculationRoot (tree, localOutGroup);
2448	}
2449
2450	listElement = treeList->first;
2451	treeList->first = listElement->next;
2452
2453	free (listElement->order);
2454	free (listElement);
2455
2456	return (NO_ERROR);
2457	}
2458
2459
2460
2461
2462
2463	/*------------------------------------------------------------------
2464	\|
2465	\| InitBrlens: This routine will set all branch lengths of a
2466	\| nonclock tree to the value given by 'v'.
2467	\|
2468	------------------------------------------------------------------*/
2469	int InitBrlens (Tree *t, MrBFlt v)
2470
2471	{
2472	int i;
2473	TreeNode *p;
2474
2475	for (i=0; i<t->nNodes; i++)
2476	{
2477	p = t->allDownPass[i];
2478	if (p->anc != NULL && !(t->isRooted == YES && p->anc->anc == NULL))
2479	p->length = v;
2480	else
2481	p->length = 0.0;
2482	}
2483
2484	return (NO_ERROR);
2485	}
2486
2487
2488
2489
2490
2491	/*
2492	@param levUp is the number of edges between the "node" and the most resent calibrated predecessor +1,
2493	for the calibrated nodes it should be 1
2494	@param calibrUp is the age of the most resent calibrated predecessor
2495	@return age of the node
2496	*/
2497	MrBFlt SetNodeCalibratedAge(TreeNode *node, unsigned levUp, MrBFlt calibrUp )
2498	{
2499	MrBFlt r,l;
2500
2501	if( node->age != -1.0 )
2502	{
2503	if(node->right != NULL)
2504	SetNodeCalibratedAge( node->right, 2, node->age );
2505	if(node->left != NULL)
2506	SetNodeCalibratedAge( node->left, 2, node->age );
2507	return node->age;
2508	}
2509
2510	r = SetNodeCalibratedAge( node->right, levUp+1, calibrUp );
2511	l = SetNodeCalibratedAge( node->left, levUp+1, calibrUp );
2512
2513	if( r>l )
2514	{
2515	assert( calibrUp - r > 0.0 );
2516	return node->age = r + ( calibrUp - r )/levUp;
2517	}
2518	else
2519	{
2520	assert( calibrUp - l> 0.0 );
2521	return node->age = l + ( calibrUp - l )/levUp;
2522	}
2523
2524	}
2525
2526
2527
2528
2529
2530	/*-------------------------------------------------------------------
2531	\|
2532	\| InitCalibratedBrlens: This routine will build a clock tree
2533	\| consistent with calibration constraints on terminal
2534	\| taxa and/or constrained interior nodes. At least one
2535	\| node should be calibrated.
2536	\| If not possible to build such a tree, ERROR
2537	\| is returned.
2538	\|
2539	--------------------------------------------------------------------*/
2540	int InitCalibratedBrlens (Tree t, MrBFlt clockRate, SafeLong seed)
2541
2542	{
2543
2544	int i, treeAgeFixed;
2545	TreeNode *p;
2546	Model *mp;
2547	MrBFlt treeAge;
2548
2549	#if 0
2550	printf ("Before initializing calibrated brlens\n");
2551	ShowNodes(t->root, 0, YES);
2552	#endif
2553
2554	if (t->isRooted == NO)
2555	{
2556	MrBayesPrint ("%s Tree is unrooted\n", spacer);
2557	return (ERROR);
2558	}
2559
2560	/* Check whether root is fixed indirectly */
2561	mp = &modelParams[t->relParts[0]];
2562	if (!strcmp(mp->clockPr, "Uniform") && !strcmp(mp->treeAgePr,"Fixed"))
2563	{
2564	treeAgeFixed = YES;
2565	treeAge = mp->treeAgeFix;
2566	}
2567	else
2568	{
2569	treeAgeFixed = NO;
2570	treeAge = -1.0;
2571	}
2572
2573	/* date all nodes from top to bottom with min. age as nodeDepth*/
2574	for (i=0; i<t->nNodes; i++)
2575	{
2576	p = t->allDownPass[i];
2577	if (p->anc != NULL)
2578	{
2579	if (p->left == NULL && p->right == NULL)
2580	{
2581	if (p->isDated == NO)
2582	{
2583	p->nodeDepth = 0.0;
2584	p->age = 0.0;
2585	}
2586	else
2587	{
2588	if (p->calibration->prior == fixed)
2589	p->nodeDepth = p->age = p->calibration->age;
2590	else if (p->calibration->prior == uniform)
2591	p->nodeDepth = p->age = p->calibration->min;
2592	else /* if (p->calibration->prior == offsetExponential) */
2593	p->nodeDepth = p->age = p->calibration->offset;
2594	}
2595	}
2596	else
2597	{
2598	if (p->left->nodeDepth > p->right->nodeDepth)
2599	p->nodeDepth = p->left->nodeDepth;
2600	else
2601	p->nodeDepth = p->right->nodeDepth;
2602	if (p->isDated == YES \|\| (p->anc->anc == NULL && treeAgeFixed))
2603	{
2604	if (p->isDated == NO)
2605	{
2606	if (treeAge <= p->nodeDepth)
2607	{
2608	MrBayesPrint ("%s Calibration inconsistency for root node\n", spacer);
2609	return (ERROR);
2610	}
2611	else
2612	p->age = p->nodeDepth = treeAge;
2613	}
2614	else if (p->calibration->prior == fixed)
2615	{
2616	if (p->calibration->age <= p->nodeDepth)
2617	{
2618	if (p->anc->anc == NULL)
2619	MrBayesPrint ("%s Calibration inconsistency for root node\n", spacer);
2620	else
2621	MrBayesPrint ("%s Calibration inconsistency for node '%s'\n", spacer, constraintNames[p->lockID]);
2622	return (ERROR);
2623	}
2624	else
2625	p->age = p->nodeDepth = p->calibration->age;
2626	}
2627	else if (p->calibration->prior == uniform)
2628	{
2629	if (p->calibration->max <= p->nodeDepth)
2630	{
2631	MrBayesPrint ("%s Calibration inconsistency for node '%s'\n", spacer, constraintNames[p->lockID]);
2632	return (ERROR);
2633	}
2634	else
2635	{
2636	if (p->calibration->min < p->nodeDepth)
2637	p->age = p->nodeDepth;
2638	else
2639	p->age = p->nodeDepth = p->calibration->min;
2640	}
2641	}
2642	else /* if (p->calibration.prior == offsetExponential) */
2643	{
2644	if (p->calibration->offset < p->nodeDepth)
2645	p->age = p->nodeDepth;
2646	else
2647	p->age = p->nodeDepth = p->calibration->offset;
2648	}
2649	}
2650	else
2651	p->age = -1.0;
2652	}
2653	}
2654	}
2655
2656
2657	/* try to make root node deeper than minimum age */
2658	p = t->root->left;
2659	if (p->calibration == NULL && treeAgeFixed == NO)
2660	{
2661	if( p->nodeDepth==0.0 ) p->nodeDepth = 1.0;
2662	p->age = p->nodeDepth *= 1.5;
2663	}
2664	else if (treeAgeFixed == YES \|\| p->calibration->prior == fixed)
2665	{
2666	/* can't do much ... */
2667	}
2668	else if (p->calibration->prior == uniform)
2669	{
2670	if (p->nodeDepth * 1.5 < p->calibration->max)
2671	p->nodeDepth = p->age = 1.5 * p->nodeDepth;
2672	else
2673	p->nodeDepth = p->age = p->calibration->max;
2674	}
2675	else /* if (t->root->calibration->prior == offsetExponential */
2676	{
2677	assert( p->calibration->prior == offsetExponential );
2678	/* Make random draw */
2679	p->age = p->calibration->offset - log (RandomNumber(seed)) / p->calibration->lambda;
2680	if (p->age > 1.5 * p->nodeDepth)
2681	p->nodeDepth = p->age;
2682	else
2683	p->age = p->nodeDepth = 1.5 * p->nodeDepth;
2684	}
2685
2686	SetNodeCalibratedAge( p, 1, p->age );
2687
2688
2689	/* Setup node depths */
2690	for (i=0; i<t->nNodes; i++)
2691	{
2692	p = t->allDownPass[i];
2693	p->nodeDepth = p->age * clockRate;
2694	assert (!(p->left == NULL && p->calibration == NULL && p->nodeDepth != 0.0));
2695	}
2696
2697	/* calculate branch lengths */
2698	for (i=0; i<t->nNodes; i++)
2699	{
2700	p = t->allDownPass[i];
2701	if (p->anc != NULL)
2702	{
2703	if (p->anc->anc != NULL)
2704	{
2705	p->length = p->anc->nodeDepth - p->nodeDepth;
2706	if( p->length < BRLENS_MIN )
2707	{
2708	//MrBayesPrint ("%s Restrictions of node calibration and clockrate makes some branch lenghts too small.\n", spacer);
2709	//return (ERROR);
2710	}
2711	if( p->length > BRLENS_MAX )
2712	{
2713	//MrBayesPrint ("%s Restrictions of node calibration and clockrate makes some branch lenghts too long.\n", spacer);
2714	//return (ERROR);
2715	}
2716	}
2717	else
2718	p->length = 0.0; //not a problem for root node.
2719	}
2720	}
2721
2722	#if 0
2723	printf ("after\n");
2724	ShowNodes (t->root, 0, YES);
2725	getchar();
2726	#endif
2727
2728	return (NO_ERROR);
2729
2730	}
2731
2732
2733
2734
2735
2736	/*-------------------------------------------------------
2737	\|
2738	\| InitClockBrlens: This routine will initialize
2739	\| a clock tree by setting the root to depth 1.0
2740	\| and then assigning node depths according to
2741	\| the relative node depth measured in terms of the
2742	\| maximum number of branches to the tip from each
2743	\| node.
2744	\|
2745	--------------------------------------------------------*/
2746	int InitClockBrlens (Tree *t)
2747
2748	{
2749
2750	int i, maxBrSegments=0;
2751	TreeNode *p;
2752
2753	if (t->isRooted == NO)
2754	{
2755	MrBayesPrint ("%s Tree is unrooted\n", spacer);
2756	return (ERROR);
2757	}
2758
2759	/* calculate maximum number of branch segments above root */
2760	for (i=0; i<t->nNodes; i++)
2761	{
2762	p = t->allDownPass[i];
2763	if (p->anc != NULL)
2764	{
2765	if (p->left == NULL && p->right == NULL)
2766	{
2767	p->x = 0;
2768	}
2769	else
2770	{
2771	if (p->left->x > p->right->x)
2772	p->x = p->left->x + 1;
2773	else
2774	p->x = p->right->x + 1;
2775	}
2776	if (p->anc->anc == NULL)
2777	maxBrSegments = p->x;
2778	}
2779	}
2780
2781	/* assign node depths */
2782	for (i=0; i<t->nNodes; i++)
2783	{
2784	p = t->allDownPass[i];
2785	if (p->anc != NULL)
2786	p->nodeDepth = (MrBFlt) (p->x) / (MrBFlt) maxBrSegments;
2787	else
2788	p->nodeDepth = 0.0;
2789	}
2790
2791	/* calculate branch lengths */
2792	for (i=0; i<t->nNodes; i++)
2793	{
2794	p = t->allDownPass[i];
2795	if (p->anc != NULL)
2796	{
2797	if (p->anc->anc != NULL)
2798	p->length = p->anc->nodeDepth - p->nodeDepth;
2799	else
2800	p->length = 0.0;
2801	}
2802	}
2803
2804	return (NO_ERROR);
2805
2806	}
2807
2808
2809
2810
2811
2812	int GetRandomEmbeddedSubtree (Tree t, int nTerminals, SafeLong seed, int *nEmbeddedTrees)
2813
2814	{
2815
2816	int i, j, k, n, ran, pP, pL, pR, nLeaves, nSubTrees;
2817	TreeNode p=NULL, *leaf;
2818
2819	/* Calculate number of leaves in subtree (number of terminals minus the root) */
2820	nLeaves = nTerminals - 1;
2821
2822	/* Initialize all flags */
2823	for (i=0; i<t->nNodes; i++)
2824	{
2825	p = t->allDownPass[i];
2826	p->marked = NO;
2827	p->x = 0;
2828	p->y = 0;
2829	}
2830
2831	/* Allocate memory */
2832	nSubTrees = (int ) SafeCalloc (nTerminals t->nNodes, sizeof(int));
2833	if (!nSubTrees)
2834	return (ERROR);
2835	leaf = (TreeNode *) SafeMalloc (nLeaves sizeof (TreeNode *));
2836	if (!leaf)
2837	{
2838	free (nSubTrees);
2839	return (ERROR);
2840	}
2841
2842	/* Calculate how many embedded trees rooted at each node */
2843	(*nEmbeddedTrees) = 0;
2844	for (i=0; i<t->nNodes-1; i++)
2845	{
2846	p = t->allDownPass[i];
2847	if (p->left == NULL)
2848	{
2849	p->x = 0;
2850	nSubTrees[p->index*nTerminals + 1] = 1;
2851	}
2852	else
2853	{
2854	pL = nSubTrees + p->left->index*nTerminals;
2855	pR = nSubTrees + p->right->index*nTerminals;
2856	pP = nSubTrees + p->index*nTerminals;
2857	pP[1] = 1;
2858	for (j=2; j<=nLeaves; j++)
2859	{
2860	for (k=1; k<j; k++)
2861	{
2862	pP[j] += pL[k] * pR[j-k];
2863	}
2864	}
2865	p->x = pP[nLeaves];
2866	(*nEmbeddedTrees) += p->x;
2867	}
2868	}
2869
2870	/* Randomly select one embedded tree of the right size */
2871	ran = (int) (RandomNumber(seed) * (*nEmbeddedTrees));
2872
2873	/* Find the interior root corresponding to this tree */
2874	for (i=j=0; i<t->nIntNodes; i++)
2875	{
2876	p = t->intDownPass[i];
2877	j += p->x;
2878	if (j>ran)
2879	break;
2880	}
2881
2882	/* Find one random embedded tree with this root */
2883	p->y = nLeaves;
2884	p->marked = YES;
2885	leaf[0] = p;
2886	n = 1;
2887	while (n < nLeaves)
2888	{
2889	/* select a node with more than one descendant */
2890	for (i=0; i<n; i++)
2891	{
2892	p = leaf[i];
2893	if (p->y > 1)
2894	break;
2895	}
2896
2897	/* break it into descendants */
2898	pL = nSubTrees + p->left->index*nTerminals;
2899	pR = nSubTrees + p->right->index*nTerminals;
2900	pP = nSubTrees + p->index*nTerminals;
2901	ran = (int) (RandomNumber (seed) * pP[p->y]);
2902	k = 0;
2903	for (j=1; j<p->y; j++)
2904	{
2905	k += pL[j] * pR[p->y-j];
2906	if (k > ran)
2907	break;
2908	}
2909
2910	p->left->y = j;
2911	p->right->y = p->y - j;
2912	p->left->marked = YES;
2913	p->right->marked = YES;
2914	leaf[i] = p->left;
2915	leaf[n++] = p->right;
2916	}
2917
2918	free (nSubTrees);
2919	free (leaf);
2920
2921	return (NO_ERROR);
2922	}
2923
2924
2925
2926
2927
2928	/*-----------------------------------------------------------------------------
2929	\|
2930	\| IsCalibratedClockSatisfied: This routine SETS(not just checks as name suggest) calibrated clock tree nodes age, depth. based on branch lengthes
2931	\| and checks that user defined brlens satisfy the specified calibration(s)
2932	\| up to tolerance tol
2933	\| TODO clock rate is devived here and used to set ages but clockrate paramiter is not updated here(make sure that it does not produce imconsitancy)
2934	\|
2935	\|------------------------------------------------------------------------------*/
2936	int IsCalibratedClockSatisfied (Tree t,MrBFlt minClockRate,MrBFlt *maxClockRate , MrBFlt tol)
2937
2938	{
2939
2940	int i, j, maxRateConstrained, minRateConstrained, isViolated;
2941	MrBFlt f, maxHeight, minRate=0, maxRate=0, ageToAdd, x, y, clockRate;
2942	TreeNode p, q, r, s;
2943
2944	/By defauult assume the tree does not have allowed range of clockrate/
2945	*minClockRate = 2.0;
2946	*maxClockRate = 1.0;
2947
2948	if (t->isRooted == NO)
2949	return (NO);
2950
2951	x = (MrBFlt ) SafeCalloc (2t->nNodes, sizeof (MrBFlt));
2952	if (x == NULL)
2953	{
2954	MrBayesPrint ("%s Out of memory in IsCalibratedClockSatisfied\n", spacer);
2955	free (x);
2956	return (NO);
2957	}
2958	y = x + t->nNodes;
2959
2960	/* reset node depth and age, and set minimum (x) and maximum (y) age of each node */
2961	for (i=0; i<t->nNodes; i++)
2962	{
2963	p = t->allDownPass[i];
2964	p->age = -1.0;
2965	p->nodeDepth = -1.0;
2966	if (p->isDated == YES)
2967	{
2968	assert(p->calibration->prior == fixed \|\| p->calibration->prior == uniform \|\| p->calibration->prior == offsetExponential);
2969	if (p->calibration->prior == fixed)
2970	x[p->index] = y[p->index] = p->calibration->age;
2971	else if (p->calibration->prior == uniform)
2972	{
2973	x[p->index] = p->calibration->min;
2974	y[p->index] = p->calibration->max;
2975	}
2976	else /* if (p->calibration->prior == offsetExponential) */
2977	{
2978	x[p->index] = p->calibration->offset;
2979	y[p->index] = -1.0;
2980	}
2981	}
2982	else if (p->left == NULL && p->right == NULL)
2983	x[p->index] = y[p->index] = 0.0;
2984	else
2985	{
2986	x[p->index] = y[p->index] = -1.0;
2987	}
2988	}
2989
2990	/* calculate node heights in branch length units */
2991	/* node depth will be set from the root for now*/
2992	p = t->root->left;
2993	p->nodeDepth = 0.0;
2994	for (i=t->nNodes-3; i>=0; i--)
2995	{
2996	p = t->allDownPass[i];
2997	p->nodeDepth = p->anc->nodeDepth + p->length;
2998	}
2999
3000	/* find maximum height of tree */
3001	maxHeight = -1.0;
3002	for (i=0; i<t->nNodes-1; i++)
3003	{
3004	p = t->allDownPass[i];
3005	if (p->left == NULL && p->right == NULL)
3006	{
3007	if (p->nodeDepth > maxHeight)
3008	{
3009	maxHeight = p->nodeDepth;
3010	}
3011	}
3012	}
3013
3014	/* calculate node depth from tip of tree */
3015	for (i=0; i<t->nNodes-1; i++)
3016	{
3017	p = t->allDownPass[i];
3018	p->nodeDepth = maxHeight - p->nodeDepth;
3019	}
3020
3021
3022	/* check potentially constraining calibrations */
3023	/* and find minimum and maximum possible rate */
3024	maxRateConstrained = NO;
3025	minRateConstrained = NO;
3026	isViolated = NO;
3027	for (i=0; i<t->nNodes-1; i++)
3028	{
3029	p = t->allDownPass[i];
3030	if (x[p->index] < 0.0 && y[p->index] < 0.0)
3031	continue;
3032	for (j=i+1; j<t->nNodes-1; j++)
3033	{
3034	q = t->allDownPass[j];
3035	if (x[q->index] < 0.0 && y[q->index] < 0.0)
3036	continue;
3037	if ( p->nodeDepth == q->nodeDepth) // becouse clock rate could be as low as possible we can not take approximate equality.
3038	{
3039	/* same depth so they must share a possible age */
3040	if ((x[p->index] != -1.0 && y[q->index] !=-1.0 && AreDoublesEqual (x[p->index], y[q->index], tol) == NO && x[p->index] > y[q->index])
3041	\|\| (y[p->index] != -1.0 && x[q->index]!=-1.0 && AreDoublesEqual (y[p->index], x[q->index], tol) == NO && y[p->index] < x[q->index]))
3042	{
3043	isViolated = YES;
3044	break;
3045	}
3046	}
3047	else
3048	{
3049	if (p->nodeDepth > q->nodeDepth)
3050	{
3051	r = p;
3052	s = q;
3053	}
3054	else
3055	{
3056	r = q;
3057	s = p;
3058	}
3059	if (x[r->index] >= 0.0 && y[s->index] >= 0.0)
3060	{
3061	f = (r->nodeDepth - s->nodeDepth) / (x[r->index] - y[s->index]);
3062	if (f <= 0.0 \|\| x[r->index] == y[s->index])
3063	{
3064	if ( AreDoublesEqual (r->nodeDepth, s->nodeDepth, 0.000001) == YES) //if defference is very very small we do not bail out. It could heppened becouse of numerical inacuracy, one node which supose to be slitly below the other one become on top.
3065	continue;
3066	isViolated = YES;
3067	break;
3068	}
3069	if (maxRateConstrained == NO)
3070	{
3071	maxRateConstrained = YES;
3072	maxRate = f;
3073	}
3074	else if (f < maxRate)
3075	maxRate = f;
3076	}
3077	if (y[r->index] >= 0.0 && x[s->index] >= 0.0)
3078	{
3079	f = (r->nodeDepth - s->nodeDepth) / (y[r->index] - x[s->index]);
3080	if (f <= 0.0 \|\| y[r->index] == x[s->index])
3081	{
3082	if ( AreDoublesEqual (r->nodeDepth, s->nodeDepth, 0.00001) == YES)
3083	continue;
3084	isViolated = YES;
3085	break;
3086	}
3087	if (minRateConstrained == NO)
3088	{
3089	minRateConstrained = YES;
3090	minRate = f;
3091	}
3092	else if (f > minRate)
3093	minRate = f;
3094	}
3095	}
3096	}
3097	if (isViolated == YES)
3098	break;
3099	}
3100
3101	/* check if outright violation */
3102	if (isViolated == YES)
3103	{
3104	MrBayesPrint ("%s Branch lengths do not satisfy the calibration(s)\n", spacer);
3105	free (x);
3106	return (NO);
3107	}
3108
3109	/Allow tollerance/
3110	if(minRateConstrained == YES && maxRateConstrained == YES && AreDoublesEqual (minRate, maxRate, tol) == YES && minRate > maxRate)
3111	{
3112	maxRate = minRate;
3113	}
3114
3115	if (minRateConstrained == YES)
3116	*minClockRate = minRate;
3117	else
3118	*minClockRate = 0.0;
3119
3120	if (maxRateConstrained == YES)
3121	*maxClockRate = maxRate;
3122	else
3123	*maxClockRate = MRBFLT_MAX;
3124
3125	/* check that minimum and maximum rates are consistent */
3126	if (minRateConstrained == YES && maxRateConstrained == YES && minRate > maxRate)
3127	{
3128	MrBayesPrint ("%s Branch lengths do not satisfy the calibration(s)\n", spacer);
3129	free (x);
3130	return (NO);
3131	}
3132
3133	/* date all nodes based on a suitable rate */
3134	if (minRateConstrained == YES)
3135	clockRate = minRate;
3136	else if (maxRateConstrained == YES)
3137	clockRate = 0.5 * maxRate;
3138	else
3139	clockRate = 1.0;
3140	for (i=0; i<t->nNodes-1; i++)
3141	{
3142	p = t->allDownPass[i];
3143	p->age = p->nodeDepth / clockRate;
3144	}
3145
3146	/* check if there is an age to add (I guess this is here because when max rate is close to minrate and we have numerical precision inacuracy)*/
3147	ageToAdd = 0.0;
3148	for (i=0; i<t->nNodes-1; i++)
3149	{
3150	p = t->allDownPass[i];
3151	if (x[p->index] > 0.0 && x[p->index] > p->age)
3152	{
3153	f = x[p->index] - p->age;
3154	if (f > ageToAdd)
3155	ageToAdd = f;
3156	}
3157	}
3158
3159	/* add extra length if any */
3160	if (AreDoublesEqual (ageToAdd, 0.0, 0.00000001) == NO)
3161	{
3162	for (i=0; i<t->nNodes-1; i++)
3163	{
3164	p = t->allDownPass[i];
3165	p->age += ageToAdd;
3166	}
3167	}
3168
3169	free (x);
3170
3171	/* reset node depths to ensure that non-dated tips have node depth 0.0 */
3172	SetNodeDepths(t);
3173
3174	return (YES);
3175
3176	}
3177
3178
3179
3180
3181
3182	int IsClockSatisfied (Tree *t, MrBFlt tol)
3183
3184	{
3185
3186	int i, foundFirstLength, isClockLike;
3187	MrBFlt firstLength=0.0, length;
3188	TreeNode p, q;
3189
3190	if (t->isRooted == NO)
3191	return (NO);
3192
3193	foundFirstLength = NO;
3194	isClockLike = YES;
3195	for (i=0; i<t->nNodes; i++)
3196	{
3197	p = t->allDownPass[i];
3198	if (p->left == NULL && p->right == NULL)
3199	{
3200	if (p->isDated == YES)
3201	{
3202	//continue;
3203	length = p->nodeDepth;
3204	}
3205	else
3206	length = 0.0;
3207	q = p;
3208	while (q->anc != NULL)
3209	{
3210	if (q->anc->anc != NULL)
3211	length += q->length;
3212	q = q->anc;
3213	}
3214	if (foundFirstLength == NO)
3215	{
3216	firstLength = length;
3217	foundFirstLength = YES;
3218	}
3219	else
3220	{
3221	if (AreDoublesEqual (firstLength, length, tol) == NO)
3222	{
3223	MrBayesPrint ("%s Node (%s) is not at the same depth as some other tip taking colibration into account. \n", spacer, p->label);
3224	isClockLike = NO;
3225	}
3226	}
3227	}
3228	}
3229	if (firstLength < BRLENS_MIN)
3230	isClockLike = NO;
3231
3232	return (isClockLike);
3233
3234	}
3235
3236
3237
3238
3239
3240	/* Check that tree obeys topology constraints and that node depths and ages are consistent */
3241	int IsTreeConsistent (Param *param, int chain, int state)
3242	{
3243	Tree *tree;
3244	TreeNode *p;
3245	int i, j;
3246	MrBFlt b, r, rAnc, clockRate;
3247
3248	if (param->paramType != P_TOPOLOGY && param->paramType != P_BRLENS && param->paramType != P_SPECIESTREE)
3249	return YES;
3250
3251	tree = GetTree(param, chain, state);
3252	if (modelSettings[param->relParts[0]].clockRate != NULL)
3253	clockRate = *GetParamVals(modelSettings[param->relParts[0]].clockRate, chain, state);
3254	else
3255	clockRate = 0.0;
3256
3257	if (CheckConstraints(tree)==ERROR) {
3258	printf ("Tree does not obey constraints\n");
3259	return NO;
3260	}
3261
3262	/* check that the last few indices are not taken in a rooted tree */
3263	if (tree->isRooted == YES && tree->root->index != tree->nNodes - 1)
3264	{
3265	printf("Problem with root index\n");
3266	return NO;
3267	}
3268	if (tree->isRooted == YES && tree->root->left->index != tree->nNodes - 2)
3269	{
3270	printf("Problem with interior root index\n");
3271	return NO;
3272	}
3273
3274	if (tree->isClock == NO)
3275	{
3276	for (i=0; i<tree->nNodes-1; i++)
3277	{
3278	p = tree->allDownPass[i];
3279	if (p->length <= 0.0)
3280	{
3281	if(p->length == 0.0)
3282	printf ("Node %d has zero branch length %f\n", p->index, p->length);
3283	else
3284	printf ("Node %d has negative branch length %f\n", p->index, p->length);
3285	return NO;
3286	}
3287	}
3288	return YES;
3289	}
3290
3291	/* Clock trees */
3292
3293	/* Check that lengths and depths are consistent */
3294	for (i=0; i<tree->nNodes-2; i++) {
3295	p = tree->allDownPass[i];
3296	if (p->length <= 0.0) {
3297	if(p->length == 0.0)
3298	printf ("Node %d has zero branch length %f\n", p->index, p->length);
3299	else
3300	printf ("Node %d has negative branch length %f\n", p->index, p->length);
3301	return NO;
3302	}
3303	if (fabs(p->anc->nodeDepth - p->nodeDepth - p->length) > 0.000001) {
3304	printf ("Node %d has length %f but nodeDepth %f and ancNodeDepth %f\n",
3305	p->index, p->length, p->nodeDepth, p->anc->nodeDepth);
3306	return NO;
3307	}
3308	if (p->left == NULL && p->isDated == NO && p->nodeDepth != 0.0) {
3309	printf ("Node %d is a nondated tip but has node depth %f\n",
3310	p->index, p->nodeDepth);
3311	//return NO;
3312	}
3313	}
3314
3315	/* Check that ages and calibrations are consistent */
3316	if (tree->isCalibrated == YES)
3317	{
3318	for (i=0; i<tree->nNodes-1; i++)
3319	{
3320	p = tree->allDownPass[i];
3321	if (p->isDated == YES) {
3322	if (fabs(p->age - p->nodeDepth/clockRate) > 0.000001)
3323	{
3324	printf ("Node %d has age %f but nodeDepth %f when clock rate is %f\n",
3325	p->index, p->age, p->nodeDepth, clockRate);
3326	return NO;
3327	}
3328	if (p->calibration->prior == fixed && fabs(p->age - p->calibration->age) > 0.000001)
3329	{
3330	printf ("Node %d has age %f but should be fixed to age %f\n",
3331	p->index, p->age, p->calibration->age);
3332	return NO;
3333	}
3334	else if (p->calibration->prior == offsetExponential && p->age < p->calibration->offset)
3335	{
3336	printf ("Node %d has age %f but should be minimally of age %f\n",
3337	p->index, p->age, p->calibration->offset);
3338	return NO;
3339	}
3340	else if (p->calibration->prior == uniform && (p->age < p->calibration->min \|\| p->age > p->calibration->max))
3341	{
3342	printf ("Node %d has age %f but should be in the interval (%f,%f)\n",
3343	p->index, p->age, p->calibration->min, p->calibration->max);
3344	return NO;
3345	}
3346	}
3347	}
3348	}
3349
3350
3351	for (i=0; i<param->nSubParams; i++)
3352	{
3353	if (param->subParams[i]->paramId == TK02BRANCHRATES)
3354	{
3355	rAnc = GetParamVals(param->subParams[i], chain, state)[tree->root->left->index];
3356	if (fabs(rAnc - 1.0) > 1E-6)
3357	{
3358	MrBayesPrint("%s TK02 relaxed clock mismatch in root rate, which is %e\n", spacer, rAnc);
3359	return NO;
3360	}
3361	for (j=0; j<tree->nNodes-2; j++)
3362	{
3363	p = tree->allDownPass[j];
3364	b = GetParamSubVals(param->subParams[i], chain, state)[p->index];
3365	r = GetParamVals(param->subParams[i], chain, state)[p->index];
3366	rAnc = GetParamVals(param->subParams[i], chain, state)[p->anc->index];
3367	if (fabs(p->length * (r + rAnc) / 2.0 - b) > 0.000001)
3368	{
3369	MrBayesPrint("%s TK02 relaxed clock mismatch in branch %d\n", spacer, p->index);
3370	return NO;
3371	}
3372	}
3373	}
3374	else if (param->subParams[i]->paramId == IGRBRANCHLENS)
3375	{
3376	for (j=0; j<tree->nNodes-2; j++)
3377	{
3378	p = tree->allDownPass[j];
3379	b = GetParamSubVals(param->subParams[i], chain, state)[p->index];
3380	r = GetParamVals(param->subParams[i], chain, state)[p->index];
3381	if (fabs(p->length * r - b) > 0.000001)
3382	{
3383	MrBayesPrint("%s Igr relaxed clock mismatch in branch %d\n", spacer, p->index);
3384	return NO;
3385	}
3386	}
3387	}
3388	}
3389
3390	if (param->paramType == P_SPECIESTREE)
3391	return (IsSpeciesTreeConsistent(GetTree(param, chain, state), chain));
3392
3393	return YES;
3394	}
3395
3396
3397
3398
3399
3400	/* LabelTree: Label tree; remove previous labels if any */
3401	int LabelTree (Tree t, char *taxonNames)
3402	{
3403	int i, nTaxa;
3404	TreeNode *p = NULL;
3405
3406	nTaxa = t->nNodes - t->nIntNodes;
3407	if (t->isRooted == YES)
3408	nTaxa--;
3409
3410	/* erase previous labels, if any */
3411	for (i=0; i<t->nNodes; i++)
3412	{
3413	p = t->allDownPass[i];
3414	p->marked = NO;
3415	t->nodes[i].label = noLabel;
3416	}
3417
3418	/* add labels */
3419	for (i=0; i<t->nNodes; i++)
3420	{
3421	p = &t->nodes[i];
3422	if (p->left == NULL \|\| (t->isRooted == NO && p->anc == NULL))
3423	{
3424	if (p->marked == YES \|\| p->index < 0 \|\| p->index >= nTaxa)
3425	{
3426	MrBayesPrint ("%s Taxon node index repeated or out of range\n", spacer);
3427	return (ERROR);
3428	}
3429	else
3430	p->label = taxonNames[p->index];
3431	p->marked = YES;
3432	}
3433	else if (p->index > 0 && p->index < nTaxa)
3434	{
3435	MrBayesPrint ("%s Terminal taxon index set for interior node\n", spacer);
3436	return (ERROR);
3437	}
3438	}
3439
3440	return (NO_ERROR);
3441	}
3442
3443
3444
3445
3446
3447	/*-------------------------------------------------------------------------------------------
3448	\|
3449	\| Mark: This routine will mark up a subtree rooted at p
3450	\|
3451	---------------------------------------------------------------------------------------------*/
3452	void Mark (TreeNode *p)
3453	{
3454	if (p != NULL)
3455	{
3456	p->marked = YES;
3457	Mark (p->left);
3458	Mark (p->right);
3459	}
3460	}
3461
3462
3463
3464
3465
3466	/*-------------------------------------------------------------------------------------------
3467	\|
3468	\| MarkDatedSubtree: This routine will mark up a subtree rooted at p with dated tips, whether
3469	\| internal or external
3470	\|
3471	---------------------------------------------------------------------------------------------*/
3472	void MarkDatedSubtree (TreeNode *p)
3473	{
3474	if (p != NULL)
3475	{
3476	p->marked = YES;
3477	MarkDatedSubtreeNodes (p->left);
3478	MarkDatedSubtreeNodes (p->right);
3479	}
3480	}
3481
3482
3483
3484
3485
3486	/* MarkDatedSubtreeNodes: Recursive function to mark parts of a dated subtree */
3487	void MarkDatedSubtreeNodes (TreeNode *p)
3488	{
3489	if (p != NULL)
3490	{
3491	p->marked = YES;
3492	if (p->isDated == NO && p->left != NULL)
3493	{
3494	MarkDatedSubtreeNodes (p->left);
3495	MarkDatedSubtreeNodes (p->right);
3496	}
3497	}
3498	}
3499
3500
3501
3502
3503
3504
3505	/*-------------------------------------------------------------------------------------------
3506	\|
3507	\| MoveCalculationRoot: This routine will move the calculation root to the terminal with
3508	\| index outgroup
3509	\|
3510	---------------------------------------------------------------------------------------------*/
3511	int MoveCalculationRoot (Tree *t, int outgroup)
3512
3513	{
3514
3515	int i;
3516	TreeNode p, q, *r;
3517
3518	if (t->isRooted == YES \|\| outgroup < 0 \|\| outgroup > t->nNodes - t->nIntNodes - (t->isRooted == YES ? 1 : 0))
3519	{
3520	MrBayesPrint ("%s Problem moving calculation root\n", spacer);
3521	return (ERROR);
3522	}
3523
3524	if (t->root->index == outgroup)
3525	return (NO_ERROR); /* nothing to do */
3526
3527	/* mark the path to the new calculation root */
3528	for (i=0; i<t->nNodes; i++)
3529	{
3530	p = t->allDownPass[i];
3531	if (p->left == NULL && p->right == NULL)
3532	{
3533	if (p->index == outgroup)
3534	p->marked = YES;
3535	else
3536	p->marked = NO;
3537	}
3538	else
3539	{
3540	if (p->left->marked == YES \|\| p->right->marked == YES)
3541	p->marked = YES;
3542	else
3543	p->marked = NO;
3544	}
3545	}
3546
3547	/* rotate the tree to use the specified calculation root */
3548	p = t->root->left;
3549	q = t->root;
3550	q->anc = p;
3551	q->left = q->right = NULL;
3552	q->length = p->length;
3553	while (p->left != NULL && p->right != NULL)
3554	{
3555	if (p->left->marked == YES)
3556	{
3557	r = p->left;
3558	p->anc = r;
3559	p->left = q;
3560	p->length = r->length;
3561	q = p;
3562	p = r;
3563	}
3564	else /* if (p->right->marked == YES) */
3565	{
3566	r = p->right;
3567	p->anc = r;
3568	p->right = q;
3569	p->length = r->length;
3570	q = p;
3571	p = r;
3572	}
3573	}
3574	p->left = p->anc;
3575	p->right = p->anc = NULL;
3576	t->root = p;
3577	p->length = 0.0;
3578
3579	GetDownPass (t);
3580
3581	return (NO_ERROR);
3582	}
3583
3584
3585
3586
3587
3588	/*-------------------------------------------------------------------------------------------
3589	\|
3590	\| MovePolyCalculationRoot: This routine will move the calculation root to the terminal with
3591	\| index outgroup and place it as the right-most descendant of the root node
3592	\|
3593	---------------------------------------------------------------------------------------------*/
3594	int MovePolyCalculationRoot (PolyTree *t, int outgroup)
3595
3596	{
3597
3598	int i;
3599	PolyNode p = NULL, q, *r;
3600
3601	/* check if tree is rooted, in which case calculation root is irrelevant */
3602	if (t->root->left->sib->sib == NULL)
3603	return (NO_ERROR);
3604
3605	if (outgroup < 0 \|\| outgroup > t->nNodes - t->nIntNodes)
3606	{
3607	MrBayesPrint ("%s Outgroup index is out of range\n", spacer);
3608	return (ERROR);
3609	}
3610
3611	if (t->root->left->sib->sib->sib != NULL)
3612	{
3613	MrBayesPrint ("%s Root has more than three descendants\n", spacer);
3614	return (ERROR);
3615	}
3616
3617	/* check if rerooting actually necessary */
3618	if (t->root->left->sib->sib->index == outgroup)
3619	return (NO_ERROR);
3620
3621	/* mark the path to the new calculation root */
3622	for (i=0; i<t->nNodes; i++)
3623	{
3624	p = t->allDownPass[i];
3625	if (p->index == outgroup)
3626	break;
3627	}
3628	if (p->left != NULL)
3629	{
3630	MrBayesPrint ("%s Outgroup index set for internal node\n", spacer);
3631	for (i=0; i<t->nNodes; i++)
3632	printf ("%d -- %d\n", i, t->allDownPass[i]->index);
3633	getchar();
3634	return (ERROR);
3635	}
3636
3637	/* mark path to current root */
3638	for (i=0; i<t->nNodes; i++)
3639	t->allDownPass[i]->mark = NO;
3640	q = p;
3641	while (q != NULL)
3642	{
3643	q->mark = YES;
3644	q = q->anc;
3645	}
3646
3647	/* rotate the tree to use the specified calculation root */
3648	p = t->root;
3649	for (;;)
3650	{
3651	/* find marked descendant */
3652	for (q=p->left; q->mark == NO; q=q->sib)
3653	;
3654	if (q->index == outgroup)
3655	break;
3656	/* add old root to descendants of that node */
3657	for (r=q->left; r->sib!=NULL; r=r->sib)
3658	;
3659	r->sib = p;
3660	p->sib = NULL; /* should not be needed */
3661	p->anc = q;
3662	p->length = q->length;
3663	/* remove that node from descendants of old root node */
3664	if (p->left == q)
3665	p->left = q->sib;
3666	else
3667	{
3668	for (r=p->left; r->sib!=q; r=r->sib)
3669	;
3670	r->sib = r->sib->sib;
3671	}
3672	/* make new node root */
3673	q->sib = NULL;
3674	q->anc = NULL;
3675	q->length = 0.0;
3676	p = q;
3677	}
3678
3679	/* p is now the new root */
3680	t->root = p;
3681
3682	/* finally make sure calculation root is last node among root's descendants */
3683	for (q=p->left; q->sib!=NULL; q=q->sib)
3684	;
3685	if (q->index != outgroup)
3686	{
3687	if (p->left->index == outgroup)
3688	{
3689	q->sib = p->left;
3690	p->left = p->left->sib;
3691	q->sib->sib = NULL;
3692	}
3693	else
3694	{
3695	for (r=p->left; r->sib->index!=outgroup; r=r->sib)
3696	;
3697	q->sib = r->sib;
3698	r->sib = r->sib->sib;
3699	q->sib->sib = NULL;
3700	}
3701	}
3702
3703	GetPolyDownPass (t);
3704
3705	return (NO_ERROR);
3706	}
3707
3708
3709
3710
3711
3712	/*
3713	@return the number of levels for the tree rooted at the "node"
3714	*/
3715	int NrSubTreeLevels(TreeNode *node)
3716	{
3717	int r,l;
3718
3719	if( node == NULL )
3720	{
3721	return -1;
3722	}
3723
3724	r = NrSubTreeLevels( node->right );
3725	l = NrSubTreeLevels( node->left );
3726
3727	return ((r>l)?(r):(l))+1;
3728	}
3729
3730
3731
3732
3733
3734	/*-------------------------------------------------------------------------------------------
3735	\|
3736	\| NumConstrainedTips: This routine will return the number of constrained tips, internal or external
3737	\|
3738	---------------------------------------------------------------------------------------------*/
3739	int NumConstrainedTips (TreeNode *p)
3740	{
3741	int i = 0;
3742
3743	if (p == NULL)
3744	return i;
3745	if (p->left == NULL)
3746	return 1;
3747
3748	i += NConstrainedTips (p->left);
3749	i += NConstrainedTips (p->right);
3750
3751	return i;
3752	}
3753
3754
3755
3756
3757
3758	/* NConstrainedTips: Recursive function to get the number of constrained tips */
3759	int NConstrainedTips (TreeNode *p)
3760	{
3761	int i=0;
3762
3763	if (p!=NULL)
3764	{
3765	if (p->left == NULL \|\| p->isLocked == YES)
3766	return 1;
3767	else
3768	{
3769	i += NConstrainedTips (p->left);
3770	i += NConstrainedTips (p->right);
3771	}
3772	}
3773	return i;
3774	}
3775
3776
3777
3778
3779
3780	/*-------------------------------------------------------------------------------------------
3781	\|
3782	\| NumDatedTips: This routine will return the number of dated tips, internal or external
3783	\|
3784	---------------------------------------------------------------------------------------------*/
3785	int NumDatedTips (TreeNode *p)
3786	{
3787	int i = 0;
3788
3789	assert (p != NULL && p->left != NULL);
3790
3791	i += NDatedTips (p->left);
3792	i += NDatedTips (p->right);
3793
3794	return i;
3795	}
3796
3797
3798
3799
3800
3801	/* NDatedTips: recursive function to get the number of dated tips */
3802	int NDatedTips (TreeNode *p)
3803	{
3804	int i=0;
3805
3806	assert(p!=NULL);
3807
3808	if (p->left == NULL \|\| p->isDated == YES)
3809	return 1;
3810	else
3811	{
3812	i += NDatedTips (p->left);
3813	i += NDatedTips (p->right);
3814	return i;
3815	}
3816	}
3817
3818
3819
3820
3821
3822	/* OrderTips: Order tips in a polytomous tree */
3823	void OrderTips (PolyTree *t)
3824	{
3825	int i, j;
3826	PolyNode p, q, r, pl, ql, rl;
3827
3828	/* label by minimum index */
3829	for (i=0; i<t->nNodes; i++)
3830	{
3831	p = t->allDownPass[i];
3832	if (p->left == NULL)
3833	{
3834	if (t->isRooted == NO && p->index == localOutGroup)
3835	p->x = -1;
3836	else
3837	p->x = p->index;
3838	}
3839	else
3840	{
3841	j = t->nNodes;
3842	for (q=p->left; q!=NULL; q=q->sib)
3843	{
3844	if (q->x < j)
3845	j = q->x;
3846	}
3847	p->x = j;
3848	}
3849	}
3850
3851	/* and rearrange */
3852	for (i=0; i<t->nNodes; i++)
3853	{
3854	p = t->allDownPass[i];
3855	if (p->left == NULL \|\| p->anc == NULL)
3856	continue;
3857	for (ql=NULL, q=p->left; q->sib!=NULL; ql=q, q=q->sib)
3858	{
3859	for (rl=q, r=q->sib; r!=NULL; rl=r, r=r->sib)
3860	{
3861	if (r->x < q->x)
3862	{
3863	if (ql == NULL)
3864	p->left = r;
3865	if (r == q->sib) /* swap adjacent q and r */
3866	{
3867	if (ql != NULL)
3868	ql->sib = r;
3869	pl = r->sib;
3870	r->sib = q;
3871	q->sib = pl;
3872	}
3873	else /* swap separated q and r */
3874	{
3875	if (ql != NULL)
3876	ql->sib = r;
3877	pl = r->sib;
3878	r->sib = q->sib;
3879	rl->sib = q;
3880	q->sib = pl;
3881	}
3882	pl = q;
3883	q = r;
3884	r = pl;
3885	}
3886	}
3887	}
3888	}
3889	GetPolyDownPass(t);
3890	}
3891
3892
3893
3894
3895
3896	/* PrintNodes: Print a list of tree nodes, pointers and length */
3897	void PrintNodes (Tree *t)
3898	{
3899	int i;
3900	TreeNode *p;
3901
3902	printf ("Node\tleft\tright\tanc\tlength\n");
3903	for (i=0; i<t->nNodes; i++)
3904	{
3905	p = &t->nodes[i];
3906	MrBayesPrint ("%d\t%d\t%d\t%d\t%f\t%f\n",
3907	p->index,
3908	p->left == NULL ? -1 : p->left->index,
3909	p->right == NULL ? -1 : p->right->index,
3910	p->anc == NULL ? -1 : p->anc->index,
3911	p->length,
3912	p->nodeDepth);
3913	}
3914
3915	if (t->root == NULL)
3916	MrBayesPrint ("root: NULL\n");
3917	else
3918	MrBayesPrint ("root: %d\n", t->root->index);
3919
3920	MrBayesPrint ("allDownPass:");
3921	for (i=0; i<t->nNodes; i++)
3922	{
3923	p = t->allDownPass[i];
3924	if (p!=NULL)
3925	MrBayesPrint (" %d", p->index);
3926	else
3927	MrBayesPrint (" NULL");
3928	}
3929	MrBayesPrint ("\nintDownPass: ");
3930	for (i=0; i<t->nIntNodes; i++)
3931	{
3932	p = t->intDownPass[i];
3933	if (p!=NULL)
3934	MrBayesPrint (" %d\t", p->index);
3935	else
3936	MrBayesPrint (" NULL\t");
3937	}
3938	MrBayesPrint ("\n");
3939	}
3940
3941
3942
3943
3944
3945	/* PrintPolyNodes: Print a list of polytomous tree nodes, pointers and length */
3946	void PrintPolyNodes (PolyTree *pt)
3947	{
3948	int i, j, k;
3949	PolyNode *p;
3950
3951	printf ("Node\tleft\tsib\tanc\tlength\tlabel\n");
3952	for (i=0; i<pt->memNodes; i++)
3953	{
3954	p = &pt->nodes[i];
3955	MrBayesPrint ("%d\t%d\t%d\t%d\t%f\t%s\n",
3956	p->index,
3957	p->left == NULL ? -1 : p->left->index,
3958	p->sib == NULL ? -1 : p->sib->index,
3959	p->anc == NULL ? -1 : p->anc->index,
3960	p->length,
3961	p->label);
3962	}
3963	MrBayesPrint ("root: %d\n", pt->root->index);
3964	fflush(stdout);
3965
3966	if (pt->nBSets > 0)
3967	{
3968	for (i=0; i<pt->nBSets; i++)
3969	{
3970	printf ("Effective branch length set '%s'\n", pt->bSetName[i]);
3971	for (j=0; j<pt->nNodes; j++)
3972	{
3973	printf ("%d:%f", j, pt->effectiveBrLen[pt->nBSets][j]);
3974	if (j != pt->nNodes-1)
3975	printf (", ");
3976	}
3977	printf ("\n");
3978	}
3979	}
3980
3981	if (pt->nESets > 0)
3982	{
3983	for (i=0; i<pt->nESets; i++)
3984	{
3985	printf ("Cpp event set '%s'\n", pt->eSetName[i]);
3986	for (j=0; j<pt->nNodes; j++)
3987	{
3988	if (pt->nEvents[i*pt->nNodes+j] > 0)
3989	{
3990	printf ("\tNode %d -- %d:(", j, pt->nEvents[i][j]);
3991	for (k=0; k<pt->nEvents[i][j]; k++)
3992	{
3993	printf ("%f %f", pt->position[i][j][k], pt->rateMult[i][j][k]);
3994	if (k != pt->nEvents[i][j]-1)
3995	printf (", ");
3996	}
3997	printf (")\n");
3998	}
3999	}
4000	printf ("\n");
4001	}
4002	}
4003
4004	fflush(stdout);
4005	}
4006
4007
4008
4009
4010
4011	/**
4012	Update relaxed clock paramiter of the branch of a node with index "b" after node with index "a" is removed.
4013	i.e. make branch of node with index "b" be a concatenation of its original branch and the branch of node with index "a"
4014	Relaxed clock paramiter of node with index "a" become invalid in the process.
4015	Note: For Non-clock models the routine has no effect.
4016
4017	\| \|
4018	\| \|
4019	a \|
4020	\| -> \|
4021	\| \|
4022	\| b
4023	b
4024
4025	*/
4026	void AppendRelaxedBranch (int a,int b,PolyTree *t){
4027
4028	int i,len;
4029
4030	for (i=0; i<t->nBSets; i++)
4031	{
4032	t->effectiveBrLen[i][b] += t->effectiveBrLen[i][a];
4033	}
4034
4035	for (i=0; i<t->nESets; i++)
4036	{
4037	len=t->nEvents[i][a]+t->nEvents[i][b];
4038	t->position[i][a] = (MrBFlt ) SafeRealloc ((void )t->position[i][a], len*sizeof(MrBFlt));
4039	t->rateMult[i][a] = (MrBFlt ) SafeRealloc ((void )t->rateMult[i][a], len*sizeof(MrBFlt));
4040	memcpy( t->position[i][a]+t->nEvents[i][a], t->position[i][b], t->nEvents[i][b]*sizeof(MrBFlt));
4041	memcpy( t->rateMult[i][a]+t->nEvents[i][a], t->rateMult[i][b], t->nEvents[i][b]*sizeof(MrBFlt));
4042	free(t->position[i][b]);
4043	free(t->rateMult[i][b]);
4044	t->position[i][b] = t->position[i][a];
4045	t->rateMult[i][b] = t->rateMult[i][a];
4046	t->position[i][a] = NULL;
4047	t->rateMult[i][a] = NULL;
4048	t->nEvents[i][a] = 0;
4049	t->nEvents[i][b] = len;
4050	}
4051
4052	}
4053
4054
4055
4056
4057
4058	/**
4059	Swap relaxed clock paramiters of the branch of nodes with index "a" and "b".
4060	*/
4061	void SwapRelaxedBranchInfo (int a,int b,PolyTree *t){
4062
4063	int i,j;
4064	MrBFlt tmp, *tmpp;
4065
4066	for (i=0; i<t->nBSets; i++)
4067	{
4068	tmp = t->effectiveBrLen[i][a];
4069	t->effectiveBrLen[i][a] = t->effectiveBrLen[i][b];
4070	t->effectiveBrLen[i][b] = tmp;
4071	}
4072	if (t->popSizeSet == YES)
4073	{
4074	tmp = t->popSize[a];
4075	t->popSize[a]=t->popSize[b];
4076	t->popSize[b] = tmp;
4077	}
4078
4079	for (i=0; i<t->nESets; i++)
4080	{
4081	tmpp = t->position[i][a];
4082	t->position[i][a] = t->position[i][b];
4083	t->position[i][b] = tmpp;
4084	tmpp = t->rateMult[i][a];
4085	t->rateMult[i][a] = t->rateMult[i][b];
4086	t->rateMult[i][b] = tmpp;
4087	j = t->nEvents[i][a];
4088	t->nEvents[i][a] = t->nEvents[i][b];
4089	t->nEvents[i][b] = j;
4090	}
4091	}
4092
4093
4094
4095
4096
4097	/*-------------------------------------------------------------------------------------------
4098	\|
4099	\| PrunePolyTree: This routine will prune a polytomous tree according to the currently
4100	\| included taxa. NB! All tree nodes cannot be accessed by cycling over the
4101	\| pt->nodes array after the deletion, because some spaces will be occupied by deleted
4102	\| nodes and pt->nNodes is no longer the length of this array
4103	\| (if proper re-arangment of pt->nodes needed then remove "#if 0" and make changes to p->x, see below).
4104	\|
4105	---------------------------------------------------------------------------------------------*/
4106	int PrunePolyTree (PolyTree *pt)
4107
4108	{
4109
4110	int i, j, numDeleted, numTermPruned, numIntPruned, index;
4111	PolyNode p = NULL, q=NULL, r=NULL, qa;
4112
4113	numDeleted = 0;
4114	for (i=0; i<pt->nNodes; i++)
4115	{
4116	p = pt->allDownPass[i];
4117	CheckString (taxaNames, numTaxa, p->label, &index);
4118	if (p->left == NULL && taxaInfo[index].isDeleted == YES)
4119	numDeleted++;
4120	}
4121
4122	if (numDeleted == 0)
4123	{
4124	/* nothing to do */
4125	return (NO_ERROR);
4126	}
4127	if (pt->nNodes - pt->nIntNodes - numDeleted < 3)
4128	{
4129	MrBayesPrint ("%s Pruned tree has less than three taxa in it\n", spacer);
4130	return (ERROR);
4131	}
4132	if (pt->nNodes - pt->nIntNodes < numLocalTaxa)
4133	{
4134	MrBayesPrint ("%s Tree to be pruned does not include all taxa\n", spacer);
4135	return (ERROR);
4136	}
4137
4138	/* prune away one node at a time */
4139	numIntPruned = 0;
4140	numTermPruned = 0;
4141	for (i=0; i<pt->nNodes; i++)
4142	{
4143	p = pt->allDownPass[i];
4144	if (p->left != NULL)
4145	continue;
4146	CheckString (taxaNames, numTaxa, p->label, &index);
4147	if (taxaInfo[index].isDeleted == YES)
4148	{
4149	numTermPruned++;
4150	for (q=p->anc->left; q!=NULL; q=q->sib)
4151	{
4152	if (q->sib == p)
4153	break;
4154	}
4155	if (q == NULL)
4156	{
4157	/* p is the left of its ancestor */
4158	assert( p->anc->left == p);
4159	p->anc->left = p->sib;
4160	}
4161	else
4162	{
4163	/* p is q->sib; this also works if p->sib is NULL */
4164	q->sib = p->sib;
4165	}
4166	/* if only one child left, delete ancestral node */
4167	j = 0;
4168	for (q=p->anc->left; q!=NULL; q=q->sib)
4169	j++;
4170	if (j == 1)
4171	{
4172	/* p->anc->left is only child left; make p->anc be p->anc->left and accommodate its length */
4173	numIntPruned++;
4174	qa= p->anc;
4175	q = qa->left;
4176	if (q->left == NULL)
4177	{
4178	AppendRelaxedBranch ( qa->index, q->index, pt );
4179	qa->index = q->index;
4180	qa->length += q->length;
4181	strcpy(qa->label, q->label);
4182	qa->left = NULL;
4183	/* To make sure that q is not treated as the representer of the tip it represented before. i.e. make condition if (p->left != NULL) true */
4184	q->left = (struct pNode*)1;
4185	}
4186	else
4187	{
4188	if (qa->anc != NULL)
4189	{
4190	AppendRelaxedBranch ( qa->index, q->index, pt );
4191	qa->length += q->length;
4192	}
4193	qa->index = q->index;
4194	qa->left = q->left;
4195	for(r=q->left; r!= NULL; r=r->sib)
4196	r->anc = qa;
4197	}
4198	}
4199	/* if unrooted, then root node has to have more then 2 children, thus the following check */
4200	if (j == 2 && pt->isRooted == NO && p->anc->anc == NULL)
4201	{
4202	numIntPruned++;
4203	r=p->anc; /r is the root with only 2 children/
4204	if ( r->left->left != NULL )
4205	{/* Make r->left new root by attaching "right" child of r to children of r->left */
4206	for (q=r->left->left; q->sib!=NULL; q=q->sib)
4207	;
4208	q->sib = r->left->sib;
4209	r->left->sib->anc = q->anc;
4210	r->left->sib->length += q->anc->length;
4211	r->left->sib = NULL;
4212	r->left->anc = NULL;
4213	pt->root = r->left;
4214	}
4215	else
4216	{/* Make "right" child of r (r->left->sib) the new root by attaching r->left to children of r->"right" */
4217	for ( q=r->left->sib->left; q->sib!=NULL; q=q->sib )
4218	;
4219	q->sib = r->left;
4220	r->left->anc = q->anc;
4221	r->left->length += q->anc->length;
4222	r->left->sib = NULL;
4223	q->anc->anc = NULL;
4224	pt->root = q->anc;
4225	}
4226	}
4227	}
4228	}
4229
4230	#if 0
4231	/* place unused space at the end of pt->nodes array. If activated this code p->x has to be set to non 0 value for all p that are deleted. */
4232	for (i=0; i<pt->nNodes; i++)
4233	{
4234	p = &pt->nodes[i];
4235	if (p->x != 0)
4236	{
4237	for (j=i+1; j<pt->nNodes; j++)
4238	{
4239	q = &pt->nodes[j];
4240	if (q->x == 0)
4241	break;
4242	}
4243	if (j != pt->nNodes)
4244	{
4245	/* swap nodes; quite difficult! */
4246	CopyPolyNodes (p, q, nLongsNeeded);
4247	p->left = q->left;
4248	p->sib = q->sib;
4249	p->anc = q->anc;
4250	for (k=0; k<pt->nNodes; k++)
4251	{
4252	r = &pt->nodes[k];
4253	if (r->left == q)
4254	r->left = p;
4255	if (r->sib == q)
4256	r->sib = p;
4257	if (r->anc == q)
4258	r->anc = p;
4259	}
4260	}
4261	}
4262	}
4263	#endif
4264
4265	/* correct number of nodes */
4266	pt->nNodes -= (numTermPruned + numIntPruned);
4267	pt->nIntNodes -= numIntPruned;
4268
4269	/* get downpass; note that the deletion procedure does not change the root in rooted case */
4270	i=j=0;
4271	GetPolyNodeDownPass ( pt, pt->root, &i, &j );
4272	assert(i==pt->nNodes);
4273	assert(j==pt->nIntNodes);
4274
4275	return (NO_ERROR);
4276
4277	}
4278
4279
4280
4281
4282
4283	/*--------------------------------------------------------------------
4284	\|
4285	\| RandPerturb: Randomly perturb a tree by nPert NNIs
4286	\|
4287	---------------------------------------------------------------------*/
4288	int RandPerturb (Tree t, int nPert, SafeLong seed)
4289	{
4290
4291	int i, whichNode;
4292	TreeNode p, q, a, b, *c;
4293
4294	if (t->nConstraints >= t->nIntNodes)
4295	{
4296	MrBayesPrint ("%s User tree cannot be perturbed because all nodes are locked\n", spacer);
4297	return (ERROR);
4298	}
4299
4300	for (i=0; i<nPert; i++)
4301	{
4302	do
4303	{
4304	whichNode = (int)(RandomNumber(seed) * (t->nIntNodes - 1));
4305	p = t->intDownPass[whichNode];
4306	} while (p->isLocked == YES);
4307
4308	q = p->anc;
4309	a = p->left;
4310	b = p->right;
4311	if (q->left == p)
4312	c = q->right;
4313	else
4314	c = q->left;
4315
4316	if (RandomNumber(seed) < 0.5)
4317	{
4318	/* swap b and c */
4319	p->right = c;
4320	c->anc = p;
4321
4322	if (q->left == c)
4323	q->left = b;
4324	else
4325	q->right = b;
4326	b->anc = q;
4327	}
4328	else
4329	{
4330	/* swap a and c */
4331	p->left = c;
4332	c->anc = p;
4333
4334	if (q->left == c)
4335	q->left = a;
4336	else
4337	q->right = a;
4338	a->anc = q;
4339	}
4340
4341	if (t->isCalibrated == YES)
4342	InitCalibratedBrlens (t, 0.0001, seed);
4343	else if (t->isClock == YES)
4344	InitClockBrlens (t);
4345	}
4346
4347	GetDownPass (t);
4348
4349	if (t->checkConstraints == YES && CheckConstraints (t) == NO_ERROR)
4350	{
4351	MrBayesPrint ("%s Broke constraints when perturbing tree\n", spacer);
4352	return (ERROR);
4353	}
4354
4355	return (NO_ERROR);
4356	}
4357
4358
4359
4360
4361
4362	/*
4363	\| Reorder array of nodes "nodeArray" such that first nodes in it could be paired with "w" to create imediat common ancestor and this ancesor node would not vialate any constraint.
4364	\|
4365	\| @param w Reference node as described
4366	\| @param nodeArray A set of node to order as described
4367	\| @param activeConstraints Array containing indeces of active constraints in the set of defined constraints
4368	\| @param nLongsNeeded Length of partition information (in SafeLong) in a node and constraint deffinition.
4369	\| @param isRooted Do constraints apply to rootet tree YES or NO
4370	\|
4371	\| @return Number of nodes in "nodeArray" that could be paired with "w" to create imediat common ancestor and this ancesor node would not vialate any constraint
4372	*/
4373	int ConstraintAllowedSet(PolyNode w, PolyNode nodeArray, int nodeArraySize, int activeConstraints, int activeConstraintsSize, int nLongsNeeded, int isRooted )
4374	{
4375
4376	int i, j, k, FirstEmpty;
4377	SafeLong **constraintPartition;
4378	PolyNode *tmp;
4379
4380	for (j=0; j<activeConstraintsSize; j++)
4381	{
4382	k=activeConstraints[j];
4383
4384	if( definedConstraintsType[k] == PARTIAL )
4385	{
4386	if( ( IsPartNested(definedConstraintPruned[k], w->partition, nLongsNeeded) == YES ) \|\|
4387	( isRooted == NO && IsPartNested(definedConstraintTwoPruned[k], w->partition, nLongsNeeded) == YES )
4388	)
4389	continue;/* all nodes are compartable because condition of the constraint has to be sutsfied in the subtree rooted at w*/
4390
4391	FirstEmpty = IsSectionEmpty(definedConstraintPruned[k], w->partition, nLongsNeeded);
4392	if( FirstEmpty == YES && IsSectionEmpty(definedConstraintTwoPruned[k], w->partition, nLongsNeeded) == YES )
4393	continue; /* all nodes are compartable becouse w does not contain any constraint taxa*/
4394
4395	assert(FirstEmpty^IsSectionEmpty(definedConstraintTwoPruned[k], w->partition, nLongsNeeded));
4396
4397	if( FirstEmpty == YES )
4398	{/w->partition has intersection with definedConstraintTwoPruned[k], thus remove all nodes from nodeArray that intersect with definedConstraintPruned[k]/
4399	constraintPartition=definedConstraintPruned;
4400	}
4401	else
4402	{/w->partition has intersection with definedConstraintPruned[k], thus remove all nodes from nodeArray that intersect with definedConstraintTwoPruned[k]/
4403	constraintPartition=definedConstraintTwoPruned;
4404	}
4405
4406	for(i=0;i<nodeArraySize;i++)
4407	{
4408	if( IsSectionEmpty(constraintPartition[k], nodeArray[i]->partition, nLongsNeeded) == NO &&
4409	( ( FirstEmpty == NO && isRooted== YES ) \|\| IsPartNested(constraintPartition[k], nodeArray[i]->partition, nLongsNeeded) == NO )
4410	)
4411	/*second part of if statment is to bail out "nodeArray[i]" when "w" contains nodes for example from definedConstraintPruned and "nodeArray[i]" have definedConstraintTwoPruned fully nested in it
4412	This bail out not applicable if t->isRooted== YES Since we should create a rooting node for the first set of taxa in the constraint.
4413	Note that such case possible because we may have hard constraint node that fully nest definedConstraintTwoPruned but also having taxa from definedConstraintPruned keeping constraint active.*/
4414	{
4415	tmp = nodeArray[i];
4416	nodeArray[i]=nodeArray[--nodeArraySize];
4417	nodeArray[nodeArraySize]=tmp;
4418	i--;
4419	}
4420	}
4421	}/end if PARTIAL/
4422	else
4423	{
4424	assert( definedConstraintsType[k] == NEGATIVE );
4425	if( isRooted == YES \|\| IsBitSet(localOutGroup, definedConstraintPruned[k])==NO )
4426	constraintPartition=definedConstraintPruned;
4427	else
4428	constraintPartition=definedConstraintTwoPruned;
4429
4430	if( IsSectionEmpty(constraintPartition[k], w->partition, nLongsNeeded)==YES )
4431	continue;
4432
4433	for(i=0;i<nodeArraySize;i++)
4434	{
4435	if( IsUnionEqThird (w->partition, nodeArray[i]->partition, constraintPartition[k], nLongsNeeded) == YES )
4436	{
4437	tmp = nodeArray[i];
4438	nodeArray[i]=nodeArray[--nodeArraySize];
4439	nodeArray[nodeArraySize]=tmp;
4440	i--;
4441	}
4442	}
4443
4444	}/end if NEGATIVE/
4445	}
4446
4447	return nodeArraySize;
4448	}
4449
4450
4451
4452
4453	/*
4454	\| Check if "partition" violate any constraint.
4455	\|
4456	\| @param partiton Partition to check
4457	\| @param activeConstraints Array containing indeces of active constraints in the set of defined constraints
4458	\| @param nLongsNeeded Length of partition information (in SafeLong) in a node and constraint deffinition
4459	\| @param isRooted Do constraints apply to rootet tree YES or NO
4460	\|
4461	\| @return Index of first violated constraint in activeConstraints array, -1 if no constraint is violated.
4462	*/
4463	int ViolatedConstraint(SafeLong partition, int activeConstraints, int activeConstraintsSize, int nLongsNeeded, int isRooted )
4464	{
4465
4466	int j, k;
4467	SafeLong **constraintPartition;
4468
4469
4470	for (j=0; j<activeConstraintsSize; j++)
4471	{
4472	k=activeConstraints[j];
4473	assert(definedConstraintsType[k] != HARD);
4474
4475	if( definedConstraintsType[k] == PARTIAL )
4476	{
4477	if( ( IsSectionEmpty(definedConstraintPruned[k], partition, nLongsNeeded) == NO ) &&
4478	( IsSectionEmpty(definedConstraintTwoPruned[k], partition, nLongsNeeded) == NO ) &&
4479	( IsPartNested(definedConstraintPruned[k], partition, nLongsNeeded) == NO) &&
4480	!( isRooted == NO && IsPartNested(definedConstraintTwoPruned[k], partition, nLongsNeeded) == YES)
4481	)
4482	return j;
4483
4484	}/end if PARTIAL/
4485	else
4486	{
4487	assert( definedConstraintsType[k] == NEGATIVE );
4488	if( isRooted == YES \|\| IsBitSet(localOutGroup, definedConstraintPruned[k])==NO )
4489	constraintPartition=definedConstraintPruned;
4490	else
4491	constraintPartition=definedConstraintTwoPruned;
4492
4493	if( IsUnionEqThird (partition, partition, constraintPartition[k], nLongsNeeded) == YES )
4494	return j;
4495
4496	}/end if NEGATIVE/
4497	}
4498
4499	return -1;
4500	}
4501
4502
4503
4504
4505
4506
4507	/*
4508	\| Remove from activeConstraints references to constraints that become satisfied if PolyNode "w" exist, i.e. they do not need to be checked furter thus become not active
4509	\|
4510	\| @param activeConstraints Array containing indeces of active constraints in the set of defined constraints
4511	\| @param nLongsNeeded Length of partition information (in SafeLong) in a node and constraint deffinition.
4512	\| @param isRooted Do constraints apply to rootet tree YES or NO
4513	\|
4514	\| @return Size of pruned "activeConstraints" array
4515	*/
4516	int PruneActiveConstraints (PolyNode w, int activeConstraints, int activeConstraintsSize, int nLongsNeeded, int isRooted )
4517	{
4518
4519	int j, k;
4520	SafeLong **constraintPartition;
4521	//PolyNode *tmp;
4522
4523	for (j=0; j<activeConstraintsSize; j++)
4524	{
4525	k=activeConstraints[j];
4526
4527	if( definedConstraintsType[k] == PARTIAL )
4528	{
4529	if( (IsPartNested(definedConstraintPruned[k], w->partition, nLongsNeeded) == YES && IsSectionEmpty(definedConstraintTwoPruned[k], w->partition, nLongsNeeded)) \|\|
4530	( isRooted == NO && IsPartNested(definedConstraintTwoPruned[k], w->partition, nLongsNeeded) == YES && IsSectionEmpty(definedConstraintPruned[k], w->partition, nLongsNeeded))
4531	)
4532	{
4533	//tmp = activeConstraints[j];
4534	activeConstraints[j]=activeConstraints[--activeConstraintsSize];
4535	//activeConstraints[activeConstraintsSize]=tmp;
4536	j--;
4537	}
4538	}/end if PARTIAL/
4539	else
4540	{
4541	assert( definedConstraintsType[k] == NEGATIVE );
4542	if( isRooted == YES \|\| IsBitSet(localOutGroup, definedConstraintPruned[k])==NO )
4543	constraintPartition=definedConstraintPruned;
4544	else
4545	constraintPartition=definedConstraintTwoPruned;
4546
4547	if( IsPartNested(constraintPartition[k], w->partition, nLongsNeeded)==NO && IsSectionEmpty(constraintPartition[k], w->partition, nLongsNeeded)==NO )
4548	{
4549	//tmp = activeConstraints[j];
4550	activeConstraints[j]=activeConstraints[--activeConstraintsSize];
4551	//activeConstraints[activeConstraintsSize]=tmp;
4552	j--;
4553	}
4554	}/end if NEGATIVE/
4555	}
4556
4557	return activeConstraintsSize;
4558	}
4559
4560
4561
4562
4563
4564	/*--------------------------------------------------------------------
4565	\|
4566	\| RandResolve: Randomly resolve a polytomous tree
4567	\|
4568	\| @param tt is a tree which contains information about applicable constraints. If it is set to NULL then no constraints will be used.
4569	If t!=NULL then partitions of nodes of polytree should be allocated for example by AllocatePolyTreePartitions (t);
4570	\| @return NO_ERROR on succes, ABORT if could not resolve a tree without vialating some consraint, ERROR if any other error occur
4571	---------------------------------------------------------------------*/
4572	int RandResolve (Tree tt, PolyTree t, SafeLong *seed, int destinationIsRooted)
4573
4574	{
4575
4576	int i, j, k, nextNode, stopNode, rand1, rand2, nLongsNeeded, tmp;
4577	PolyNode p=NULL, q, r, u, w1, w2;
4578	int nodeArrayAllowedSize, nodeArraySize, activeConstraintsSize;
4579	PolyNode **nodeArray;
4580	int *activeConstraints;
4581
4582	assert(tt==NULL \|\| t->bitsets!=NULL); /* partition fields of t nodes need to be allocated if constraints are used*/
4583	assert( numLocalTaxa <= t->memNodes/2); /* allocated tree has to be big enough*/
4584	nLongsNeeded = (numLocalTaxa - 1) / nBitsInALong + 1; /* allocated lenght of partitions is t->memNodes/2 bits but only first numLocalTaxa bits are used */
4585
4586	nodeArray = t->allDownPass; /temporary use t->allDownPass for different purpose. It get properly reset at the end. /
4587	activeConstraints = tempActiveConstraints;
4588	activeConstraintsSize = 0;
4589
4590	/* collect constraints to consider if applicable*/
4591	if( tt!=NULL && tt->constraints!=NULL )
4592	{
4593	for (k=0; k<numDefinedConstraints; k++)
4594	{
4595	if( tt->constraints[k] == YES && definedConstraintsType[k] != HARD )
4596	activeConstraints[activeConstraintsSize++]=k;
4597	}
4598	}
4599
4600	/* count immediate descendants */
4601	GetPolyDownPass(t);
4602	for (i=0; i<t->nIntNodes; i++)
4603	{
4604	p = t->intDownPass[i];
4605	tmp=ViolatedConstraint(p->partition, activeConstraints, activeConstraintsSize, nLongsNeeded, t->isRooted);
4606	if ( tmp != -1)
4607	{
4608	assert( p->isLocked == YES );
4609	MrBayesPrint ("%s Could not build a constraint tree since hard constraint \"%s\" and constraint \"%s\" are incompatible\n", spacer, constraintNames[p->lockID], constraintNames[activeConstraints[tmp]]);
4610	return (ERROR);
4611	}
4612	activeConstraintsSize = PruneActiveConstraints (p, activeConstraints, activeConstraintsSize, nLongsNeeded, t->isRooted);
4613	j = 0;
4614	for (q=p->left; q!=NULL; q=q->sib)
4615	j++;
4616	p->x = j;
4617	}
4618
4619	/* add one node at a time */
4620	if (destinationIsRooted == NO)
4621	stopNode = 2*numLocalTaxa - 2;
4622	else
4623	stopNode = 2*numLocalTaxa - 1;
4624	for (nextNode=t->nNodes; nextNode < stopNode; nextNode++)
4625	{
4626	/* find a polytomy to break */
4627	for (i=0; i<t->nIntNodes; i++)
4628	{
4629	p = t->intDownPass[i];
4630	if (destinationIsRooted == YES && p->x > 2)
4631	break;
4632	if (destinationIsRooted == NO && ((p->anc != NULL && p->x > 2) \|\| (p->anc == NULL && p->x > 3)))
4633	break;
4634	}
4635
4636	/* if we can't find one, there's an error */
4637	if (i == t->nIntNodes)
4638	{
4639	return ERROR;
4640	}
4641
4642	nodeArraySize=0;
4643	/Collect initial list of condidate nodes to join/
4644	for (q = p->left; q!= NULL; q = q->sib)
4645	{
4646	nodeArray[nodeArraySize++]=q;
4647	}
4648	assert(nodeArraySize==p->x);
4649
4650	/* identify two descendants randomly */
4651	/* make sure we do not select outgroup if it is an unrooted tree */
4652	if (p->anc == NULL && destinationIsRooted == NO)
4653	nodeArraySize--;
4654
4655	do
4656	{
4657	/Pick first node/
4658	rand1 = (int) (RandomNumber(seed) * nodeArraySize);
4659	w1 = nodeArray[rand1];
4660	nodeArray[rand1] = nodeArray[--nodeArraySize];
4661
4662	if( nodeArraySize==0 )
4663	return ABORT; /Potentaily here we could instead revert by removing last added node and try again. /
4664
4665	/Move all nodes in nodeArray which can be paired with w to the begining of array/
4666	nodeArrayAllowedSize=ConstraintAllowedSet(w1, nodeArray, nodeArraySize, activeConstraints, activeConstraintsSize, nLongsNeeded, t->isRooted );
4667	/TODO optimization for Maxim(if not Maxim remove it if you still see it): if nodeArrayAllowedSize==0 then set w1->y /
4668	}while( nodeArrayAllowedSize==0 );
4669
4670	rand2 = (int) (RandomNumber(seed) *nodeArrayAllowedSize);
4671	w2 = nodeArray[rand2];
4672
4673
4674	/* create a new node */
4675	u = &t->nodes[nextNode];
4676	u->anc = p;
4677	u->x = 2;
4678	p->x--;
4679
4680	if( tt!=NULL ){
4681	for (j=0; j<nLongsNeeded; j++)
4682	u->partition[j] = w1->partition[j] \| w2->partition[j] ;
4683	activeConstraintsSize = PruneActiveConstraints (u, activeConstraints, activeConstraintsSize, nLongsNeeded, t->isRooted );
4684	}
4685
4686	u->left = w1;
4687	t->nNodes++;
4688	t->nIntNodes++;
4689
4690
4691	/* connect tree together */
4692	r = u;
4693	for (q = p->left; q!= NULL; q = q->sib)
4694	{
4695	if ( q != w1 && q!=w2)
4696	{
4697	r->sib=q;
4698	r = q;
4699	}
4700	}
4701	r->sib = NULL;
4702	w1->sib = w2;
4703	w2->sib = NULL;
4704	w1->anc = u;
4705	w2->anc = u;
4706	p->left = u;
4707
4708	/* update tree */
4709	GetPolyDownPass (t);
4710	}
4711
4712	/* relabel interior nodes (important that last indices are at the bottom!) */
4713	for (i=0; i<t->nIntNodes; i++)
4714	{
4715	p = t->intDownPass[i];
4716	p->index = numLocalTaxa + i;
4717	}
4718	return NO_ERROR;
4719	}
4720
4721
4722
4723
4724
4725	/* ResetTreeNode: Reset tree node except for memory index */
4726	void ResetTreeNode (TreeNode *p)
4727	{
4728	/* do not change memoryIndex; that is set once and for all when tree is allocated */
4729	p->index = 0;
4730	p->scalerNode = NO;
4731	p->upDateCl = NO;
4732	p->upDateTi = NO;
4733	p->marked = NO;
4734	p->length = 0.0;
4735	p->nodeDepth = 0.0;
4736	p->x = 0;
4737	p->y = 0;
4738	p->index = 0;
4739	p->isDated = NO;
4740	p->calibration = NULL;
4741	p->age = -1.0;
4742	p->isLocked = NO;
4743	p->lockID = -1;
4744	p->label = noLabel;
4745	p->d = 0.0;
4746	p->partition = NULL;
4747	}
4748
4749
4750
4751
4752
4753	/* ResetPolyNode: Reset all values of one node in a polytree */
4754	void ResetPolyNode (PolyNode *p)
4755	{
4756	/* we reset everything here except memoryIndex, which should be immutable */
4757	p->length = 0.0;
4758	p->depth = 0.0;
4759	p->age = 0.0;
4760	p->anc = p->left = p->sib = NULL;
4761	p->calibration = NULL;
4762	p->f = 0.0;
4763	p->index = 0;
4764	p->isDated = NO;
4765	p->isLocked = NO;
4766	strcpy (p->label,"");
4767	p->lockID = 0;
4768	p->partition = NULL;
4769	p->partitionIndex = 0;
4770	p->support = 0.0;
4771	p->x = p->y = 0;
4772	}
4773
4774
4775
4776
4777
4778	/* ResetPolyTree: Reset polytomous tree to pristine state but keep relevant memory. */
4779	void ResetPolyTree (PolyTree *pt)
4780	{
4781	int i, maxTaxa, nLongsNeeded;
4782
4783	/* clear nodes */
4784	for (i=0; i<pt->memNodes; i++)
4785	ResetPolyNode (&pt->nodes[i]);
4786
4787	/* empty node arrays and tree properties but keep space */
4788	for (i=0; i<pt->nNodes; i++)
4789	pt->allDownPass[i] = NULL;
4790	for (i=0; i<pt->nIntNodes; i++)
4791	pt->intDownPass[i] = NULL;
4792	pt->nNodes = 0;
4793	pt->nIntNodes = 0;
4794	pt->root = NULL;
4795	pt->brlensDef = NO;
4796	pt->isRooted = NO;
4797	pt->isClock = NO;
4798	pt->isRelaxed = NO;
4799	pt->clockRate = 0.0;
4800
4801	/* empty bitsets but keep space and pointers */
4802	if (pt->bitsets)
4803	{
4804	maxTaxa = pt->memNodes / 2;
4805	nLongsNeeded = (maxTaxa - 1) / nBitsInALong + 1;
4806	for (i=0; i<pt->memNodes*nLongsNeeded; i++)
4807	pt->bitsets[i] = 0;
4808	for (i=0; i<pt->memNodes; i++)
4809	pt->nodes[i].partition = pt->bitsets + i*nLongsNeeded;
4810	}
4811
4812	/* empty relaxed clock parameters */
4813	FreePolyTreeRelClockParams (pt);
4814
4815	/* empty population size set parameters */
4816	FreePolyTreePopSizeParams (pt);
4817	}
4818
4819
4820
4821
4822
4823	/* ResetPolyTreePartitions: Reset and set bit patterns describing partitions */
4824	void ResetPolyTreePartitions (PolyTree *pt)
4825	{
4826	int i, j, numTaxa, nLongsNeeded;
4827	PolyNode *pp;
4828
4829	/* get some handy numbers */
4830	numTaxa = pt->memNodes/2;
4831	nLongsNeeded = (numTaxa - 1) / nBitsInALong + 1;
4832
4833	/* reset bits describing partitions */
4834	for (i=0; i<pt->memNodes*nLongsNeeded; i++)
4835	{
4836	pt->bitsets[i] = 0;
4837	}
4838
4839	/* set bits describing partitions */
4840	for (i=0; i<pt->nNodes; i++)
4841	{
4842	assert (pt->allDownPass != NULL && pt->allDownPass[i] != NULL);
4843	assert (pt->allDownPass[i]->partition != NULL);
4844
4845	pp = pt->allDownPass[i];
4846	if (pp->left == NULL)
4847	{
4848	SetBit (pp->index, pp->partition);
4849	}
4850	if (pp->anc != NULL)
4851	{
4852	for (j=0; j<nLongsNeeded; j++)
4853	pp->anc->partition[j] \|= pp->partition[j];
4854	}
4855	}
4856	}
4857
4858
4859
4860
4861
4862	/*----------------------------------------------
4863	\|
4864	\| ResetRootHeight: Reset node heights in a clock
4865	\| tree to fit a new root height. Assumes
4866	\| node depths and lengths set correctly.
4867	\|
4868	-----------------------------------------------*/
4869	int ResetRootHeight (Tree *t, MrBFlt rootHeight)
4870	{
4871	int i;
4872	TreeNode *p;
4873	MrBFlt factor, x, y;
4874
4875	if (t->isClock == NO)
4876	return ERROR;
4877
4878	/* make sure node depths are set */
4879	for (i=0; i<t->nNodes-1; i++)
4880	{
4881	p = t->allDownPass[i];
4882	if (p->left == NULL)
4883	p->nodeDepth = 0.0;
4884	else
4885	{
4886	x = p->left->nodeDepth + p->left->length;
4887	y = p->right->nodeDepth + p->right->length;
4888	if (x > y)
4889	p->nodeDepth = x;
4890	else
4891	p->nodeDepth = y;
4892	}
4893	}
4894	for (i=t->nNodes-3; i>=0; i--)
4895	{
4896	p = t->allDownPass[i];
4897	p->nodeDepth = p->anc->nodeDepth - p->length;
4898	}
4899
4900	/* now reset node depths and branch lengths */
4901	factor = rootHeight / t->root->left->nodeDepth;
4902	t->root->left->nodeDepth = rootHeight;
4903	for (i=t->nNodes-2; i>=0; i--)
4904	{
4905	p = t->allDownPass[i];
4906	p->nodeDepth *= factor;
4907	p->length *= factor;
4908	}
4909
4910	return NO_ERROR;
4911	}
4912
4913
4914
4915
4916
4917	/*----------------------------------------------
4918	\|
4919	\| ResetTipIndices: reset tip indices to be from
4920	\| 0 to number of included taxa, in same order
4921	\| as in the original taxon set.
4922	\|
4923	-----------------------------------------------*/
4924	void ResetTipIndices(PolyTree *pt)
4925	{
4926	int i, j, k, m;
4927	PolyNode *p;
4928
4929
4930	for (i=j=0; i<numTaxa; i++)
4931	{
4932	for (k=0; k<pt->nNodes; k++)
4933	{
4934	p = pt->allDownPass[k];
4935	if (StrCmpCaseInsensitive(p->label,taxaNames[i]) == 0)
4936	break;
4937	}
4938	if (k < pt->nNodes)
4939	{
4940	assert (p->left == NULL);
4941	if(p->index!=j){
4942	SwapRelaxedBranchInfo( p->index, j, pt);
4943	for(m=0; m<pt->nNodes; m++)
4944	{
4945	if(pt->allDownPass[m]->index==j)
4946	{
4947	pt->allDownPass[m]->index=p->index;
4948	break;
4949	}
4950	}
4951	p->index = j;
4952	}
4953	j++;
4954	}
4955	}
4956	}
4957
4958
4959
4960
4961
4962	/*----------------------------------------------
4963	\|
4964	\| ResetTopology: rebuild the tree t to fit the
4965	\| Newick string s. Everyting except topology
4966	\| is left in the same state in t.
4967	\|
4968	-----------------------------------------------*/
4969	int ResetTopology (Tree t, char s)
4970	{
4971	TreeNode p, q;
4972	int i, j, k, inLength;
4973	char temp[30];
4974
4975	/* set all pointers to NULL */
4976	for (i=0; i<t->memNodes; i++)
4977	{
4978	p = &t->nodes[i];
4979	p->anc = p->right = p->left = NULL;
4980	p->index = -1;
4981	}
4982	p = &t->nodes[0];
4983
4984	/* start out assuming that the tree is rooted; we will detect below if it is not */
4985	t->isRooted = YES;
4986	inLength = NO;
4987	for (i=0, j=1; *s!='\0'; s++)
4988	{
4989	if (s == ',' \|\| s == ')' \|\| *s == ':')
4990	{
4991	if (p->right == NULL && inLength == NO)
4992	{
4993	temp[i] = '\0';
4994	k = atoi (temp);
4995	p->index = k-1;
4996	i = 0;
4997	}
4998	else
4999	inLength = NO;
5000	}
5001	if (*s == '(')
5002	{
5003	q = p;
5004	p = &t->nodes[j++];
5005	q->left = p;
5006	p->anc = q;
5007	}
5008	else if (*s == ',')
5009	{
5010	if (p->anc->right == NULL)
5011	{
5012	q = p->anc;
5013	p = &t->nodes[j++];
5014	p->anc = q;
5015	q->right = p;
5016	}
5017	else /* if p->anc->right == p (near 'root' of unrooted trees) */
5018	{
5019	q = p->anc;
5020	p = &t->nodes[j++];
5021	q->anc = p;
5022	p->left = q;
5023	t->isRooted = NO;
5024	}
5025	}
5026	else if (*s == ')')
5027	{
5028	p = p->anc;
5029	}
5030	else if (*s == ':')
5031	{
5032	inLength = YES;
5033	}
5034	else if (inLength == NO)
5035	{
5036	temp[i++] = *s;
5037	}
5038	}
5039
5040	/* attach root to rooted tree */
5041	if (t->isRooted == YES)
5042	{
5043	p = &t->nodes[0];
5044	q = &t->nodes[j++];
5045	q->left = p;
5046	p->anc = q;
5047	}
5048
5049	/* relabel interior nodes, find number of nodes and root */
5050	t->nNodes = j;
5051	t->nIntNodes = t->nNodes/2 - 1;
5052
5053	if (t->isRooted == NO)
5054	j = t->nNodes - t->nIntNodes;
5055	else
5056	j = t->nNodes - t->nIntNodes - 1;
5057
5058	for (i=0; i<t->nNodes; i++)
5059	{
5060	p = &t->nodes[i];
5061	if (p->index == -1)
5062	p->index = j++;
5063	if (p->anc == NULL)
5064	t->root = p;
5065	}
5066
5067	GetDownPass (t);
5068
5069	return NO_ERROR;
5070	}
5071
5072
5073
5074
5075
5076	/*-----------------------------------------------------------------
5077	\|
5078	\| ResetBrlensFromTree: copies brlens and depths from second tree (vTree) to
5079	\| first tree (used to initialize brlen sets for same topology)
5080	\|
5081	-----------------------------------------------------------------*/
5082	int ResetBrlensFromTree (Tree tree, Tree vTree)
5083
5084	{
5085
5086	int i, j, k, nLongsNeeded, numTips;
5087	MrBFlt d1, d2;
5088	TreeNode p, q;
5089
5090	if (tree->isRooted != vTree->isRooted)
5091	return (ERROR);
5092
5093	if (AreTopologiesSame (tree, vTree) == NO)
5094	return (ERROR);
5095
5096	/* allocate and set up partitions */
5097	AllocateTreePartitions (vTree);
5098	AllocateTreePartitions (tree);
5099	numTips = tree->nNodes - tree->nIntNodes - (tree->isRooted == YES ? 1 : 0);
5100	nLongsNeeded = (int) ((numTips - 1) / nBitsInALong) + 1;
5101
5102	/copy lengths and nodeDepthes/
5103	for (i=0; i<vTree->nNodes; i++)
5104	{
5105	p = vTree->allDownPass[i];
5106	for (j=0; j<tree->nNodes; j++)
5107	{
5108	q = tree->allDownPass[j];
5109	for (k=0; k<nLongsNeeded; k++)
5110	if (p->partition[k] != q->partition[k])
5111	break;
5112	if (k==nLongsNeeded)
5113	{
5114	q->length = p->length;
5115	if (tree->isRooted == YES)
5116	q->nodeDepth = p->nodeDepth;
5117	}
5118	}
5119	}
5120
5121	if (tree->isRooted == YES)
5122	{
5123	/Next compute height for the root. /
5124	for (i=0; i<tree->nNodes-1; i++)
5125	{
5126	p = tree->allDownPass[i];
5127	if (p->left == NULL)
5128	p->nodeDepth = 0.0;
5129	else
5130	{
5131	d1 = p->left->nodeDepth + p->left->length;
5132	d2 = p->right->nodeDepth + p->right->length;
5133	if (d1 > d2)
5134	p->nodeDepth = d1;
5135	else
5136	p->nodeDepth = d2;
5137	}
5138	}
5139	for (i=tree->nNodes-3; i>=0; i--)
5140	{
5141	p = tree->allDownPass[i];
5142	if (p->left==NULL && p->calibration==NULL)
5143	continue; /* leave at 0.0 */
5144	p->nodeDepth = p->anc->nodeDepth - p->length;
5145	}
5146	}
5147
5148	FreeTreePartitions(tree);
5149	FreeTreePartitions(vTree);
5150
5151	return (NO_ERROR);
5152
5153	}
5154
5155
5156
5157
5158
5159	/* ResetIntNodeIndices: Set int node indices in downpass order from numTaxa to 2numTaxa-2 /
5160	void ResetIntNodeIndices (PolyTree *t)
5161	{
5162	int i, m, index;
5163
5164	index = t->nNodes - t->nIntNodes;
5165
5166	for (i=0; i<t->nIntNodes; i++)
5167	{
5168	if(t->intDownPass[i]->index != index)
5169	{
5170	SwapRelaxedBranchInfo( t->intDownPass[i]->index, index, t);
5171	for(m=0; m<t->nIntNodes; m++)
5172	{
5173	if(t->intDownPass[m]->index==index)
5174	{
5175	t->intDownPass[m]->index=t->intDownPass[i]->index;
5176	break;
5177	}
5178	}
5179	t->intDownPass[i]->index = index;
5180	}
5181	index++;
5182	}
5183	}
5184
5185
5186
5187
5188
5189	/* ResetTopologyFromTree: use top to set topology in tree */
5190	int ResetTopologyFromTree (Tree tree, Tree top)
5191	{
5192	int i, j, k;
5193	TreeNode p, q, r, p1;
5194
5195	/* adopt rooting */
5196	tree->isRooted = top->isRooted;
5197	tree->nNodes = top->nNodes;
5198	tree->nIntNodes = top->nIntNodes;
5199
5200	/* set all pointers to NULL */
5201	for (i=0; i<tree->nNodes; i++)
5202	{
5203	p = &tree->nodes[i];
5204	p->anc = p->right = p->left = NULL;
5205	}
5206
5207	/* now copy topology */
5208	for (i=0; i<top->nIntNodes; i++)
5209	{
5210	p1 = top->intDownPass[i];
5211
5212	k = p1->index;
5213	for (j=0; j<tree->nNodes; j++)
5214	if (tree->nodes[j].index == k)
5215	break;
5216	p = &tree->nodes[j];
5217
5218	k = p1->left->index;
5219	for (j=0; j<tree->nNodes; j++)
5220	if (tree->nodes[j].index == k)
5221	break;
5222	q = &tree->nodes[j];
5223
5224	k = p1->right->index;
5225	for (j=0; j<tree->nNodes; j++)
5226	if (tree->nodes[j].index == k)
5227	break;
5228	r = &tree->nodes[j];
5229
5230	p->left = q;
5231	p->right= r;
5232	q->anc = r->anc = p;
5233	}
5234
5235	/* arrange the root */
5236	k = top->root->index;
5237	for (j=0; j<tree->nNodes; j++)
5238	if (tree->nodes[j].index == k)
5239	break;
5240	p = &tree->nodes[j];
5241
5242	k = top->root->left->index;
5243	for (j=0; j<tree->nNodes; j++)
5244	if (tree->nodes[j].index == k)
5245	break;
5246	q = &tree->nodes[j];
5247	p->left = q;
5248	q->anc = p;
5249	p->right = p->anc = NULL;
5250	tree->root = p;
5251
5252	GetDownPass(tree);
5253
5254	return (NO_ERROR);
5255	}
5256
5257
5258
5259
5260
5261	/* ResetTopologyFromPolyTree: use polytree top to set topology in tree */
5262	int ResetTopologyFromPolyTree (Tree tree, PolyTree top)
5263	{
5264	int i, j, k;
5265	TreeNode p, q, *r;
5266	PolyNode *p1;
5267
5268	if (tree->isRooted != top->isRooted)
5269	return (ERROR);
5270
5271	/* set all pointers to NULL */
5272	for (i=0; i<tree->nNodes; i++)
5273	{
5274	p = &tree->nodes[i];
5275	p->anc = p->right = p->left = NULL;
5276	}
5277
5278	/* now copy topology */
5279	for (i=0; i<top->nIntNodes; i++)
5280	{
5281	p1 = top->intDownPass[i];
5282
5283	k = p1->index;
5284	for (j=0; j<tree->nNodes; j++)
5285	if (tree->nodes[j].index == k)
5286	break;
5287	p = &tree->nodes[j];
5288
5289	k = p1->left->index;
5290	for (j=0; j<tree->nNodes; j++)
5291	if (tree->nodes[j].index == k)
5292	break;
5293	q = &tree->nodes[j];
5294
5295	k = p1->left->sib->index;
5296	for (j=0; j<tree->nNodes; j++)
5297	if (tree->nodes[j].index == k)
5298	break;
5299	r = &tree->nodes[j];
5300
5301	p->left = q;
5302	p->right= r;
5303	q->anc = r->anc = p;
5304	}
5305
5306	/* arrange the root */
5307	if (top->isRooted == YES)
5308	{
5309	k = top->root->index;
5310	for (j=0; j<tree->nNodes; j++)
5311	if (tree->nodes[j].index == k)
5312	break;
5313	p = &tree->nodes[j];
5314
5315	k = top->nNodes;
5316	for (j=0; j<tree->nNodes; j++)
5317	if (tree->nodes[j].index == k)
5318	break;
5319	q = &tree->nodes[j];
5320
5321	q->left = p;
5322	q->anc = NULL;
5323	q->right = NULL;
5324	tree->root = q;
5325	}
5326	else /* if (top->isRooted == NO) */
5327	{
5328	k = top->root->index;
5329	for (j=0; j<tree->nNodes; j++)
5330	if (tree->nodes[j].index == k)
5331	break;
5332	p = &tree->nodes[j];
5333
5334	k = localOutGroup;
5335	for (p1=top->root->left; p1!=NULL; p1=p1->sib)
5336	if (p1->index == k)
5337	break;
5338
5339	assert (p1 != NULL);
5340	if (p1 == NULL)
5341	return (ERROR);
5342
5343	q = &tree->nodes[p1->index];
5344	k = p1->anc->left->sib->sib->index; /* index of missing child */
5345	if (p->left == q)
5346	p->left = &tree->nodes[k];
5347	else if (p->right == q)
5348	p->right = &tree->nodes[k];
5349
5350	q->anc = q->right = NULL;
5351	p->anc = q;
5352	q->left = p;
5353	}
5354
5355	GetDownPass(tree);
5356
5357	return (NO_ERROR);
5358	}
5359
5360
5361
5362
5363
5364	/* ResetTreePartitions: Reset bitsets describing tree partitions */
5365	void ResetTreePartitions (Tree *t)
5366	{
5367	int i, j, numTaxa, nLongsNeeded;
5368	TreeNode *p;
5369
5370	/* get some handy numbers */
5371	numTaxa = t->nNodes - t->nIntNodes - (t->isRooted == YES ? 1 : 0);
5372	nLongsNeeded = (numTaxa - 1) / nBitsInALong + 1;
5373
5374	/* reset bits describing partitions */
5375	for (i=0; i<t->nNodes; i++)
5376	{
5377	assert (t->allDownPass != NULL && t->allDownPass[i] != NULL);
5378	assert (t->allDownPass[i]->partition != NULL);
5379
5380	p = t->allDownPass[i];
5381	for (j=0; j<nLongsNeeded; j++)
5382	p->partition[j] = 0;
5383	}
5384
5385	/* set bits describing partitions */
5386	for (i=0; i<t->nNodes; i++)
5387	{
5388	p = t->allDownPass[i];
5389	if (p->left == NULL \|\| (p->anc == NULL && t->isRooted == NO))
5390	SetBit (p->index, p->partition);
5391	else if (p->anc != NULL)
5392	{
5393	for (j=0; j<nLongsNeeded; j++)
5394	p->partition[j] = p->left->partition[j] \| p->right->partition[j];
5395	}
5396	}
5397	}
5398
5399
5400
5401
5402
5403	/*-------------------------------------------------------
5404	\|
5405	\| RetrieveRTopology: This routine will rebuild a rooted
5406	\| tree from the order array created by StoreRTopology.
5407	\| All tree information except the structure will
5408	\| remain unaltered.
5409	\|
5410	--------------------------------------------------------*/
5411	int RetrieveRTopology (Tree t, int order)
5412
5413	{
5414	int i, numTaxa;
5415	TreeNode p, q, *r;
5416
5417	numTaxa = t->nNodes - t->nIntNodes - 1;
5418
5419	/* sort the tips in the t->allDownPass array */
5420	p = t->nodes;
5421	for (i=0; i<t->nNodes; i++, p++)
5422	t->allDownPass[p->index] = p;
5423
5424	/* make sure the root has index 2numTaxa-1 /
5425	q = t->allDownPass[t->nNodes-1];
5426	q->anc = q->right = NULL;
5427	t->root = q;
5428
5429	/* connect the first two tips */
5430	p = t->allDownPass[numTaxa];
5431	p->anc = q;
5432	q->left = p;
5433	p->length = 0.0;
5434	q = t->allDownPass[0];
5435	r = t->allDownPass[1];
5436	p->left = q;
5437	p->right = r;
5438	q->anc = r->anc = p;
5439
5440	/* add one tip at a time */
5441	for (i=2; i<numTaxa; i++)
5442	{
5443	p = t->allDownPass[i];
5444	q = t->allDownPass[numTaxa-1+i];
5445	r = t->allDownPass[*(order++)];
5446	p->anc = q;
5447	q->left = p;
5448	q->right = r;
5449	q->anc = r->anc;
5450	if (r->anc->left == r)
5451	r->anc->left = q;
5452	else
5453	r->anc->right = q;
5454	r->anc = q;
5455	}
5456
5457	/* get downpass */
5458	GetDownPass (t);
5459
5460	/* relabel interior nodes (root is correctly labeled already) */
5461	for (i=0; i<t->nIntNodes; i++)
5462	t->intDownPass[i]->index = i+numTaxa;
5463
5464	return (NO_ERROR);
5465	}
5466
5467
5468
5469
5470
5471	/*-------------------------------------------------------
5472	\|
5473	\| RetrieveRTree: This routine will rebuild a rooted
5474	\| tree from the arrays created by StoreRTree.
5475	\| All tree information except the structure and
5476	\| branch lengths will remain unaltered.
5477	\|
5478	--------------------------------------------------------*/
5479	int RetrieveRTree (Tree t, int order, MrBFlt *brlens)
5480
5481	{
5482	int i, numTaxa;
5483	TreeNode p, q, *r;
5484
5485	numTaxa = t->nNodes - t->nIntNodes - 1;
5486
5487	/* sort the tips in the t->allDownPass array */
5488	p = t->nodes;
5489	for (i=0; i<t->nNodes; i++, p++)
5490	t->allDownPass[p->index] = p;
5491
5492	/* make sure that root has index 2numTaxa-1 /
5493	q = t->allDownPass[t->nNodes-1];
5494	q->anc = q->right = NULL;
5495	q->length = 0.0;
5496	t->root = q;
5497
5498	/* connect the first three tips */
5499	p = t->allDownPass[numTaxa];
5500	p->anc = q;
5501	q->left = p;
5502	p->length = 0.0;
5503	q = t->allDownPass[0];
5504	r = t->allDownPass[1];
5505	p->left = q;
5506	p->right = r;
5507	q->anc = r->anc = p;
5508	q->length = *(brlens++);
5509	r->length = *(brlens++);
5510
5511	/* add one tip at a time */
5512	for (i=2; i<numTaxa; i++)
5513	{
5514	p = t->allDownPass[i];
5515	q = t->allDownPass[numTaxa-1+i];
5516	r = t->allDownPass[*(order++)];
5517	p->anc = q;
5518	q->left = p;
5519	q->right = r;
5520	q->anc = r->anc;
5521	if (r->anc->left == r)
5522	r->anc->left = q;
5523	else
5524	r->anc->right = q;
5525	r->anc = q;
5526	if (q->anc->anc != NULL)
5527	q->length = *(brlens++);
5528	else
5529	{
5530	r->length = *(brlens++);
5531	q->length = 0.0;
5532	}
5533	p->length = *(brlens++);
5534	}
5535
5536	/* get downpass */
5537	GetDownPass (t);
5538
5539	/* relabel interior nodes (root is correctly labeled already) */
5540	for (i=0; i<t->nIntNodes; i++)
5541	t->intDownPass[i]->index = i+numTaxa;
5542
5543	/* set the node depths */
5544	SetNodeDepths (t);
5545
5546	return (NO_ERROR);
5547	}
5548
5549
5550
5551
5552
5553	/*-------------------------------------------------------
5554	\|
5555	\| RetrieveRTreeWithIndices: This routine will rebuild a rooted
5556	\| tree from the arrays created by StoreRTreeWithIndices.
5557	\| All tree information except the structure, branch lengths
5558	\| and node indices will remain unaltered.
5559	\|
5560	--------------------------------------------------------*/
5561	int RetrieveRTreeWithIndices (Tree t, int order, MrBFlt *brlens)
5562
5563	{
5564	int i, numTaxa;
5565	TreeNode p, q, *r;
5566
5567	extern void ShowNodes (TreeNode *, int, int);
5568
5569	numTaxa = t->nNodes - t->nIntNodes - 1;
5570
5571	/* sort the tips in the t->allDownPass array */
5572	p = t->nodes;
5573	for (i=0; i<t->nNodes; i++, p++)
5574	t->allDownPass[p->index] = p;
5575
5576	/* make sure that root has index 2numTaxa-1 /
5577	q = t->allDownPass[t->nNodes-1];
5578	q->anc = q->right = NULL;
5579	q->length = 0.0;
5580	t->root = q;
5581
5582	/* connect the first three 'tips' with interior node, index from order array */
5583	p = t->allDownPass[numTaxa];
5584	p->x = *(order++);
5585	p->anc = q;
5586	q->left = p;
5587	p->length = 0.0;
5588	q = t->allDownPass[0];
5589	r = t->allDownPass[1];
5590	p->left = q;
5591	p->right = r;
5592	q->anc = r->anc = p;
5593	q->length = *(brlens++);
5594	r->length = *(brlens++);
5595
5596	/* add one tip at a time */
5597	for (i=2; i<numTaxa; i++)
5598	{
5599	p = t->allDownPass[i];
5600	assert (order >= numTaxa && order < 2*numTaxa - 1);
5601	q = t->allDownPass[numTaxa-1+i];
5602	q->x = *(order++);
5603	r = t->allDownPass[*(order++)];
5604	p->anc = q;
5605	q->left = p;
5606	q->right = r;
5607	q->anc = r->anc;
5608	if (r->anc->left == r)
5609	r->anc->left = q;
5610	else
5611	r->anc->right = q;
5612	r->anc = q;
5613	if (q->anc->anc != NULL)
5614	q->length = *(brlens++);
5615	else
5616	{
5617	r->length = *(brlens++);
5618	q->length = 0.0;
5619	}
5620	p->length = *(brlens++);
5621	}
5622
5623	/* get downpass */
5624	GetDownPass (t);
5625
5626	/* relabel interior nodes using labels in scratch variable x */
5627	for (i=0; i<t->nIntNodes; i++)
5628	{
5629	p = t->intDownPass[i];
5630	p->index = p->x;
5631	}
5632
5633	/* set the node depths */
5634	SetNodeDepths (t);
5635
5636	return (NO_ERROR);
5637	}
5638
5639
5640
5641
5642
5643	/*-------------------------------------------------------
5644	\|
5645	\| RetrieveUTopology: This routine will rebuild an unrooted
5646	\| tree from the order array created by StoreUTopology.
5647	\| All tree information except the structure
5648	\| will remain unaltered.
5649	\|
5650	--------------------------------------------------------*/
5651	int RetrieveUTopology (Tree t, int order)
5652
5653	{
5654	int i, numTips;
5655	TreeNode p, q, *r;
5656
5657	/* preliminaries */
5658	numTips = t->nNodes - t->nIntNodes;
5659	for (i=0; i<t->nNodes; i++)
5660	t->nodes[i].left = t->nodes[i].right = t->nodes[i].anc = NULL;
5661
5662	/* sort the tips in the t->allDownPass array */
5663	p = t->nodes;
5664	for (i=0; i<t->nNodes; i++, p++)
5665	t->allDownPass[p->index] = p;
5666
5667	/* make sure root has index 0 */
5668	q = t->allDownPass[0];
5669	q->anc = q->right = NULL;
5670	t->root = q;
5671
5672	/* connect the first three tips */
5673	p = t->allDownPass[numTips];
5674	p->anc = q;
5675	q->left = p;
5676	q = t->allDownPass[1];
5677	r = t->allDownPass[2];
5678	p->left = q;
5679	p->right = r;
5680	q->anc = r->anc = p;
5681
5682	/* add one tip at a time */
5683	for (i=3; i<numTips; i++)
5684	{
5685	p = t->allDownPass[i];
5686	q = t->allDownPass[numTips-2+i];
5687	r = t->allDownPass[order[i-3]];
5688	p->anc = q;
5689	q->left = p;
5690	q->right = r;
5691	q->anc = r->anc;
5692	if (r->anc->left == r)
5693	r->anc->left = q;
5694	else
5695	r->anc->right = q;
5696	r->anc = q;
5697	}
5698
5699	/* get downpass */
5700	GetDownPass (t);
5701
5702	/* relabel interior nodes (root is correctly labeled already) */
5703	for (i=0; i<t->nIntNodes; i++)
5704	t->intDownPass[i]->index = i+numTips;
5705
5706	return (NO_ERROR);
5707	}
5708
5709
5710
5711
5712
5713	/*-------------------------------------------------------
5714	\|
5715	\| RetrieveUTree: This routine will rebuild an unrooted
5716	\| tree from the arrays created by StoreUTree.
5717	\| All tree information except the structure and
5718	\| branch lengths will remain unaltered.
5719	\|
5720	--------------------------------------------------------*/
5721	int RetrieveUTree (Tree t, int order, MrBFlt *brlens)
5722
5723	{
5724	int i, numTips;
5725	TreeNode p, q, *r;
5726
5727	/* preliminaries */
5728	numTips = t->nNodes - t->nIntNodes;
5729	for (i=0; i<t->nNodes; i++)
5730	t->nodes[i].left = t->nodes[i].right = t->nodes[i].anc = NULL;
5731
5732	/* sort the tips in the t->allDownPass array */
5733	p = t->nodes;
5734	for (i=0; i<t->nNodes; i++, p++)
5735	t->allDownPass[p->index] = p;
5736
5737	/* make sure that root has index 0 */
5738	q = t->allDownPass[0];
5739	q->anc = q->right = NULL;
5740	t->root = q;
5741
5742	/* connect the first three tips */
5743	p = t->allDownPass[numTips];
5744	p->anc = q;
5745	q->left = p;
5746	p->length = *(brlens++);
5747	q = t->allDownPass[1];
5748	r = t->allDownPass[2];
5749	p->left = q;
5750	p->right = r;
5751	q->anc = r->anc = p;
5752	q->length = *(brlens++);
5753	r->length = *(brlens++);
5754
5755	/* add one tip at a time */
5756	for (i=3; i<numTips; i++)
5757	{
5758	p = t->allDownPass[i];
5759	q = t->allDownPass[numTips-2+i];
5760	r = t->allDownPass[order[i-3]];
5761	p->anc = q;
5762	q->left = p;
5763	q->right = r;
5764	q->anc = r->anc;
5765	if (r->anc->left == r)
5766	r->anc->left = q;
5767	else
5768	r->anc->right = q;
5769	r->anc = q;
5770	q->length = *(brlens++);
5771	p->length = *(brlens++);
5772	}
5773
5774	/* get downpass */
5775	GetDownPass (t);
5776
5777	/* relabel interior nodes (root is correctly labeled already) */
5778	for (i=0; i<t->nIntNodes; i++)
5779	t->intDownPass[i]->index = i+numTips;
5780
5781	return (NO_ERROR);
5782	}
5783
5784
5785
5786
5787
5788	void SetDatedNodeAges (Param *param, int chain, int state)
5789
5790	{
5791
5792	int i;
5793	MrBFlt clockRate;
5794	ModelInfo *m;
5795	TreeNode *p;
5796	Tree *t;
5797
5798	extern void ShowNodes(TreeNode *,int,int);
5799
5800	t = GetTree (param, chain, state);
5801	m = &modelSettings[t->relParts[0]];
5802
5803	if (m->clockRate == NULL)
5804	clockRate = 1.0;
5805	else
5806	clockRate = *GetParamVals(m->clockRate, chain, state);
5807
5808	for (i=0; i<t->nNodes-1; i++)
5809	{
5810	p = t->allDownPass[i];
5811	if (p->isDated == YES)
5812	p->age = p->nodeDepth / clockRate;
5813	else
5814	p->age = -1.0;
5815	}
5816	}
5817
5818
5819
5820
5821
5822	void SetNodeDepths (Tree *t)
5823
5824	{
5825
5826	int i;
5827	MrBFlt d1, d2;
5828	TreeNode *p;
5829
5830	extern void ShowNodes(TreeNode *,int,int);
5831
5832	for (i=0; i<t->nNodes-1; i++)
5833	{
5834	p = t->allDownPass[i];
5835	if (p->left == NULL)
5836	p->nodeDepth = 0.0;
5837	else
5838	{
5839	d1 = p->left->nodeDepth + p->left->length;
5840	d2 = p->right->nodeDepth + p->right->length;
5841	//assert (!(t->isCalibrated == NO && AreDoublesEqual(d1,d2,0.00001)==NO)); // may not work if we set startval topology of strict clock tree by non clock tree.
5842	if (d1 > d2)
5843	p->nodeDepth = d1;
5844	else
5845	p->nodeDepth = d2;
5846	}
5847	}
5848
5849	for (i=t->nNodes-3; i>=0; i--)
5850	{
5851	p = t->allDownPass[i];
5852	if (p->left == NULL && p->calibration == NULL)
5853	p->nodeDepth = 0.0;
5854	else
5855	p->nodeDepth = p->anc->nodeDepth - p->length;
5856	}
5857	}
5858
5859
5860
5861
5862
5863	/* Set ages of a clock tree according to depth and clockrate. Check that resulting ages are consistant with calibration.
5864	\| return YES if tree is age consistent, No otherwise.
5865	*/
5866	int SetTreeNodeAges (Param *param, int chain, int state)
5867	{
5868	Tree *tree;
5869	TreeNode *p;
5870	int i;
5871	MrBFlt clockRate;
5872
5873	if (param->paramType != P_TOPOLOGY && param->paramType != P_BRLENS && param->paramType != P_SPECIESTREE)
5874	return YES;
5875
5876	tree = GetTree(param, chain, state);
5877	if (modelSettings[param->relParts[0]].clockRate != NULL)
5878	clockRate = *GetParamVals(modelSettings[param->relParts[0]].clockRate, chain, state);
5879	else
5880	return YES;
5881
5882
5883	/* Clock trees */
5884
5885	/* Check that lengths and depths are consistant. That would work for the case when we set up branch lenght from starting tree */
5886	for (i=0; i<tree->nNodes-1; i++) {
5887	p = tree->allDownPass[i];
5888	p->age = p->nodeDepth / clockRate;
5889	}
5890
5891	/* Check that ages and calibrations are consistent */
5892	if (tree->isCalibrated == YES)
5893	{
5894	for (i=0; i<tree->nNodes-1; i++)
5895	{
5896	p = tree->allDownPass[i];
5897	if (p->isDated == YES) {
5898	if (p->calibration->prior == fixed && fabs(p->age - p->calibration->age) > 0.000001)
5899	{
5900	printf ("Node %d has age %f but should be fixed to age %f\n",
5901	p->index, p->age, p->calibration->age);
5902	return NO;
5903	}
5904	else if (p->calibration->prior == offsetExponential && p->age < p->calibration->offset)
5905	{
5906	printf ("Node %d has age %f but should be minimally of age %f\n",
5907	p->index, p->age, p->calibration->offset);
5908	return NO;
5909	}
5910	else if (p->calibration->prior == uniform && (p->age < p->calibration->min \|\| p->age > p->calibration->max))
5911	{
5912	printf ("Node %d has age %f but should be in the interval (%f,%f)\n",
5913	p->index, p->age, p->calibration->min, p->calibration->max);
5914	return NO;
5915	}
5916	}
5917	}
5918	}
5919
5920
5921	return YES;
5922	}
5923
5924
5925
5926
5927
5928
5929	int ShowPolyNodes (PolyTree *pt)
5930
5931	{
5932
5933	int i;
5934	PolyNode *p;
5935
5936	/* this is the tree, on a node-by-node basis */
5937	printf (" memnodes = %d nNodes = %d nIntNodes = %d root = %d\n", pt->memNodes, pt->nNodes, pt->nIntNodes, pt->root->index);
5938	printf (" isRooted = %d\n", pt->isRooted);
5939	printf (" no. index (left sib anc) -- locked/free -- label (p->x)\n");
5940	for (i=0; i<pt->memNodes; i++)
5941	{
5942	p = &pt->nodes[i];
5943	if (!(p->left == NULL && p->sib == NULL && p->anc == NULL))
5944	{
5945	printf ("%4d -- %4d ", i, p->index);
5946	if (p->left != NULL)
5947	printf ("(%4d ", p->left->index);
5948	else
5949	printf ("(null ");
5950
5951	if (p->sib != NULL)
5952	printf ("%4d ", p->sib->index);
5953	else
5954	printf ("null ");
5955
5956	if (p->anc != NULL)
5957	printf ("%4d)", p->anc->index);
5958	else
5959	printf ("null)");
5960
5961	if (p->isLocked == YES)
5962	printf ("-- locked -- ");
5963	else
5964	printf ("-- free --");
5965
5966	if (p->left == NULL && p->anc != NULL)
5967	printf (" \"%s\" (%d)\n", p->label, p->x);
5968	else
5969	printf (" \"\" (%d)\n", p->x);
5970	}
5971	}
5972
5973	return NO_ERROR;
5974	}
5975
5976
5977
5978
5979
5980	/* ShowTree: Show tree on screen */
5981	int ShowTree (Tree *t)
5982
5983	{
5984
5985	int i, j, k, x, nLines, nLevels, levelDepth, from, to;
5986	char treeLine[SCREENWIDTH2], labelLine[100];
5987	TreeNode *p;
5988
5989	/* get coordinates */
5990	x = 0;
5991	nLines = 0;
5992	for (i=0; i<t->nNodes; i++)
5993	{
5994	p = t->allDownPass[i];
5995	if (p->left == NULL && p->right == NULL)
5996	{
5997	p->x = x;
5998	x += 2;
5999	p->y = 0;
6000	nLines += 2;
6001	}
6002	else if (p->left != NULL && p->right != NULL && p->anc != NULL)
6003	{
6004	p->x = p->left->x + (p->right->x - p->left->x) / 2;
6005	if (p->left->y > p->right->y)
6006	p->y = p->left->y + 1;
6007	else
6008	p->y = p->right->y + 1;
6009	}
6010	else
6011	{
6012	p->x = x;
6013	x += 2;
6014	p->y = 0;
6015	}
6016	}
6017
6018	/* print tree out, line-by-line */
6019	levelDepth = SCREENWIDTH / t->root->left->y;
6020	nLevels = t->root->left->y;
6021	for (j=0; j<=nLines-2; j++)
6022	{
6023	for (i=0; i<SCREENWIDTH2-2; i++)
6024	treeLine[i] = ' ';
6025	treeLine[SCREENWIDTH-1] = '\n';
6026	if (j % 2 == 0)
6027	{
6028	for (i=0; i<t->nNodes; i++)
6029	{
6030	p = t->allDownPass[i];
6031	if (p->left == NULL && p->x == j)
6032	{
6033	strcpy (labelLine, p->label);
6034	}
6035	}
6036	}
6037	for (i=0; i<t->nNodes; i++)
6038	{
6039	p = t->allDownPass[i];
6040	if (p->anc != NULL)
6041	{
6042	if (p->anc->anc != NULL)
6043	{
6044	if (p->x == j)
6045	{
6046	from = (nLevels - p->anc->y) * levelDepth;
6047	to = (nLevels - p->y) * levelDepth;
6048	if (p->y == 0)
6049	to = SCREENWIDTH-1;
6050	if (to >= SCREENWIDTH)
6051	to = SCREENWIDTH-1;
6052
6053	for (k=from; k<to; k++)
6054	treeLine[k] = '-';
6055	if (p->anc->left == p)
6056	treeLine[from] = '/';
6057	else
6058	treeLine[from] = '\\';
6059	if (p->left != NULL)
6060	{
6061	treeLine[to] = '+';
6062	}
6063	if (p->anc->anc == t->root && p->anc->right == p)
6064	{
6065	if (t->isRooted == NO)
6066	treeLine[to] = '+';
6067	else
6068	treeLine[from] = '\\';
6069	}
6070	}
6071	else
6072	{
6073	if (p->left != NULL && p->right != NULL)
6074	{
6075	if (j < p->x && j > p->left->x)
6076	{
6077	from = (nLevels - p->y) * levelDepth;
6078	treeLine[from] = '\|';
6079	}
6080	else if (j > p->x && j < p->right->x && p->left != NULL)
6081	{
6082	from = (nLevels - p->y) * levelDepth;
6083	treeLine[from] = '\|';
6084	}
6085	}
6086	}
6087	}
6088	else
6089	{
6090	if (p->x == j)
6091	{
6092	treeLine[0] = '\|'; /* temp */
6093	}
6094	else if (j < p->x && j > p->left->x)
6095	{
6096	treeLine[0] = '\|';
6097	}
6098	else if (j > p->x && j < p->right->x)
6099	{
6100	treeLine[0] = '\|';
6101	}
6102	if (t->isRooted == NO)
6103	{
6104	if (j > p->x && j <= nLines-2)
6105	treeLine[0] = '\|';
6106	if (j == p->right->x)
6107	treeLine[0] = '+';
6108	}
6109	else
6110	{
6111	if (j == p->x)
6112	treeLine[0] = '+';
6113	}
6114	}
6115	}
6116	}
6117	treeLine[SCREENWIDTH-1] = '\0';
6118	if (j % 2 == 0)
6119	MrBayesPrint (" %s %s\n", treeLine, labelLine);
6120	else
6121	MrBayesPrint (" %s \n", treeLine);
6122	}
6123
6124	if (t->isRooted == NO)
6125	{
6126	for (i=0; i<SCREENWIDTH; i++)
6127	treeLine[i] = ' ';
6128	treeLine[SCREENWIDTH-1] = '\0';
6129	MrBayesPrint (" \|\n");
6130	for (k=0; k<SCREENWIDTH; k++)
6131	treeLine[k] = '-';
6132	treeLine[SCREENWIDTH-1] = '\0';
6133	treeLine[0] = '\\';
6134	strcpy (labelLine, t->root->label);
6135	labelLine[19] = '\0';
6136	MrBayesPrint (" %s %s\n", treeLine, labelLine);
6137	}
6138
6139	#if defined (DEBUG_CONSTRAINTS)
6140	for (i=0; i<t->nNodes; i++)
6141	printf ("%d -- %s\n", t->allDownPass[i]->index + 1, t->allDownPass[i]->isLocked == YES ? "locked" : "free");
6142	#endif
6143
6144	return (NO_ERROR);
6145
6146	}
6147
6148
6149
6150
6151
6152	/*-------------------------------------------------------
6153	\|
6154	\| StoreRPolyTopology: Same as StoreRTopology but for
6155	\| binary polytree source trees.
6156	\|
6157	--------------------------------------------------------*/
6158	int StoreRPolyTopology (PolyTree t, int order)
6159
6160	{
6161	int i, numTaxa;
6162	PolyNode p, q;
6163
6164	/* find number of taxa */
6165	numTaxa = t->nNodes - t->nIntNodes;
6166
6167	/* first get the terminal taxon positions and store
6168	them in the order array. */
6169	for (i=0; i<t->nNodes; i++)
6170	{
6171	p = t->allDownPass[i];
6172	/* we do not need to worry about the first two taxa */
6173	if (p->index > 1 && p->index < numTaxa)
6174	order[p->index-2] = i;
6175	}
6176
6177	/* label the interior nodes with the correct index */
6178	for (i=0; i<t->nNodes; i++)
6179	{
6180	p = t->allDownPass[i];
6181	if (p->left == NULL)
6182	p->x = p->y = p->index;
6183	else
6184	{
6185	if (p->left->y < p->left->sib->y)
6186	{
6187	p->y = p->left->y;
6188	p->x = p->left->sib->y + numTaxa - 1;
6189	}
6190	else
6191	{
6192	p->y = p->left->sib->y;
6193	p->x = p->left->y + numTaxa - 1;
6194	}
6195	}
6196	}
6197
6198	/* break the tree into pieces */
6199	for (i=0; i<numTaxa-2; i++)
6200	{
6201	/* find the next node to remove */
6202	p = t->allDownPass[order[numTaxa-3-i]];
6203	q = p->anc;
6204	if (q->left == p)
6205	{
6206	order[numTaxa-3-i] = q->left->sib->x;
6207	p->sib->anc = q->anc;
6208	if (q->anc == NULL)
6209	{
6210	p->sib->left->sib->sib = p->sib->sib;
6211	p->sib->sib = NULL;
6212	}
6213	else if (q->anc->left == q)
6214	{
6215	q->anc->left = q->left->sib;
6216	p->sib->sib = q->sib;
6217	}
6218	else
6219	q->anc->left->sib = q->left->sib;
6220	}
6221	else
6222	{
6223	order[numTaxa-3-i] = q->left->x;
6224	q->left->anc = q->anc;
6225	if (q->anc == NULL)
6226	{
6227	q->left->left->sib->sib = p->sib;
6228	q->left->sib = NULL;
6229	}
6230	else if (q->anc->left == q)
6231	{
6232	q->anc->left = q->left;
6233	q->anc->left->sib = q->sib;
6234	}
6235	else
6236	{
6237	q->anc->left->sib = q->left;
6238	q->left->sib = NULL;
6239	}
6240	}
6241	}
6242
6243	return (NO_ERROR);
6244	}
6245
6246
6247
6248
6249
6250	/*-------------------------------------------------------
6251	\|
6252	\| StoreRPolyTree: Same as StoreRTree but for
6253	\| binary rooted polytree source trees.
6254	\|
6255	--------------------------------------------------------*/
6256	int StoreRPolyTree (PolyTree t, int order, MrBFlt *brlens)
6257
6258	{
6259	int i, j, numTaxa;
6260	PolyNode p, q;
6261
6262	/* find number of taxa */
6263	numTaxa = t->nNodes - t->nIntNodes;
6264
6265	/* first get the terminal taxon positions and store
6266	them in the order array. */
6267	for (i=0; i<t->nNodes; i++)
6268	{
6269	p = t->allDownPass[i];
6270	/* we do not need to worry about the first two taxa */
6271	if (p->index > 1 && p->index < numTaxa)
6272	order[p->index-2] = i;
6273	}
6274
6275	/* label the interior nodes with the correct index */
6276	for (i=0; i<t->nNodes; i++)
6277	{
6278	p = t->allDownPass[i];
6279	if (p->left == NULL)
6280	p->x = p->y = p->index;
6281	else
6282	{
6283	if (p->left->y < p->left->sib->y)
6284	{
6285	p->y = p->left->y;
6286	p->x = p->left->sib->y + numTaxa - 1;
6287	}
6288	else
6289	{
6290	p->y = p->left->sib->y;
6291	p->x = p->left->y + numTaxa - 1;
6292	}
6293	}
6294	}
6295
6296	/* break the tree into pieces */
6297	j = t->nNodes - 2; /* index of first branch length */
6298	for (i=0; i<numTaxa-2; i++)
6299	{
6300	/* find the next node to remove */
6301	p = t->allDownPass[order[numTaxa-3-i]];
6302	q = p->anc;
6303	brlens[j--] = p->length;
6304	brlens[j--] = q->length;
6305	if (q->left == p)
6306	{
6307	order[numTaxa-3-i] = q->left->sib->x;
6308	p->sib->anc = q->anc;
6309	if (q->anc == NULL)
6310	{
6311	p->sib->left->sib->sib = p->sib->sib;
6312	p->sib->sib = NULL;
6313	}
6314	else if (q->anc->left == q)
6315	{
6316	q->anc->left = q->left->sib;
6317	p->sib->sib = q->sib;
6318	}
6319	else
6320	q->anc->left->sib = q->left->sib;
6321	}
6322	else
6323	{
6324	order[numTaxa-3-i] = q->left->x;
6325	q->left->anc = q->anc;
6326	if (q->anc == NULL)
6327	{
6328	q->left->left->sib->sib = p->sib;
6329	q->left->sib = NULL;
6330	}
6331	else if (q->anc->left == q)
6332	{
6333	q->anc->left = q->left;
6334	q->anc->left->sib = q->sib;
6335	}
6336	else
6337	{
6338	q->anc->left->sib = q->left;
6339	q->left->sib = NULL;
6340	}
6341	}
6342	}
6343
6344	/* store the last two lengths; index 0 and 1 */
6345	p = t->root;
6346	brlens[p->left->index] = p->left->length;
6347	brlens[p->left->sib->index] = p->left->sib->length;
6348
6349	return (NO_ERROR);
6350	}
6351
6352
6353
6354
6355
6356	/*-------------------------------------------------------
6357	\|
6358	\| StoreRTopology: This routine will break a rooted tree
6359	\| into an array of ints describing the structure
6360	\| of the tree. The tree will be destroyed
6361	\| in the process (the node pointers, that is).
6362	\| However, the tree is not deleted.
6363	\|
6364	--------------------------------------------------------*/
6365	int StoreRTopology (Tree t, int order)
6366
6367	{
6368	int i, numTaxa;
6369	TreeNode p, q;
6370
6371	/* find number of taxa */
6372	numTaxa = t->nNodes - t->nIntNodes - 1;
6373
6374	/* first get the terminal taxon positions and store
6375	them in the order array. */
6376	for (i=0; i<t->nNodes; i++)
6377	{
6378	p = t->allDownPass[i];
6379	/* we do not need to worry about the first two taxa */
6380	if (p->index > 1 && p->index < numTaxa)
6381	order[p->index-2] = i;
6382	}
6383
6384	/* label the interior nodes with the correct index */
6385	for (i=0; i<t->nNodes; i++)
6386	{
6387	p = t->allDownPass[i];
6388	if (p->left == NULL)
6389	p->x = p->y = p->index;
6390	else if (p->right != NULL)
6391	{
6392	if (p->left->y < p->right->y)
6393	{
6394	p->y = p->left->y;
6395	p->x = p->right->y + numTaxa - 1;
6396	}
6397	else
6398	{
6399	p->y = p->right->y;
6400	p->x = p->left->y + numTaxa - 1;
6401	}
6402	}
6403	}
6404
6405	/* break the tree into pieces */
6406	for (i=0; i<numTaxa-2; i++)
6407	{
6408	/* find the next node to remove */
6409	p = t->allDownPass[order[numTaxa-3-i]];
6410	q = p->anc;
6411	if (q->left == p)
6412	{
6413	order[numTaxa-3-i] = q->right->x;
6414	q->right->anc = q->anc;
6415	if (q->anc->left == q)
6416	q->anc->left = q->right;
6417	else
6418	q->anc->right = q->right;
6419	}
6420	else
6421	{
6422	order[numTaxa-3-i] = q->left->x;
6423	q->left->anc = q->anc;
6424	if (q->anc->left == q)
6425	q->anc->left = q->left;
6426	else
6427	q->anc->right = q->left;
6428	}
6429	}
6430
6431	return (NO_ERROR);
6432	}
6433
6434
6435
6436
6437
6438	/*-------------------------------------------------------
6439	\|
6440	\| StoreRTree: This routine will break a rooted tree
6441	\| into an array of ints describing the structure
6442	\| of the tree and an array of doubles storing
6443	\| the branch lengths. The tree will be
6444	\| destroyed in the process (the node pointers,
6445	\| that is). However, the tree is not deleted.
6446	\|
6447	--------------------------------------------------------*/
6448	int StoreRTree (Tree t, int order, MrBFlt *brlens)
6449
6450	{
6451	int i, j, numTaxa;
6452	TreeNode p, q;
6453
6454	extern void ShowNodes (TreeNode *p, int indent, int isRooted);
6455
6456	/* find number of taxa */
6457	numTaxa = t->nNodes - t->nIntNodes - 1;
6458
6459	/* first get the terminal taxon positions and store
6460	them in the order array. */
6461	for (i=0; i<t->nNodes; i++)
6462	{
6463	p = t->allDownPass[i];
6464	/* we do not need to worry about the first two taxa */
6465	if (p->index > 1 && p->index < numTaxa)
6466	order[p->index-2] = i;
6467	}
6468
6469	/* label the interior nodes with the correct index */
6470	for (i=0; i<t->nNodes; i++)
6471	{
6472	p = t->allDownPass[i];
6473	if (p->left == NULL)
6474	p->x = p->y = p->index;
6475	else if (p->right != NULL)
6476	{
6477	if (p->left->y < p->right->y)
6478	{
6479	p->y = p->left->y;
6480	p->x = p->right->y + numTaxa - 1;
6481	}
6482	else
6483	{
6484	p->y = p->right->y;
6485	p->x = p->left->y + numTaxa - 1;
6486	}
6487	}
6488	}
6489
6490	/* break the tree into pieces */
6491	j = 2 * numTaxa - 3;
6492	for (i=0; i<numTaxa-2; i++)
6493	{
6494	/* find the next node to remove */
6495	p = t->allDownPass[order[numTaxa-3-i]];
6496	q = p->anc;
6497	brlens[j--] = p->length;
6498	if (q->left == p)
6499	{
6500	if (q->anc->anc != NULL)
6501	brlens[j--] = q->length;
6502	else
6503	brlens[j--] = q->right->length;
6504	order[numTaxa-3-i] = q->right->x;
6505	q->right->anc = q->anc;
6506	if (q->anc->left == q)
6507	q->anc->left = q->right;
6508	else
6509	q->anc->right = q->right;
6510	}
6511	else
6512	{
6513	if (q->anc->anc != NULL)
6514	brlens[j--] = q->length;
6515	else
6516	brlens[j--] = q->left->length;
6517	order[numTaxa-3-i] = q->left->x;
6518	q->left->anc = q->anc;
6519	if (q->anc->left == q)
6520	q->anc->left = q->left;
6521	else
6522	q->anc->right = q->left;
6523	}
6524	}
6525
6526	/* store the final two branch lengths in the right order; they have indices 0 and 1 */
6527	p = t->root->left;
6528	brlens[p->left->index] = p->left->length;
6529	brlens[p->right->index] = p->right->length;
6530
6531	return (NO_ERROR);
6532	}
6533
6534
6535
6536
6537
6538	/*-------------------------------------------------------
6539	\|
6540	\| StoreRTreeWithIndices: This routine will break a rooted
6541	\| tree into an array of ints describing the structure
6542	\| of the tree and the interior node indices, and an array
6543	\| of doubles storing the branch lengths. The tree will be
6544	\| destroyed in the process (the node pointers,
6545	\| that is). However, the tree is not deleted.
6546	\|
6547	--------------------------------------------------------*/
6548	int StoreRTreeWithIndices (Tree t, int order, MrBFlt *brlens)
6549
6550	{
6551	int i, j, k, numTaxa;
6552	TreeNode p, q;
6553
6554	extern void ShowNodes (TreeNode *p, int indent, int isRooted);
6555
6556	/* find number of taxa */
6557	numTaxa = t->nNodes - t->nIntNodes - 1;
6558
6559	/* first get the terminal taxon positions and store
6560	them in the order array. */
6561	for (i=0; i<t->nNodes; i++)
6562	{
6563	p = t->allDownPass[i];
6564	/* we do not need to worry about the first two taxa */
6565	if (p->index > 1 && p->index < numTaxa)
6566	order[p->index-2] = i;
6567	}
6568
6569	/* label the interior nodes with the correct index */
6570	for (i=0; i<t->nNodes; i++)
6571	{
6572	p = t->allDownPass[i];
6573	if (p->left == NULL)
6574	p->x = p->y = p->index;
6575	else if (p->right != NULL)
6576	{
6577	if (p->left->y < p->right->y)
6578	{
6579	p->y = p->left->y;
6580	p->x = p->right->y + numTaxa - 1;
6581	}
6582	else
6583	{
6584	p->y = p->right->y;
6585	p->x = p->left->y + numTaxa - 1;
6586	}
6587	}
6588	}
6589
6590	/* break the tree into pieces */
6591	j = 2 * numTaxa - 3;
6592	k = 2*(numTaxa - 2);
6593	for (i=0; i<numTaxa-2; i++)
6594	{
6595	/* find the next node to remove */
6596	p = t->allDownPass[order[numTaxa-3-i]];
6597	q = p->anc;
6598	brlens[j--] = p->length;
6599	if (q->left == p)
6600	{
6601	if (q->anc->anc != NULL)
6602	brlens[j--] = q->length;
6603	else
6604	brlens[j--] = q->right->length;
6605	order[k--] = q->right->x;
6606	order[k--] = q->index;
6607	q->right->anc = q->anc;
6608	if (q->anc->left == q)
6609	q->anc->left = q->right;
6610	else
6611	q->anc->right = q->right;
6612	}
6613	else
6614	{
6615	if (q->anc->anc != NULL)
6616	brlens[j--] = q->length;
6617	else
6618	brlens[j--] = q->left->length;
6619	order[k--] = q->left->x;
6620	order[k--] = q->index;
6621	q->left->anc = q->anc;
6622	if (q->anc->left == q)
6623	q->anc->left = q->left;
6624	else
6625	q->anc->right = q->left;
6626	}
6627	}
6628
6629	/* store the final two branch lengths in the right order; they have indices 0 and 1 */
6630	p = t->root->left;
6631	order[k] = p->index;
6632	brlens[p->left->index] = p->left->length;
6633	brlens[p->right->index] = p->right->length;
6634
6635	return (NO_ERROR);
6636	}
6637
6638
6639
6640
6641
6642	/*-------------------------------------------------------
6643	\|
6644	\| StoreUPolyTopology: Same as StoreUTopology but for
6645	\| binary polytree source.
6646	\|
6647	--------------------------------------------------------*/
6648	int StoreUPolyTopology (PolyTree t, int order)
6649
6650	{
6651	int i, numTips;
6652	PolyNode p, q;
6653
6654	/* check if the tree is rooted on taxon 0 */
6655	if (t->root->left->sib->sib->index != 0)
6656	MovePolyCalculationRoot (t, 0);
6657
6658	/* rearrange the root */
6659	t->root->anc = t->root->left->sib->sib;
6660	t->root->left->sib->sib = NULL;
6661	t->root->anc->left = t->root;
6662	t->root->anc->sib = NULL;
6663	t->root->anc->anc = NULL;
6664	t->root = t->root->anc;
6665
6666	/* find number of tips */
6667	numTips = t->nNodes - t->nIntNodes;
6668
6669	/* first get the terminal taxon positions and store
6670	them in the order array. */
6671	for (i=0; i<t->nNodes; i++)
6672	{
6673	p = t->allDownPass[i];
6674	/* we do not need to worry about the first three taxa */
6675	if (p->index > 2 && p->index < numTips)
6676	order[p->index-3] = i;
6677	}
6678
6679	/* label the interior nodes with the correct index */
6680	for (i=0; i<t->nNodes; i++)
6681	{
6682	p = t->allDownPass[i];
6683	if (p->left == NULL \|\| p->anc == NULL)
6684	p->x = p->y = p->index;
6685	else
6686	{
6687	if (p->left->y < p->left->sib->y)
6688	{
6689	p->y = p->left->y;
6690	p->x = p->left->sib->y + numTips - 2;
6691	}
6692	else
6693	{
6694	p->y = p->left->sib->y;
6695	p->x = p->left->y + numTips - 2;
6696	}
6697	}
6698	}
6699
6700	/* break the tree into pieces */
6701	for (i=0; i<numTips-3; i++)
6702	{
6703	/* find the next node to remove */
6704	p = t->allDownPass[order[numTips-4-i]];
6705	q = p->anc;
6706	if (q->left == p)
6707	{
6708	order[numTips-4-i] = q->left->sib->x;
6709	p->sib->anc = q->anc;
6710	if (q->anc->left == q)
6711	{
6712	q->anc->left = p->sib;
6713	p->sib->sib = q->sib;
6714	}
6715	else
6716	{
6717	q->anc->left->sib = p->sib;
6718	p->sib->sib = q->sib;
6719	}
6720	}
6721	else
6722	{
6723	order[numTips-4-i] = q->left->x;
6724	q->left->anc = q->anc;
6725	if (q->anc->left == q)
6726	{
6727	q->anc->left = q->left;
6728	q->left->sib = q->sib;
6729	}
6730	else
6731	{
6732	q->anc->left->sib = q->left;
6733	q->left->sib = q->sib;
6734	}
6735	}
6736	}
6737
6738	return (NO_ERROR);
6739	}
6740
6741
6742
6743
6744
6745	/*-------------------------------------------------------
6746	\|
6747	\| StoreUPolyTree: Same as StoreUTopology but for
6748	\| binary polytree source.
6749	\|
6750	--------------------------------------------------------*/
6751	int StoreUPolyTree (PolyTree t, int order, MrBFlt *brlens)
6752
6753	{
6754	int i, j, numTips;
6755	PolyNode p, q;
6756
6757	/* check if the tree is rooted on taxon 0 */
6758	if (t->root->left->sib->sib->index != 0)
6759	MovePolyCalculationRoot (t, 0);
6760
6761	/* rearrange the root */
6762	t->root->anc = t->root->left->sib->sib;
6763	t->root->left->sib->sib = NULL;
6764	t->root->anc->left = t->root;
6765	t->root->anc->sib = NULL;
6766	t->root->anc->anc = NULL;
6767	t->root = t->root->anc;
6768
6769	/* find number of tips */
6770	numTips = t->nNodes - t->nIntNodes;
6771
6772	/* first get the terminal taxon positions and store
6773	them in the order array. */
6774	for (i=0; i<t->nNodes; i++)
6775	{
6776	p = t->allDownPass[i];
6777	/* we do not need to worry about the first three taxa */
6778	if (p->index > 2 && p->index < numTips)
6779	order[p->index-3] = i;
6780	}
6781
6782	/* label the interior nodes with the correct index */
6783	for (i=0; i<t->nNodes; i++)
6784	{
6785	p = t->allDownPass[i];
6786	if (p->left == NULL \|\| p->anc == NULL)
6787	p->x = p->y = p->index;
6788	else
6789	{
6790	if (p->left->y < p->left->sib->y)
6791	{
6792	p->y = p->left->y;
6793	p->x = p->left->sib->y + numTips - 2;
6794	}
6795	else
6796	{
6797	p->y = p->left->sib->y;
6798	p->x = p->left->y + numTips - 2;
6799	}
6800	}
6801	}
6802
6803	/* break the tree into pieces */
6804	j = 2*numTips - 4;
6805	for (i=0; i<numTips-3; i++)
6806	{
6807	/* find the next node to remove */
6808	p = t->allDownPass[order[numTips-4-i]];
6809	assert (p->index > 2 && p->index < numTips);
6810	assert (p->anc->anc != NULL);
6811	q = p->anc;
6812	brlens[j--] = p->length;
6813	brlens[j--] = q->length;
6814	if (q->left == p)
6815	{
6816	order[numTips-4-i] = q->left->sib->x;
6817	p->sib->anc = q->anc;
6818	if (q->anc->left == q)
6819	{
6820	q->anc->left = p->sib;
6821	p->sib->sib = q->sib;
6822	}
6823	else
6824	{
6825	q->anc->left->sib = p->sib;
6826	p->sib->sib = q->sib;
6827	}
6828	}
6829	else
6830	{
6831	order[numTips-4-i] = q->left->x;
6832	q->left->anc = q->anc;
6833	if (q->anc->left == q)
6834	{
6835	q->anc->left = q->left;
6836	q->left->sib = q->sib;
6837	}
6838	else
6839	{
6840	q->anc->left->sib = q->left;
6841	q->left->sib = q->sib;
6842	}
6843	}
6844	}
6845
6846	/* store last three branch lengths, index 0, 1, 2 */
6847	q = t->root;
6848	assert (q->index == 0);
6849	brlens[q->index] = q->length;
6850	q = q->left->left;
6851	assert (q->index == 1 \|\| q->index == 2);
6852	brlens[q->index] = q->length;
6853	q = q->sib;
6854	assert (q->index == 1 \|\| q->index == 2);
6855	brlens[q->index] = q->length;
6856
6857	return (NO_ERROR);
6858	}
6859
6860
6861
6862
6863
6864	/*-------------------------------------------------------
6865	\|
6866	\| StoreUTopology: This routine will break an unrooted tree
6867	\| into an array of ints describing the structure
6868	\| of the tree. The tree will be destroyed
6869	\| in the process (the node pointers, that is).
6870	\| However, the tree is not deleted.
6871	\|
6872	--------------------------------------------------------*/
6873	int StoreUTopology (Tree t, int order)
6874
6875	{
6876	int i, numTips;
6877	TreeNode p, q;
6878
6879	/* check if the tree is rooted on taxon 0 */
6880	if (t->root->index != 0)
6881	MoveCalculationRoot (t, 0);
6882
6883	/* find number of tips */
6884	numTips = t->nNodes - t->nIntNodes;
6885
6886	/* first get the terminal taxon positions and store
6887	them in the order array. */
6888	for (i=0; i<t->nNodes; i++)
6889	{
6890	p = t->allDownPass[i];
6891	/* we do not need to worry about the first three taxa */
6892	if (p->index > 2 && p->index < numTips)
6893	order[p->index-3] = i;
6894	}
6895
6896	/* label the interior nodes with the correct index */
6897	for (i=0; i<t->nNodes; i++)
6898	{
6899	p = t->allDownPass[i];
6900	if (p->left == NULL)
6901	p->x = p->y = p->index;
6902	else if (p->right != NULL)
6903	{
6904	if (p->left->y < p->right->y)
6905	{
6906	p->y = p->left->y;
6907	p->x = p->right->y + numTips - 2;
6908	}
6909	else
6910	{
6911	p->y = p->right->y;
6912	p->x = p->left->y + numTips - 2;
6913	}
6914	}
6915	}
6916
6917	/* break the tree into pieces */
6918	for (i=0; i<numTips-3; i++)
6919	{
6920	/* find the next node to remove */
6921	p = t->allDownPass[order[numTips-4-i]];
6922	q = p->anc;
6923	if (q->left == p)
6924	{
6925	order[numTips-4-i] = q->right->x;
6926	q->right->anc = q->anc;
6927	if (q->anc->left == q)
6928	q->anc->left = q->right;
6929	else
6930	q->anc->right = q->right;
6931	}
6932	else
6933	{
6934	order[numTips-4-i] = q->left->x;
6935	q->left->anc = q->anc;
6936	if (q->anc->left == q)
6937	q->anc->left = q->left;
6938	else
6939	q->anc->right = q->left;
6940	}
6941	}
6942
6943	return (NO_ERROR);
6944	}
6945
6946
6947
6948
6949
6950	/*-------------------------------------------------------
6951	\|
6952	\| StoreUTree: This routine will break an unrooted tree
6953	\| into an array of ints describing the structure
6954	\| of the tree and an array of doubles storing
6955	\| the branch lengths. The tree will be
6956	\| destroyed in the process (the node pointers,
6957	\| that is). However, the tree is not deleted.
6958	\|
6959	--------------------------------------------------------*/
6960	int StoreUTree (Tree t, int order, MrBFlt *brlens)
6961
6962	{
6963	int i, j, numTips;
6964	TreeNode p, q;
6965
6966	/* check if the tree is rooted on taxon 0 */
6967	if (t->root->index != 0)
6968	MoveCalculationRoot(t, 0);
6969
6970	/* find number of tips */
6971	numTips = t->nNodes - t->nIntNodes;
6972
6973	/* first get the terminal taxon positions and store
6974	them in the order array. */
6975	for (i=0; i<t->nNodes; i++)
6976	{
6977	p = t->allDownPass[i];
6978	/* we do not need to worry about the first three taxa */
6979	if (p->index > 2 && p->index < numTips)
6980	order[p->index-3] = i;
6981	}
6982
6983	/* label the interior nodes with the correct index */
6984	for (i=0; i<t->nNodes; i++)
6985	{
6986	p = t->allDownPass[i];
6987	if (p->left == NULL)
6988	p->x = p->y = p->index;
6989	else if (p->right != NULL)
6990	{
6991	if (p->left->y < p->right->y)
6992	{
6993	p->y = p->left->y;
6994	p->x = p->right->y + numTips - 2;
6995	}
6996	else
6997	{
6998	p->y = p->right->y;
6999	p->x = p->left->y + numTips - 2;
7000	}
7001	}
7002	}
7003
7004	/* break the tree into pieces */
7005	j = 2 * numTips - 4;
7006	for (i=0; i<numTips-3; i++)
7007	{
7008	/* find the next node to remove */
7009	p = t->allDownPass[order[numTips-4-i]];
7010	q = p->anc;
7011	brlens[j--] = p->length;
7012	brlens[j--] = q->length;
7013	if (q->left == p)
7014	{
7015	order[numTips-4-i] = q->right->x;
7016	q->right->anc = q->anc;
7017	if (q->anc->left == q)
7018	q->anc->left = q->right;
7019	else
7020	q->anc->right = q->right;
7021	}
7022	else
7023	{
7024	order[numTips-4-i] = q->left->x;
7025	q->left->anc = q->anc;
7026	if (q->anc->left == q)
7027	q->anc->left = q->left;
7028	else
7029	q->anc->right = q->left;
7030	}
7031	}
7032
7033	/* store the final three branch lengths */
7034	/* we need to check the rotation of the tree to
7035	store the brlens in the right order (after node index) */
7036	p = t->root->left;
7037	if (p->right->index == 2)
7038	{
7039	brlens[j--] = p->right->length;
7040	brlens[j--] = p->left->length;
7041	}
7042	else
7043	{
7044	brlens[j--] = p->left->length;
7045	brlens[j--] = p->right->length;
7046	}
7047	brlens[j--] = p->length;
7048
7049	return (NO_ERROR);
7050	}
7051
7052
7053
7054
7055
7056	/* TreeLength: Calculate tree length */
7057	MrBFlt TreeLen (Tree *t)
7058	{
7059	int i, numLenNodes;
7060	MrBFlt len = 0.0;
7061
7062	if (t->isRooted == NO)
7063	numLenNodes = t->nNodes - 1;
7064	else
7065	numLenNodes = t->nNodes - 2;
7066
7067	for (i=0; i<numLenNodes; i++)
7068	len += t->allDownPass[i]->length;
7069
7070	return len;
7071	}
7072
7073
7074
7075
7076
7077	/*-------------------------------------------------------------------------------------------
7078	\|
7079	\| Unmark: This routine will unmark a subtree rooted at p
7080	\|
7081	---------------------------------------------------------------------------------------------*/
7082	void Unmark (TreeNode *p)
7083	{
7084	if (p != NULL)
7085	{
7086	p->marked = NO;
7087	Unmark (p->left);
7088	Unmark (p->right);
7089	}
7090	}
7091
7092
7093
7094
7095
7096	void WriteEventTree (TreeNode p, int chain, Param param)
7097
7098	{
7099	int j, nEvents;
7100	MrBFlt brlen, position, rateMult;
7101
7102	if (p != NULL)
7103	{
7104	if (p->left == NULL && p->right == NULL)
7105	{
7106	printf ("%d:%s", p->index + 1, MbPrintNum(p->length));
7107	if (param->paramType == P_CPPEVENTS)
7108	{
7109	nEvents = param->nEvents[2*chain+state[chain]][p->index];
7110	if (nEvents > 0)
7111	{
7112	printf ("[&E %s %d: (", param->name, nEvents);
7113	position = param->position[2*chain+state[chain]][p->index];
7114	rateMult = param->rateMult[2*chain+state[chain]][p->index];
7115	for (j=0; j<nEvents; j++)
7116	{
7117	printf ("%s", MbPrintNum(position[j]) );
7118	printf (" %s", MbPrintNum(rateMult[j]) );
7119	if (j != nEvents-1)
7120	printf (", ");
7121	}
7122	printf (")]");
7123	}
7124	else
7125	printf ("[&E %s 0]", param->name);
7126	}
7127	brlen = (GetParamSubVals (param, chain, state[chain])[p->index] + GetParamVals (param, chain, state[chain])[p->anc->index]) / 2.0;
7128	printf ("[&B %s %s]", param->name, MbPrintNum(brlen));
7129	}
7130	else
7131	{
7132	if (p->anc != NULL)
7133	printf ("(");
7134	WriteEventTree(p->left, chain, param);
7135	printf (",");
7136	WriteEventTree(p->right, chain, param);
7137	if (p->anc != NULL)
7138	{
7139	if (p->anc->anc != NULL)
7140	{
7141	printf ("):%s", MbPrintNum(p->length));
7142	if (param->paramType == P_CPPEVENTS)
7143	{
7144	nEvents = param->nEvents[2*chain+state[chain]][p->index];
7145	if (nEvents > 0)
7146	{
7147	printf ("[&E %s %d: (", param->name, nEvents);
7148	position = param->position[2*chain+state[chain]][p->index];
7149	rateMult = param->rateMult[2*chain+state[chain]][p->index];
7150	for (j=0; j<nEvents; j++)
7151	{
7152	printf ("%s", MbPrintNum(position[j]) );
7153	printf (" %s", MbPrintNum(rateMult[j]) );
7154	if (j != nEvents-1)
7155	printf (", ");
7156	}
7157	printf (")]");
7158	}
7159	else
7160	printf ("[&E %s 0]", param->name);
7161	}
7162	brlen = (GetParamSubVals (param, chain, state[chain])[p->index] + GetParamVals (param, chain, state[chain])[p->anc->index]) / 2.0;
7163	printf ("[&B %s %s]", param->name, MbPrintNum(brlen));
7164	}
7165	else
7166	printf(")");
7167	}
7168	}
7169	}
7170	}
7171
7172
7173
7174
7175
7176	void WriteEventTreeToPrintString (TreeNode p, int chain, Param param, int printAll)
7177
7178	{
7179	char *tempStr;
7180	int i, j, nEvents, tempStrSize = TEMPSTRSIZE;
7181	MrBFlt brlen, position, rateMult;
7182
7183	tempStr = (char ) SafeMalloc((size_t) (tempStrSize sizeof(char)));
7184	if (!tempStr)
7185	MrBayesPrint ("%s Problem allocating tempString (%d)\n", spacer, tempStrSize * sizeof(char));
7186
7187	if (p != NULL)
7188	{
7189	if (p->left == NULL && p->right == NULL)
7190	{
7191	SafeSprintf (&tempStr, &tempStrSize, "%d:%s", p->index + 1, MbPrintNum(p->length));
7192	AddToPrintString (tempStr);
7193	for (i=0; i<param->nSubParams; i++)
7194	{
7195	if (param->subParams[i]->paramType == P_CPPEVENTS)
7196	{
7197	nEvents = param->subParams[i]->nEvents[2*chain+state[chain]][p->index];
7198	if (nEvents > 0)
7199	{
7200	SafeSprintf (&tempStr, &tempStrSize, "[&E %s %d", param->subParams[i]->name, nEvents);
7201	AddToPrintString (tempStr);
7202	position = param->subParams[i]->position[2*chain+state[chain]][p->index];
7203	rateMult = param->subParams[i]->rateMult[2*chain+state[chain]][p->index];
7204	if (printAll == YES)
7205	{
7206	SafeSprintf (&tempStr, &tempStrSize, ": (");
7207	AddToPrintString (tempStr);
7208	for (j=0; j<nEvents; j++)
7209	{
7210	SafeSprintf (&tempStr, &tempStrSize, "%s", MbPrintNum(position[j]));
7211	AddToPrintString (tempStr);
7212	SafeSprintf (&tempStr, &tempStrSize, " %s", MbPrintNum(rateMult[j]));
7213	AddToPrintString (tempStr);
7214	if (j != nEvents-1)
7215	AddToPrintString (",");
7216	else
7217	AddToPrintString (")");
7218	}
7219	}
7220	AddToPrintString ("]");
7221	}
7222	else
7223	{
7224	SafeSprintf (&tempStr, &tempStrSize, "[&E %s 0]", param->subParams[i]->name);
7225	AddToPrintString (tempStr);
7226	}
7227	}
7228	else if (param->subParams[i]->paramType != P_CPPEVENTS)
7229	{
7230	/* other relaxed clock models */
7231	brlen = GetParamSubVals (param->subParams[i], chain, state[chain])[p->index];
7232	SafeSprintf (&tempStr, &tempStrSize, "[&B %s %s]", param->subParams[i]->name, MbPrintNum(brlen));
7233	AddToPrintString (tempStr);
7234	}
7235	}
7236	}
7237	else
7238	{
7239	if (p->anc != NULL)
7240	AddToPrintString ("(");
7241	WriteEventTreeToPrintString (p->left, chain, param, printAll);
7242	AddToPrintString (",");
7243	WriteEventTreeToPrintString (p->right, chain, param, printAll);
7244	if (p->anc != NULL)
7245	{
7246	if (p->anc->anc != NULL)
7247	{
7248	SafeSprintf (&tempStr, &tempStrSize, "):%s", MbPrintNum(p->length));
7249	AddToPrintString (tempStr);
7250	for (i=0; i<param->nSubParams; i++)
7251	{
7252	if (param->subParams[i]->paramType == P_CPPEVENTS)
7253	{
7254	nEvents = param->subParams[i]->nEvents[2*chain+state[chain]][p->index];
7255	if (nEvents > 0)
7256	{
7257	SafeSprintf (&tempStr, &tempStrSize, "[&E %s %d", param->subParams[i]->name, nEvents);
7258	AddToPrintString (tempStr);
7259	position = param->subParams[i]->position[2*chain+state[chain]][p->index];
7260	rateMult = param->subParams[i]->rateMult[2*chain+state[chain]][p->index];
7261	if (printAll == YES)
7262	{
7263	SafeSprintf (&tempStr, &tempStrSize, ": (");
7264	AddToPrintString (tempStr);
7265	for (j=0; j<nEvents; j++)
7266	{
7267	SafeSprintf (&tempStr, &tempStrSize, "%s", MbPrintNum(position[j]));
7268	AddToPrintString (tempStr);
7269	SafeSprintf (&tempStr, &tempStrSize, " %s", MbPrintNum(rateMult[j]));
7270	AddToPrintString (tempStr);
7271	if (j != nEvents-1)
7272	AddToPrintString (",");
7273	else
7274	AddToPrintString (")");
7275	}
7276	}
7277	AddToPrintString ("]");
7278	}
7279	else
7280	{
7281	SafeSprintf (&tempStr, &tempStrSize, "[&E %s 0]", param->subParams[i]->name);
7282	AddToPrintString (tempStr);
7283	}
7284	}
7285	else if (param->subParams[i]->paramType != P_CPPEVENTS)
7286	{
7287	/* other relaxed clock models */
7288	brlen = GetParamSubVals (param->subParams[i], chain, state[chain])[p->index];
7289	SafeSprintf (&tempStr, &tempStrSize, "[&B %s %s]", param->subParams[i]->name, MbPrintNum(brlen));
7290	AddToPrintString (tempStr);
7291	}
7292	}
7293	}
7294	else
7295	AddToPrintString(")");
7296	}
7297	}
7298	}
7299	free (tempStr);
7300	}
7301
7302
7303
7304
7305
7306	void WriteEvolTree (TreeNode p, int chain, Param param)
7307
7308	{
7309	MrBFlt *length;
7310
7311	if (p != NULL)
7312	{
7313	length = GetParamSubVals(param, chain, state[chain]);
7314	if (p->left == NULL && p->right == NULL)
7315	{
7316	printf ("%d:%s", p->index + 1, MbPrintNum(length[p->index]));
7317	}
7318	else
7319	{
7320	if (p->anc != NULL)
7321	printf ("(");
7322	WriteEvolTree(p->left, chain, param);
7323	printf (",");
7324	WriteEvolTree(p->right, chain, param);
7325	if (p->anc != NULL)
7326	{
7327	if (p->anc->anc != NULL)
7328	printf ("):%s", MbPrintNum(length[p->index]));
7329	else
7330	printf(")");
7331	}
7332	}
7333	}
7334	}
7335
7336
7337
7338
7339
7340	void WriteTreeToPrintString (Param param, int chain, TreeNode p, int showBrlens, int isRooted)
7341
7342	{
7343	char *tempStr;
7344	int i, tempStrSize = TEMPSTRSIZE, nEvents;
7345	MrBFlt brlen, N;
7346
7347	tempStr = (char ) SafeMalloc((size_t) (tempStrSize sizeof(char)));
7348	if (!tempStr)
7349	MrBayesPrint ("%s Problem allocating tempString (%d)\n", spacer, tempStrSize * sizeof(char));
7350
7351	if (p != NULL)
7352	{
7353	if (p->left == NULL && p->right == NULL)
7354	{
7355	if (showBrlens == YES)
7356	{
7357	SafeSprintf (&tempStr, &tempStrSize, "%d:%s", p->index + 1, MbPrintNum(p->length));
7358	}
7359	else
7360	SafeSprintf (&tempStr, &tempStrSize, "%d", p->index + 1);
7361	AddToPrintString (tempStr);
7362	if (param->paramType == P_BRLENS)
7363	{
7364	for (i=0; i<param->nSubParams; i++)
7365	{
7366	if (param->subParams[i]->paramType == P_CPPEVENTS)
7367	{
7368	nEvents = param->subParams[i]->nEvents[2*chain+state[chain]][p->index];
7369	SafeSprintf (&tempStr, &tempStrSize, "[&E %s %d]", param->subParams[i]->name, nEvents);
7370	AddToPrintString (tempStr);
7371	}
7372	brlen = GetParamSubVals (param->subParams[i], chain, state[chain])[p->index];
7373	SafeSprintf (&tempStr, &tempStrSize, "[&B %s %s]", param->subParams[i]->name, MbPrintNum(brlen));
7374	AddToPrintString (tempStr);
7375	}
7376	}
7377	else if (param->paramType == P_SPECIESTREE && modelSettings[param->relParts[0]].popSize->nValues > 1)
7378	{
7379	N = GetParamVals (modelSettings[param->relParts[0]].popSize, chain, state[chain])[p->index];
7380	SafeSprintf (&tempStr, &tempStrSize, "[&N %s %s]", modelSettings[param->relParts[0]].popSize->name, MbPrintNum(N));
7381	AddToPrintString (tempStr);
7382	}
7383	}
7384	else
7385	{
7386	if (p->anc != NULL)
7387	AddToPrintString ("(");
7388	WriteTreeToPrintString (param, chain, p->left, showBrlens, isRooted);
7389	if (p->anc != NULL)
7390	AddToPrintString (",");
7391	WriteTreeToPrintString (param, chain, p->right, showBrlens, isRooted);
7392	if (p->anc != NULL)
7393	{
7394	if (p->anc->anc == NULL && isRooted == NO)
7395	{
7396	if (showBrlens == YES)
7397	SafeSprintf (&tempStr, &tempStrSize, ",%d:%s)", p->anc->index + 1, MbPrintNum(p->length));
7398	else
7399	SafeSprintf (&tempStr, &tempStrSize, ",%d)", p->anc->index + 1);
7400	AddToPrintString (tempStr);
7401	}
7402	else if (p->anc->anc != NULL)
7403	{
7404	if (showBrlens == YES)
7405	SafeSprintf (&tempStr, &tempStrSize, "):%s", MbPrintNum(p->length));
7406	else
7407	SafeSprintf (&tempStr, &tempStrSize, ")");
7408	AddToPrintString (tempStr);
7409	if (param->paramType == P_BRLENS)
7410	{
7411	for (i=0; i<param->nSubParams; i++)
7412	{
7413	if (param->subParams[i]->paramType == P_CPPEVENTS)
7414	{
7415	nEvents = param->subParams[i]->nEvents[2*chain+state[chain]][p->index];
7416	SafeSprintf (&tempStr, &tempStrSize, "[&E %s %d]", param->subParams[i]->name, nEvents);
7417	AddToPrintString (tempStr);
7418	}
7419	brlen = GetParamSubVals (param->subParams[i], chain, state[chain])[p->index];
7420	SafeSprintf (&tempStr, &tempStrSize, "[&B %s %s]", param->subParams[i]->name, MbPrintNum(brlen));
7421	AddToPrintString (tempStr);
7422	}
7423	}
7424	else if (param->paramType == P_SPECIESTREE && modelSettings[param->relParts[0]].popSize->nValues > 1)
7425	{
7426	N = GetParamVals (modelSettings[param->relParts[0]].popSize, chain, state[chain])[p->index];
7427	SafeSprintf (&tempStr, &tempStrSize, "[&N %s %s]", modelSettings[param->relParts[0]].popSize->name, MbPrintNum(N));
7428	AddToPrintString (tempStr);
7429	}
7430	}
7431	else if (param->paramType == P_SPECIESTREE && modelSettings[param->relParts[0]].popSize->nValues > 1)
7432	{
7433	N = GetParamVals (modelSettings[param->relParts[0]].popSize, chain, state[chain])[p->index];
7434	SafeSprintf (&tempStr, &tempStrSize, ")[&N %s %s]", modelSettings[param->relParts[0]].popSize->name, MbPrintNum(N));
7435	AddToPrintString (tempStr);
7436	}
7437	else
7438	AddToPrintString(")");
7439	}
7440	}
7441	}
7442	free (tempStr);
7443	}
7444
7445
7446
7447
7448
7449	/* WriteTopologyToFile: Simply write topology to file */
7450	void WriteTopologyToFile (FILE fp, TreeNode p, int isRooted)
7451
7452	{
7453	if (p != NULL)
7454	{
7455	if (p->left == NULL && p->right == NULL)
7456	fprintf (fp, "%d", p->index + 1);
7457	else
7458	{
7459	if (p->anc != NULL)
7460	fprintf (fp, "(");
7461	WriteTopologyToFile (fp, p->left, isRooted);
7462	if (p->anc != NULL)
7463	fprintf (fp, ",");
7464	WriteTopologyToFile (fp, p->right, isRooted);
7465	if (p->anc != NULL)
7466	{
7467	if (p->anc->anc == NULL && isRooted == NO)
7468	fprintf (fp, ",%d", p->anc->index + 1);
7469	fprintf (fp, ")");
7470	}
7471	}
7472	}
7473	}

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: trunk/GDE/MrBAYES/mrbayes_3.2.1/tree.c

Download in other formats: