1 /*------------------------------------------------------
2 Maximum likelihood estimation
3 of migration rate and effective population size
4 using a Metropolis-Hastings Monte Carlo algorithm
5 -------------------------------------------------------
6 M C M C R O U T I N E S
7
8 Markov Monte Carlo stuff: treechange, acceptance
9
10
11 Peter Beerli 1996, Seattle
12 beerli@fsu.edu
13
14 Copyright 1996-2002 Peter Beerli and Joseph Felsenstein, Seattle WA
15 Copyright 2003-2013 Peter Beerli, Tallahassee FL
16
17 Permission is hereby granted, free of charge, to any person obtaining
18 a copy of this software and associated documentation files (the
19 "Software"), to deal in the Software without restriction, including
20 without limitation the rights to use, copy, modify, merge, publish,
21 distribute, sublicense, and/or sell copies of the Software, and to
22 permit persons to whom the Software is furnished to do so, subject
23 to the following conditions:
24
25 The above copyright notice and this permission notice shall be
26 included in all copies or substantial portions of the Software.
27
28 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
29 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
30 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
31 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
32 ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
33 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
34 WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
35
36 $Id: mcmc1.c 2169 2013-08-24 19:02:04Z beerli $
37 -------------------------------------------------------*/
38 /* \file mcmc1.c
39
40 Tree changer and acceptance rejection scheme
41
42 */
43 #include "migration.h"
44 #include "sighandler.h"
45 #include "tools.h"
46 #include "random.h"
47 #include "tree.h"
48 #include "mcmc2.h"
49 #ifdef UEP
50 #include "uep.h"
51 #endif
52
53 #include "bayes.h"
54
55 #ifdef BEAGLE
56 #include "calculator.h"
57 #endif
58
59 #ifdef DMALLOC_FUNC_CHECK
60 #include <dmalloc.h>
61 #endif
62
63 #define MIGRATION_AIR (boolean) 1
64 #define MIGRATION_IN_TREE (boolean) 0
65 #define NO_MIGR_NODES 0
66 #define WITH_MIGR_NODES 1
67
68 boolean debugtip;
69
70 extern int myID;
71 extern MYREAL probg_treetimesX(world_fmt *world, vtlist *tl, long T);
72
73 extern void zero_xseq(xarray_fmt *x, long sites, long categs);
74
75 /* prototypes ------------------------------------------- */
76 // metropolize over trees
77 long tree_update (world_fmt * world, long g);
78 /* private functions */
79 void new_localtimelist (timelist_fmt ** ntl, timelist_fmt * otl, long numpop);
80 void new_proposal (proposal_fmt ** proposal, timelist_fmt * tl,
81 world_fmt * world);
82 void set_new_proposal (proposal_fmt ** proposal, timelist_fmt * tl, world_fmt * world);
83 void chooseOrigin (proposal_fmt * proposal);
84 void construct_localtimelist (timelist_fmt * timevector,
85 proposal_fmt * proposal);
86 void traverseAllNodes (node * theNode, node *** nodelist, long *node_elem,
87 long *oldnode_elem, int include_migration);
88 void chooseTarget (proposal_fmt * proposal, timelist_fmt * timevector,
89 node ** bordernodes, long *bordernum);
90 void findbordernodes (node * theNode, proposal_fmt * proposal, long pop,
91 node *** bordernodes, long *allocsize, long *bordernum, vtlist ** tyme,
92 long gte);
93 void free_masterproposal (proposal_fmt * proposal);
94 void free_timevector (timelist_fmt * timevector);
95 void prune_timelist (timelist_fmt * oldtv, timelist_fmt * newtv, register node ** __restrict ptr, proposal_fmt * proposal, long numpop);
96
97 long xor (node ** ptrl1, node ** ptrl2);
98 long rmigrcount (proposal_fmt * proposal);
99
100 int migrate_old (proposal_fmt * proposal, node * up,
101 long *old_migr_table_counter, boolean air);
102 int migrate (proposal_fmt * proposal, node * up);
103 int migrateb (proposal_fmt * proposal, node * up);
104
105 int pre_population (proposal_fmt * proposal, vtlist * ltime, long gte,
106 long *slider);
107
108 boolean acceptlike (world_fmt * world, proposal_fmt * proposal, long g,
109 timelist_fmt * tyme);
110 MYREAL eventtime (proposal_fmt * proposal, long pop, vtlist * tentry,
111 char *event);
112 node *showsister (node * theNode);
113 void count_migrations (node * p, long *count);
114 long migration_from (long to, proposal_fmt * proposal);
115
116 MYREAL prob_tree (world_fmt * world, timelist_fmt * tyme);
117 void traverse_check(node *theNode);
118 void reset_proposal (proposal_fmt ** proposal, world_fmt *world);
119
120
121 /* Functions ++++++++++++++++++++++++++++++++++++++++++++++++*/
122 void
set_tree_dirty(node * p)123 set_tree_dirty (node * p)
124 {
125 switch (p->type)
126 {
127 case 'm':
128 set_dirty (p);
129 set_tree_dirty (p->next->back);
130 break;
131 case 't':
132 break;
133 case 'i':
134 set_dirty (p);
135 set_tree_dirty (p->next->back);
136 set_tree_dirty (p->next->next->back);
137 break;
138 case 'r':
139 set_dirty (p);
140 set_tree_dirty (p->next->back);
141 break;
142 }
143 }
144
145 /*=======================================================*/
146 long
tree_update(world_fmt * world,long g)147 tree_update (world_fmt * world, long g)
148 { /*return 1 if tree was accepted, 0 otherwise */
149 // static long treefilepos; /* write position in the treefile */
150 #ifdef BEAGLE
151 int beagle_control=1;
152 #endif
153 //short show=0;
154 boolean coalesced;
155 // boolean test;
156 char event;
157 long slider;
158 long bordernum;
159 long actualpop = -99, zz;
160 MYREAL endtime, nexttime, age;
161 #ifdef UEP
162
163 boolean uepsuccess = FALSE;
164 #endif
165 #ifdef TESTING2
166 proposal_fmt *proposal = world->proposal;
167 #else
168 proposal_fmt *proposal=NULL;
169 #endif
170 /*scratchpad in which all involved nodes
171 are recorded and help arrays, e.g. migration arrays,
172 are stored */
173 timelist_fmt *timevector = NULL; /* local timelist */
174 vtlist *tentry = NULL; /*pointer into timeslice */
175
176 /* ---------------------------------------------------------
177 initialize local timelist and construct residual timelist
178 find the start node and snip of the branch and node below */
179 #ifdef __MWERKS__
180
181 eventloop ();
182 #endif
183 new_localtimelist (&timevector, &world->treetimes[0], world->numpop);
184 #ifdef TESTING2
185 if(world->has_proposal_details)
186 reset_proposal (&proposal, world);//, &world->treetimes[0], world);
187 else
188 {
189 set_new_proposal (&proposal, &world->treetimes[0], world);
190 world->has_proposal_details = TRUE;
191 }
192 #else
193 new_proposal (&proposal, &world->treetimes[0], world);
194 #endif
195 chooseOrigin (proposal);
196 //
197 if(world->options->bayes_infer)
198 {
199 world->bayes->starttime = -proposal->origin->tyme; // we start with negative value to
200 // indicate that we do not know yet whether we will accept this move or not
201 // if the move is accepted then the time will be negated (->positiv), we record
202 // all start and stop times one could easily imagine that some parameter values
203 // never gets accepted by the genealogy, and when we record this for bayesallfile we would
204 // like to know this. see further doen at the acceptance for the proposals
205 }
206 //
207 construct_localtimelist (timevector, proposal);
208 tentry = &(*timevector).tl[0];
209 age = proposal->origin->tyme;
210 zz = 0;
211 while ((tentry->age < age || tentry->age - age < SMALLEPSILON)&& zz < (*timevector).T)
212 {
213 tentry = &(*timevector).tl[zz];
214 zz++;
215 }
216 zz--;
217 nexttime = tentry->age;
218 if ((*timevector).T > 1)
219 endtime = (*timevector).tl[(*timevector).T - 2].age;
220 else
221 endtime = 0.0;
222 proposal->time = 0.0;
223 #ifdef TREEDEBUG
224 fprintf(stdout,"%i> begin zz=%li nextime=%f (pt=%f+age=%f)=%f,pot=%f\n",myID, zz, nexttime, proposal->time, age, age+proposal->time,proposal->origin->tyme);
225 #endif
226 proposal->time = age;
227 coalesced = FALSE;
228 /*------------------------------------
229 main loop: sliding down the tree */
230 slider = 0;
231 #ifdef DEBUG
232 //printf("\nTime Pop Type Lineages [mutrate=%f]\n",proposal->world->options->mu_rates[world->locus]);
233 #endif
234 while (nexttime <= endtime)
235 {
236 actualpop =
237 (proposal->migr_table_counter >
238 0) ? proposal->migr_table[proposal->migr_table_counter -
239 1].from : proposal->origin->pop;
240 #ifdef TREEDEBUG
241 if(age == 0.0 && proposal->origin->tyme>0.0)
242 {
243 printf("\n\n%f\n\n",age);
244 }
245 #endif
246 proposal->time = eventtime (proposal, actualpop, tentry, &event);
247 #ifdef TREEDEBUG
248 fprintf(stdout,"%i> zz=%li (pt=%f+age=%f)=%f,pot=%f\n",myID, zz, proposal->time, age, age+proposal->time,proposal->origin->tyme);
249 #endif
250 proposal->time += age;
251 // age = proposal->time;
252
253 if(proposal->time < proposal->origin->tyme)
254 {
255 fprintf(stdout,"%i> Proposal failed because of unordered entry in time list, abort this sample\nProposed time=%f origin time=%f nexttime=%f\n", myID, proposal->time, proposal->origin->tyme, nexttime);
256 // we end up here when the migration events exceed the upper limit
257 #ifndef TESTING2
258 free_masterproposal (proposal);
259 #endif
260 free_timevector (timevector);
261 return 0;
262 // error("Proposal of time failed\n");
263 }
264 if (proposal->time < nexttime)
265 {
266 if (event == 'm')
267 {
268 if (!migrate (proposal, proposal->origin))
269 {
270 // we end up here when the migration events exceed the upper limit
271 #ifndef TESTING2
272 free_masterproposal (proposal);
273 #endif
274 free_timevector (timevector);
275 //warning("Event was migration event but could not match a lineage to it\n");
276 return 0;
277 }
278 age = proposal->time;
279 continue;
280 }
281 else
282 { /*coalesce */
283 chooseTarget (proposal, timevector, proposal->bordernodes,
284 &bordernum);
285 if(bordernum == 0)
286 {
287 if (!migrate (proposal, proposal->origin))
288 {
289 // we end up here when the migration events exceed the upper limit
290 #ifndef TESTING2
291 free_masterproposal (proposal);
292 #endif
293 free_timevector (timevector);
294 return 0;
295 }
296 age = proposal->time;
297 continue;
298 // warning("chooseTarget failed at age %f\n", proposal->time);
299 // free_proposal (proposal);
300 // free_timevector (timevector);
301 // return 0;
302 }
303 pretendcoalesce1p (proposal);
304 #ifdef UEP
305 uepsuccess = is_success_pseudo_uep (proposal);
306 #endif
307
308 coalesced = TRUE;
309 break;
310 }
311 } /*end if proposal->time < nextime */
312 age = nexttime;
313 #ifdef TREEDEBUG
314 if(age < proposal->origin->tyme)
315 {
316 printf("darn: nextime=%f\n",nexttime);
317 }
318 #endif
319 zz++;
320 if(zz >= timevector->T)
321 {
322 break;
323 }
324 tentry = &(*timevector).tl[zz]; /*next entry in timelist */
325 nexttime = tentry->age;
326 }
327 if (!coalesced)
328 {
329 if (!pre_population(proposal, (*timevector).tl, (*timevector).T - 1, &slider))
330 {
331 #ifndef TESTING2
332 free_masterproposal (proposal);
333 #endif
334 free_timevector (timevector);
335 return 0;
336 }
337 // fprintf(stderr,"%i> last time=%f, end time = %f (%c)\n", myID, age,(*timevector).tl[(*timevector).T-1].eventnode->tyme,(*timevector).tl[(*timevector).T-1].eventnode->type );
338 pretendcoalesce1p (proposal);
339 #ifdef UEP
340
341 uepsuccess = is_success_pseudo_uep (proposal);
342 #endif
343
344 }
345 if (
346 #ifdef UEP
347 ((!world->options->uep && !uepsuccess)
348 || (world->options->uep && uepsuccess)) &&
349 #endif
350 (acceptlike (world, proposal, g, timevector)))
351 {
352 #ifdef BEAGLE
353 printf("Accepted\n");
354 change_beagle(world->root->next->back,world->beagle,world->sumtips);
355 #endif
356 if (proposal->time > world->root->tyme)
357 { /*saveguard */
358 world->root->tyme += proposal->time;
359 }
360 coalesce1p (proposal);
361 //DEBUG
362 //traverse_check(crawlback (proposal->root->next));printf("*");
363 #ifdef DEBUG
364 //traverse_adjust(world->root->next->back, 1.0);
365 #endif
366 // record the time interval that was used for the lineage.
367 if(world->options->bayes_infer)
368 {
369 world->bayes->starttime = -world->bayes->starttime;
370 world->bayes->stoptime = proposal->time;
371 world->treelen = 0.0;
372 calc_treelength (world->root->next->back, &world->treelen);
373 }
374 //
375 #ifdef UEP
376 world->likelihood[g] = treelikelihood (world);
377 if (world->options->uep)
378 {
379 update_uep (world->root->next->back, world);
380 check_uep_root (world->root->next->back, world);
381 world->treelen = 0.0;
382 calc_treelength (world->root->next->back, &world->treelen);
383 world->ueplikelihood = ueplikelihood (world);
384 world->likelihood[g] = /*world->ueplikelihood +*/ world->likelihood[g];
385 }
386 #else
387 #ifdef BEAGLE
388 // set_beagle_dirty(proposal->origin,proposal->target,showtop(world->root->next->back));
389 //reset_beagle(world->beagle);
390 //smooth (world->root->next, crawlback(world->root->next), world, world->locus);
391 //if(world->beagle->numbranches>0)
392 // {
393 world->likelihood[g] = proposal->likelihood; //treelikelihood (world);
394 if(beagle_control==1)
395 {
396 printf("Recalculate LnL using treelikelihood()\n");
397 world->likelihood[g] = force_beagle_recalculate(world,world->locus);
398 printf("%i> %f = newLnL=%f + recalcLnL=%f\n",myID, proposal->likelihood - world->likelihood[g], proposal->likelihood,world->likelihood[g]);
399 }
400 // }
401 #else
402 world->likelihood[g] = treelikelihood (world);
403 #endif /*BEAGLE*/
404 #endif /*UEP*/
405 /* create a new timelist */
406 construct_tymelist (world, &world->treetimes[0]);
407 // if (world->options->treeprint != _NONE && world->cold)
408 // print_tree (world, g, &treefilepos);
409 world->migration_counts = 0;
410 /* report the number of migration on the tree */
411 count_migrations (world->root->next->back, &world->migration_counts);
412 #ifndef TESTING2
413 free_masterproposal (proposal);
414 #endif
415 free_timevector (timevector);
416 return 1; /* new tree accepted */
417 }
418 else
419 {
420 // record the time interval that was used for the lineage.
421 if(world->options->bayes_infer)
422 {
423 world->bayes->stoptime = -proposal->time;
424 }
425 //traverse_check(crawlback (proposal->root->next));printf(".");
426 //
427 }
428 #ifndef TESTING2
429 free_masterproposal (proposal);
430 #endif
431 free_timevector (timevector);
432 return 0; /* not accepted */
433 }
434
435 #ifdef SLATKIN_IMPORTANCE
slatkin_importance(world_mt * world,long g)436 long slatkin_importance(world_mt *world, long g)
437 {
438 long slice=0;
439 // create timelist with times drawn using the parameters
440 new_localtimelist (&timevector, &world->treetimes[0], world->numpop);
441 while (slice < sumlineages)
442 {
443 actualpop =
444 (proposal->migr_table_counter >
445 0) ? proposal->migr_table[proposal->migr_table_counter -
446 1].from : proposal->origin->pop;
447 age = (*timevector).tl[slice].age;
448 (*timevector).tl[slice].age = age + eventtime (proposal, actualpop, tentry, &event);
449 }
450 // calculate which two datapoints to join give the distance
451 }
452 #endif
453
454 /*=======================================================*/
455 ///
456 /// allocate the timelist, contains times and lineages at that time
457 /// lineages have a large storage device that is accessed from the
458 /// tl[i].lineages.
459 void
new_localtimelist(timelist_fmt ** ntl,timelist_fmt * otl,long numpop)460 new_localtimelist (timelist_fmt ** ntl, timelist_fmt * otl, long numpop)
461 {
462 // long i;
463 (*ntl) = (timelist_fmt *) mycalloc (1, sizeof (timelist_fmt));
464 (*ntl)->tl = (vtlist *) mymalloc ((*otl).allocT * sizeof (vtlist));
465 (*ntl)->allocT = otl->allocT;
466 (*ntl)->T = otl->T;
467 (*ntl)->oldT = otl->oldT;
468 //memcpy ((*ntl)->tl, otl->tl, otl->allocT * sizeof (vtlist));
469 allocate_lineages (ntl, 0, numpop);
470 //memcpy ((*ntl)->lineages, otl->lineages, (*ntl)->allocT * numpop * sizeof (long));
471 }
472
473 ///
474 /// reuse a global timelist to store values, this should save
475 /// malloc/free calls
476 void
new_localtimelist_new(timelist_fmt ** ntl,timelist_fmt * otl,long numpop)477 new_localtimelist_new (timelist_fmt ** ntl, timelist_fmt * otl, long numpop)
478 {
479 // UNFINISHED
480 // long i;
481 // (*ntl) = (timelist_fmt *) mycalloc (1, sizeof (timelist_fmt));
482 if((*ntl)->allocT < otl->allocT)
483 {
484 (*ntl)->tl = (vtlist *) myrealloc ((*ntl)->tl,(*otl).allocT * sizeof (vtlist));
485 // memset((*ntl)->tl+((*ntl)->allocT),0,sizeof(vtlist)*((*otl)->allocT-(*ntl)->allocT));
486 (*ntl)->allocT = otl->allocT;
487 }
488 (*ntl)->T = otl->T;
489 memcpy ((*ntl)->tl, otl->tl, otl->allocT * sizeof (vtlist));
490 allocate_lineages (ntl, 0, numpop);
491 memcpy ((*ntl)->lineages, otl->lineages, (*ntl)->allocT * numpop * sizeof (long));
492 }
493
494
495
496 void
new_proposal(proposal_fmt ** proposal,timelist_fmt * tl,world_fmt * world)497 new_proposal (proposal_fmt ** proposal, timelist_fmt * tl, world_fmt * world)
498 {
499 #ifdef UEP
500 long j;
501 #endif
502 long mal = world->data->maxalleles[world->locus];
503 //long listsize = ((*tl).allocT + (*tl).T + 5);
504 long listsize = 2*(world->sumtips + 2);
505 long sumtips = world->sumtips;
506
507 // allocate the scratchpad (contains a shadow of the tree)
508 (*proposal) = (proposal_fmt *) mycalloc (1, sizeof (proposal_fmt));
509 (*proposal)->likelihood = -HUGE;
510 // pointers and values to outside structures
511 (*proposal)->world = world;
512 (*proposal)->datatype = world->options->datatype;
513 (*proposal)->sumtips = world->sumtips;
514 (*proposal)->numpop = world->numpop;
515 (*proposal)->endsite = world->data->seq[0]->endsite;
516 (*proposal)->fracchange = world->data->seq[0]->fracchange;
517 (*proposal)->param0 = world->param0;
518 (*proposal)->root = world->root;
519 (*proposal)->migration_model = world->options->migration_model;
520 // precalculated values
521 (*proposal)->mig0list = world->mig0list;
522 (*proposal)->design0list = world->design0list;
523
524 (*proposal)->listsize = listsize;
525
526 // nodes above the picked node + migration nodes
527 (*proposal)->aboveorigin =
528 (node **) mycalloc (listsize, sizeof (node *));
529
530 // node data holding vector for bordernodes, line_f, line_t
531 (*proposal)->nodedata =
532 (node **) mycalloc (3 * listsize, sizeof (node *));
533 // adjacent nodes of the picked nodes
534 (*proposal)->bordernodes = (*proposal)->nodedata;
535 // line_f nodes
536 (*proposal)->line_f = (*proposal)->bordernodes + listsize;
537 (*proposal)->line_t = (*proposal)->line_f + listsize;
538 // mf holds also mt array
539 (*proposal)->mf = (MYREAL *) mycalloc (2* (*proposal)->endsite, sizeof (MYREAL));
540 (*proposal)->mt = (*proposal)->mf + (*proposal)->endsite;
541
542 if (strchr (SEQUENCETYPES, (*proposal)->datatype))
543 {
544 allocate_xseq(&(*proposal)->xf, world->data->seq[0]->endsite, world->options->rcategs);
545 allocate_xseq(&(*proposal)->xt, world->data->seq[0]->endsite, world->options->rcategs);
546 }
547 else
548 {
549 (*proposal)->xf.a = (MYREAL *) mycalloc (1, sizeof (MYREAL) * mal);
550 (*proposal)->xt.a = (MYREAL *) mycalloc (1, sizeof (MYREAL) * mal);
551 }
552 (*proposal)->old_migr_table_counter = 4 * sumtips /* 100 */ ;
553 (*proposal)->old_migr_table_counter2 = 4 * sumtips /* 100 */ ;
554 (*proposal)->migr_table =
555 (migr_table_fmt *) mycalloc (1,
556 sizeof (migr_table_fmt) *
557 (*proposal)->old_migr_table_counter);
558 (*proposal)->migr_table2 =
559 (migr_table_fmt *) mycalloc ((*proposal)->old_migr_table_counter2,
560 sizeof (migr_table_fmt));
561 (*proposal)->migr_table_counter = 0;
562 (*proposal)->migr_table_counter2 = 0;
563 #ifdef UEP
564 if (world->options->uep)
565 {
566 (*proposal)->ueplike =
567 (MYREAL **) mycalloc (world->data->uepsites, sizeof (MYREAL *));
568 (*proposal)->ueplike[0] =
569 (MYREAL *) mycalloc (world->numpop * world->data->uepsites,
570 sizeof (MYREAL));
571 for (j = 1; j < world->data->uepsites; ++j)
572 (*proposal)->ueplike[j] = (*proposal)->ueplike[0] + j * world->numpop;
573
574 (*proposal)->ut.s = (pair *) mycalloc (world->data->uepsites, sizeof (pair));
575 (*proposal)->uf.s = (pair *) mycalloc (world->data->uepsites, sizeof (pair));
576 (*proposal)->umt = (MYREAL *) mycalloc (world->data->uepsites, sizeof (MYREAL));
577 (*proposal)->umf = (MYREAL *) mycalloc (world->data->uepsites, sizeof (MYREAL));
578 }
579 #endif
580
581 #ifdef BEAGLE
582 (*proposal)->leftid = 0;
583 (*proposal)->rightid = 0;
584 (*proposal)->parentid = 0;
585 reset_beagle(world->beagle);
586 #endif
587 }
588
589 ///
590 /// sets new proposal structure using template gproposal, this function will be called most of the time instead of new_proposal()
591 void
set_new_proposal(proposal_fmt ** proposal,timelist_fmt * tl,world_fmt * world)592 set_new_proposal (proposal_fmt ** proposal, timelist_fmt * tl, world_fmt * world)
593 {
594 long mal = world->data->maxalleles[world->locus];
595 //long oldsize=0;
596 long newsize =0;
597 long sumtips = world->sumtips;
598 #ifdef UEP
599 long j;
600 #endif
601 long listsize = 2*(world->sumtips + 2);
602 (*proposal)->likelihood = -HUGE;
603 (*proposal)->sumtips = sumtips;
604 // pointers and values to outside structures
605 (*proposal)->world = world;
606 (*proposal)->datatype = world->options->datatype;
607 (*proposal)->numpop = world->numpop;
608 (*proposal)->endsite = world->data->seq[0]->endsite;
609 (*proposal)->fracchange = world->data->seq[0]->fracchange;
610 (*proposal)->param0 = world->param0;
611 (*proposal)->root = world->root;
612 (*proposal)->migration_model = world->options->migration_model;
613 // precalculated values
614 (*proposal)->mig0list = world->mig0list;
615 (*proposal)->design0list = world->design0list;
616 newsize = 4 * listsize;
617 if((*proposal)->nodedata == NULL)
618 (*proposal)->nodedata = (node **) mycalloc (newsize, sizeof (node *));
619 else
620 (*proposal)->nodedata = (node **) myrealloc ((*proposal)->nodedata, newsize * sizeof (node *));
621 //xcode oldsize = (*proposal)->listsize;
622 memset((*proposal)->nodedata, 0, newsize * sizeof (node *));
623 (*proposal)->aboveorigin = (*proposal)->nodedata;
624 (*proposal)->bordernodes = (*proposal)->nodedata + listsize;
625 (*proposal)->line_f = (*proposal)->bordernodes + listsize;
626 (*proposal)->line_t = (*proposal)->line_f + listsize;
627 (*proposal)->listsize = listsize;
628 // mf holds also mt array
629 (*proposal)->mf = (MYREAL *) mycalloc (2* (*proposal)->endsite, sizeof (MYREAL));
630 (*proposal)->mt = (*proposal)->mf + (*proposal)->endsite;
631 if (strchr (SEQUENCETYPES, (*proposal)->datatype))
632 {
633 allocate_xseq(&(*proposal)->xf, world->data->seq[0]->endsite, world->options->rcategs);
634 allocate_xseq(&(*proposal)->xt, world->data->seq[0]->endsite, world->options->rcategs);
635 }
636 else
637 {
638 (*proposal)->xf.a = (MYREAL *) mycalloc (1, sizeof (MYREAL) * mal);
639 (*proposal)->xt.a = (MYREAL *) mycalloc (1, sizeof (MYREAL) * mal);
640 }
641 (*proposal)->old_migr_table_counter = 4 * sumtips /* 100 */ ;
642 (*proposal)->old_migr_table_counter2 = 4 * sumtips /* 100 */ ;
643 (*proposal)->migr_table =
644 (migr_table_fmt *) mycalloc (1,
645 sizeof (migr_table_fmt) *
646 (*proposal)->old_migr_table_counter);
647 (*proposal)->migr_table2 =
648 (migr_table_fmt *) mycalloc ((*proposal)->old_migr_table_counter2,
649 sizeof (migr_table_fmt));
650 (*proposal)->migr_table_counter = 0;
651 (*proposal)->migr_table_counter2 = 0;
652 #ifdef UEP
653 if (world->options->uep)
654 {
655 (*proposal)->ueplike =
656 (MYREAL **) mycalloc (world->data->uepsites, sizeof (MYREAL *));
657 (*proposal)->ueplike[0] =
658 (MYREAL *) mycalloc (world->numpop * world->data->uepsites,
659 sizeof (MYREAL));
660 for (j = 1; j < world->data->uepsites; ++j)
661 (*proposal)->ueplike[j] = (*proposal)->ueplike[0] + j * world->numpop;
662
663 (*proposal)->ut.s = (pair *) mycalloc (world->data->uepsites, sizeof (pair));
664 (*proposal)->uf.s = (pair *) mycalloc (world->data->uepsites, sizeof (pair));
665 (*proposal)->umt = (MYREAL *) mycalloc (world->data->uepsites, sizeof (MYREAL));
666 (*proposal)->umf = (MYREAL *) mycalloc (world->data->uepsites, sizeof (MYREAL));
667 }
668 #endif
669
670 #ifdef BEAGLE
671 (*proposal)->leftid = 0;
672 (*proposal)->rightid = 0;
673 (*proposal)->parentid = 0;
674 reset_beagle(world->beagle);
675 #endif
676 }
677
678
679 void
jumblenodes(node ** s,long n)680 jumblenodes (node ** s, long n)
681 {
682 node **temp;
683
684 long i, rr, tn = n;
685
686 temp = (node **) mycalloc (1, sizeof (node *) * n);
687 memcpy (temp, s, sizeof (node *) * n);
688 for (i = 0; i < n && tn > 0; i++)
689 {
690 s[i] = temp[rr = RANDINT (0L, tn - 1)];
691 temp[rr] = temp[tn - 1];
692 tn--;
693 }
694 myfree(temp);
695 }
696
traverse_check(node * theNode)697 void traverse_check(node *theNode)
698 {
699 if (theNode != NULL)
700 {
701 if (theNode->type != 't')
702 {
703 if (theNode->next->back != NULL)
704 {
705 if(theNode->tyme < showtop(theNode->next->back)->tyme)
706 {
707 printf("Problem in traverse_check: id=%li time=%f type=%c time-up-next=%f type_up=%c\n",
708 theNode->id, theNode->tyme, theNode->type, theNode->next->back->tyme, theNode->next->back->type);
709 error("time conflict");
710 }
711 traverse_check (theNode->next->back);
712 }
713 if (theNode->type != 'm' && theNode->next->next->back != NULL)
714 {
715 if(theNode->tyme < theNode->next->next->back->tyme)
716 {
717 printf("Problem in traverse_check: id=%li time=%f type=%c time-up-next-next=%f\n",
718 theNode->id, theNode->tyme, theNode->type, theNode->next->next->back->tyme);
719 error("time conflict");
720 }
721
722 traverse_check (theNode->next->next->back);
723 }
724 }
725 }
726 else
727 {
728 error("Node is NULL????");
729 }
730 }
731
732 ///
733 /// tags the lineage that was added the last time the tree was changed.
traverse_tagnew(node * theNode,node * origin)734 void traverse_tagnew(node *theNode, node *origin)
735 {
736 if (theNode != NULL)
737 {
738 theNode->visited = FALSE;
739 if (theNode->type != 't' && theNode->next->next->back != NULL)
740 {
741 if(theNode->top != 1)
742 theNode = showtop(theNode);
743 //error("");
744
745 traverse_tagnew (theNode->next->back, origin);
746 if (theNode->type != 'm' && theNode->next->next->back != NULL)
747 {
748 traverse_tagnew (theNode->next->next->back, origin);
749 }
750 }
751 if(theNode==origin)
752 {
753 while(theNode != NULL && theNode->type != 'r')
754 {
755 theNode = showtop(showtop(theNode)->back);
756 theNode->visited = TRUE;
757 }
758 }
759 }
760 else
761 {
762 error("a missing node pointer encountered, aborted");
763 }
764 }
765
766 ///
767 /// pick a node at random from the available list of nodes, this code will even
768 /// work with datasets where there are only two nodes.
769 void
chooseOrigin(proposal_fmt * proposal)770 chooseOrigin (proposal_fmt * proposal)
771 {
772 long elem = 0;
773 // oldelem = (proposal->sumtips * 2.);
774 node *tmp=NULL, **goal;
775 //032110 goal = (node **) mycalloc (oldelem, sizeof (node *));
776 //#ifdef DEBUG
777 //traverse_check(crawlback (proposal->root->next));
778 //#endif
779 //032110 traverseAllNodes (crawlback (proposal->root->next), &goal, &elem, &oldelem,
780 //032110 NO_MIGR_NODES);
781 goal = proposal->world->nodep;
782 elem = proposal->world->sumtips * 2;
783 switch(elem)
784 {
785 case 0:
786 error("problem with choosing a node for the tree change [choosOrigin()]");
787 break;
788 case 1:
789 case 2:
790 tmp = goal[0];
791 break;
792 default:
793 tmp = goal[RANDINT (0L, elem - 2)];
794 while(tmp->back->type =='r')
795 tmp = goal[RANDINT (0L, elem - 2)];
796 break;
797 }
798 //032110 myfree(goal);
799 proposal->origin = tmp;
800 if (proposal->origin != showtop (crawlback (proposal->root->next)))
801 {
802 proposal->oback = showtop (crawlback (proposal->origin));
803 proposal->osister = showsister (proposal->origin);
804 if (proposal->oback != showtop (crawlback (proposal->root->next)))
805 {
806 proposal->ocousin = showsister (proposal->oback);
807 }
808 else
809 {
810 proposal->ocousin = NULL;
811 }
812 }
813 if (proposal->origin == NULL)
814 error ("Designation of origin for branch removal failed");
815 }
816
817
818 ///
819 /// construct a list of ordered times with pointers to the tree
820 /// using the origin to prune the old time list
821 void
construct_localtimelist(timelist_fmt * timevector,proposal_fmt * proposal)822 construct_localtimelist (timelist_fmt * timevector, proposal_fmt * proposal)
823 {
824 #ifdef TREEDEBUG
825 double dif;
826 long ii;
827 #endif
828 long z = 0;
829 long oz = proposal->listsize;
830 // long numpop = timevector->numpop = proposal->numpop;
831 traverseAllNodes (crawlback (proposal->origin)->back,
832 &proposal->aboveorigin, &z, &oz, WITH_MIGR_NODES);
833 proposal->aboveorigin[z++] = proposal->oback;
834 proposal->aboveorigin[z] = NULL;
835 prune_timelist(&proposal->world->treetimes[0],timevector, proposal->aboveorigin, proposal, proposal->world->numpop);
836 add_partlineages(proposal->numpop, &timevector);
837 #ifdef TREEDEBUG
838 for(ii=0;ii<timevector->T-2;ii++)
839 {
840 printf("%i> ii=%li %0.5f (%li->%li) %c %3li",myID, ii, timevector->tl[ii].age,
841 timevector->tl[ii].from, timevector->tl[ii].to,
842 timevector->tl[ii].eventnode->type,timevector->tl[ii].lineages[0]);
843 long jj;
844 for(jj = 1;jj < proposal->numpop; jj++)
845 printf(" %3li",timevector->tl[ii].lineages[jj]);
846 printf(" backtyme:%f dif:%f\n",showtop(timevector->tl[ii].eventnode->back)->tyme,dif=showtop(timevector->tl[ii].eventnode->back)->tyme-timevector->tl[ii].eventnode->tyme);
847 if(dif<0.0)
848 {
849 error("shit");
850 }
851 }
852 #endif
853 }
854
855 /*----------------------------------------------------------------------------
856 finds all nodes in a tree starting at the root node and crawling up
857 to the tips in a recursive fashion, writing nodeptrs in the nodelist vector
858 the flag include_migration is 1 if we want to touch the migration nodes too,
859 otherwise =0 -> jump over the migration nodes. for convenience we define the
860 the macros NO_MIGR_NODES=0 and WITH_MIGR_NODES=1 in the treesetup.h file
861 PB 1995
862 */
863 void
traverseAllNodes(node * theNode,node *** nodelist,long * node_elem,long * oldnode_elem,int include_migration)864 traverseAllNodes (node * theNode, node *** nodelist, long *node_elem,
865 long *oldnode_elem, int include_migration)
866 {
867 static long counter=0;
868 long elem;
869 counter++;
870 if(theNode->type != 'r' && theNode->tyme > showtop(theNode->back)->tyme)
871 {
872 printf("%i> TTRAVERSETRAVERSETRAVERSETRAVERSE\nfunction was called with %li times: with node (%c) time %f and node (%c) time %f\n",
873 myID, counter, theNode->type, theNode->tyme, showtop(theNode->back)->type, showtop(theNode->back)->tyme);
874 error("traverseAllNodes() failed");
875 }
876
877 if (include_migration == NO_MIGR_NODES)
878 {
879 // nodelist needs to be at least twice as long as sumtips
880 // otherwise it will run out of reserved memory
881 if (theNode->type != 't')
882 {
883 if (crawlback (theNode->next) != NULL)
884 traverseAllNodes (crawlback (theNode->next), nodelist, node_elem,
885 oldnode_elem, NO_MIGR_NODES);
886 if (theNode->type != 'm' && crawlback (theNode->next->next) != NULL)
887 traverseAllNodes (crawlback (theNode->next->next), nodelist,
888 node_elem, oldnode_elem, NO_MIGR_NODES);
889 (*nodelist)[(*node_elem)] = theNode;
890 (*node_elem) += 1;
891 if (theNode->type == 'm')
892 {
893 error ("Migration node encountered?! and died!");
894 }
895 if(*node_elem > *oldnode_elem)
896 error("die on the spot, this should not happen. Failed in traverseAllNodes without migration nodes");
897 }
898 else
899 {
900 (*nodelist)[(*node_elem)] = theNode;
901 (*node_elem) += 1;
902 if(*node_elem > *oldnode_elem)
903 error("die on the spot, this should not happen. Failed in traverseAllNodes without migration nodes at tip");
904 }
905 }
906 else
907 {
908 if (theNode->type != 't')
909 {
910 if (theNode->next->back != NULL)
911 traverseAllNodes (theNode->next->back, nodelist, node_elem,
912 oldnode_elem, WITH_MIGR_NODES);
913 if (theNode->type != 'm' && theNode->next->next->back != NULL)
914 traverseAllNodes (theNode->next->next->back, nodelist, node_elem,
915 oldnode_elem, WITH_MIGR_NODES);
916 if ((*node_elem) == (*oldnode_elem - 2))
917 {
918 elem = (*oldnode_elem) + (*oldnode_elem);
919 (*nodelist) =
920 (node **) myrealloc ((*nodelist), sizeof (node *) * elem);
921 memset ((*nodelist) + (*oldnode_elem), 0,
922 sizeof (node *) * (elem - (*oldnode_elem)));
923 *oldnode_elem = elem;
924 }
925 (*nodelist)[(*node_elem)++] = theNode;
926 }
927 else
928 {
929 if ((*node_elem) == (*oldnode_elem - 2))
930 {
931 elem = (*oldnode_elem) + (*oldnode_elem);
932 (*nodelist) =
933 (node **) myrealloc ((*nodelist), sizeof (node *) * elem);
934 memset ((*nodelist) + (*oldnode_elem), 0,
935 sizeof (node *) * (elem - (*oldnode_elem)));
936 *oldnode_elem = elem;
937 }
938 (*nodelist)[(*node_elem)++] = theNode;
939 }
940 }
941 }
942
943 ///
944 /// copies elements from the old timelist (oldtv)
945 /// into the new local timeslist (newtv) checking whether they
946 /// are used after the removal of the residual tree that is above the origin (ptr)
947 void
prune_timelist(timelist_fmt * oldtv,timelist_fmt * newtv,register node ** __restrict ptr,proposal_fmt * proposal,long numpop)948 prune_timelist (timelist_fmt * oldtv, timelist_fmt * newtv, register node ** __restrict ptr, proposal_fmt * proposal, long numpop)
949 {
950 register long i = 0;
951 register long j = 0;
952 register node * thenode;
953 long slot = 0;
954 vtlist * tls;
955
956 slot=0;
957 // go through all slices in old time list and
958 // extract node pointers
959 for (i = 0; i < (*oldtv).T; i++)
960 {
961 j = 0;
962 thenode = (*oldtv).tl[i].eventnode;
963 while ((thenode != ptr[j]) && (ptr[j] != NULL))
964 {
965 j++;
966 }
967 // if the comparison reaches the end of the ptr list that holds an NULL element at the end
968 // then the node needs to be present in the new time list.
969 if (ptr[j] == NULL)
970 {
971 tls = &((*newtv).tl[slot]);
972 tls->eventnode = thenode ;
973 tls->age = (*oldtv).tl[i].age ;
974 tls->interval = (*oldtv).tl[i].interval ;
975 if(thenode!=NULL)
976 {
977 tls->from = thenode->pop;
978 tls->to = thenode->actualpop;
979 }
980 tls->slice = slot;
981 slot++;
982 }
983 }
984 (*newtv).T = slot;
985 // if ((*newtv).tl[(*newtv).T - 1].eventnode->type != 'r')
986 //{
987 // error ("Root not at the end of local timelist\n");
988 //}
989 // sort seems to be necessary with dated samples
990 qsort ((void *) (*newtv).tl, (*newtv).T, sizeof (vtlist), agecmp);
991 }
992
993 /* replaces nodepointers in list 1 with NULL if they are present in list 2
994 returns the first NULL slot in the array.
995 */
996 long
xor(node ** ptrl1,node ** ptrl2)997 xor (node ** ptrl1, node ** ptrl2)
998 {
999 long i = 0, j = 0, slot = -1;
1000 /* assumes that there is an NULL element at the end */
1001 for (i = 0; ptrl1[i] != NULL; j = 0, i++)
1002 {
1003 while ((ptrl1[i] != ptrl2[j]) && (ptrl2[j] != NULL))
1004 j++;
1005 if (ptrl2[j] != NULL)
1006 {
1007 if (slot == -1)
1008 slot = i;
1009 ptrl1[i] = NULL;
1010 }
1011 }
1012 return slot;
1013 }
1014
1015 /* migrate() fills the PROPOSAL->MIGRATION_TABLE in the tree or
1016 PROPOSAL->MIGRATION_TABLE2 when at the bottom of the tree
1017
1018 PROPOSAL proposal-scratchpad
1019 UP node above (younger)
1020 MIGR_TABLE_COUNTER migration array counter,
1021 will increase by one during execution
1022 AIR if true standard execution, if false updating the last
1023 lineage in the residual tree.
1024 */
1025 int
migrate(proposal_fmt * proposal,node * up)1026 migrate (proposal_fmt * proposal, node * up)
1027 {
1028 long tmp;
1029 long numpop = proposal->numpop;
1030 long i = proposal->migr_table_counter;
1031 migr_table_fmt *array = proposal->migr_table;
1032 if (i > MIGRATION_LIMIT * numpop || numpop < 2)
1033 {
1034 return 0;
1035 }
1036 if (i > 0)
1037 array[i].to = array[i - 1].from;
1038 else
1039 array[i].to = up->pop;
1040 tmp = migration_from (array[i].to, proposal);
1041
1042 if(tmp < numpop)
1043 array[i].from = tmp;
1044 else
1045 return 0;
1046
1047 //DEBUG
1048 // printf("i: %3li %li %li\n",i,proposal->migr_table[i].from,proposal->migr_table[i].to);
1049 array[i++].time = proposal->time;
1050 if(proposal->time < proposal->origin->tyme)
1051 {
1052 error("in migrate() wrong time found");
1053 }
1054 if (i >= proposal->old_migr_table_counter)
1055 {
1056 proposal->old_migr_table_counter += 10;
1057 proposal->migr_table =
1058 (migr_table_fmt *) myrealloc (proposal->migr_table,
1059 sizeof (migr_table_fmt) *
1060 (proposal->old_migr_table_counter));
1061 }
1062 proposal->migr_table_counter = i;
1063 return 1;
1064 }
1065
1066 int
migrateb(proposal_fmt * proposal,node * up)1067 migrateb (proposal_fmt * proposal, node * up)
1068 {
1069 long i = proposal->migr_table_counter2;
1070 migr_table_fmt *array = proposal->migr_table2;
1071 if (i > MIGRATION_LIMIT * proposal->numpop)
1072 {
1073 return 0;
1074 }
1075 if (i > 0)
1076 array[i].to = array[i - 1].from;
1077 else
1078 array[i].to = up->pop;
1079 array[i].from = migration_from (array[i].to, proposal);
1080 // printf("t: %3li %li %li\n",i,proposal->migr_table2[i].from,proposal->migr_table2[i].to);
1081 array[i].time = proposal->time;
1082 i++;
1083 if (i >= proposal->old_migr_table_counter2)
1084 {
1085 proposal->old_migr_table_counter2 += 10;
1086 proposal->migr_table2 =
1087 (migr_table_fmt *) myrealloc (proposal->migr_table2,
1088 sizeof (migr_table_fmt) *
1089 (proposal->old_migr_table_counter2));
1090 }
1091 proposal->migr_table_counter2 = i;
1092 return 1;
1093 }
1094
1095 int
migrate_old(proposal_fmt * proposal,node * up,long * old_migr_table_counter,boolean air)1096 migrate_old (proposal_fmt * proposal, node * up, long *old_migr_table_counter,
1097 boolean air)
1098 {
1099 migr_table_fmt *array;
1100 long i;
1101 if (air)
1102 {
1103 array = proposal->migr_table;
1104 i = proposal->migr_table_counter;
1105 }
1106 else
1107 {
1108 array = proposal->migr_table2;
1109 i = proposal->migr_table_counter2;
1110 }
1111 if (i > MIGRATION_LIMIT * proposal->numpop)
1112 {
1113 // FPRINTF (stdout, "migration limit reached\n");
1114 return 0;
1115 }
1116 switch (proposal->migration_model)
1117 {
1118 case ISLAND:
1119 case ISLAND_VARTHETA:
1120 case MATRIX:
1121 case MATRIX_SAMETHETA:
1122 case MATRIX_ARBITRARY:
1123 if (i > 0)
1124 array[i].to = array[i - 1].from;
1125 else
1126 array[i].to = up->pop;
1127 array[i].from = migration_from (array[i].to, proposal);
1128 // printf("i: %3li %li %li\n",i,proposal->migr_table[i].from,proposal->migr_table[i].to);
1129 // printf("b: %3li %li %li\n",i,proposal->migr_table2[i].from,proposal->migr_table2[i].to);
1130 break;
1131 case CONTINUUM:
1132 case STEPSTONE:
1133 error ("not yet implemented\n");
1134 break;
1135 default:
1136 break;
1137 }
1138 array[i++].time = proposal->time;
1139 if (i >= (*old_migr_table_counter))
1140 {
1141 (*old_migr_table_counter) += 10;
1142 if (air)
1143 {
1144 proposal->migr_table =
1145 (migr_table_fmt *) myrealloc (proposal->migr_table,
1146 sizeof (migr_table_fmt) *
1147 (*old_migr_table_counter));
1148 //xcode array = proposal->migr_table;
1149 }
1150 else
1151 {
1152 proposal->migr_table2 =
1153 (migr_table_fmt *) myrealloc (proposal->migr_table2,
1154 sizeof (migr_table_fmt) *
1155 (*old_migr_table_counter));
1156 //xcode array = proposal->migr_table;
1157 }
1158 }
1159 if (air)
1160 {
1161 proposal->migr_table_counter = i;
1162 }
1163 else
1164 {
1165 proposal->migr_table_counter2 = i;
1166 }
1167 return 1;
1168 }
1169
1170 /* migration_from() returns the FROM population when there was a migration
1171 TO population to migrate to
1172 PROPOSAL proposal-scratchpad
1173 */
1174 long
migration_from_old(long to,proposal_fmt * proposal)1175 migration_from_old (long to, proposal_fmt * proposal)
1176 {
1177 // not completely fixed for change if Nm tracking!!!! //@@
1178 long j, ii, msta, msto;
1179 MYREAL *geo = proposal->world->data->geo;
1180 MYREAL *r, rr = UNIF_RANDUM ();
1181 boolean usem = proposal->world->options->usem;
1182 r = (MYREAL *) mycalloc (1, sizeof (MYREAL) * proposal->numpop);
1183 msta = mstart (to, proposal->numpop);
1184 msto = mend (to, proposal->numpop);
1185 if (usem)
1186 {
1187 r[0] = proposal->param0[msta]/proposal->param0[to] * geo[msta];//@@
1188 for (j = 1, ii = msta + 1; ii < msto; j++, ii++)
1189 {
1190 r[j] = r[j - 1] + geo[ii] * proposal->param0[ii]/proposal->param0[to];//@@
1191 }
1192 }
1193 else
1194 {
1195 r[0] = proposal->param0[msta] * geo[msta];
1196 for (j = 1, ii = msta + 1; ii < msto; j++, ii++)
1197 {
1198 r[j] = r[j - 1] + geo[ii] * proposal->param0[ii];
1199 }
1200 }
1201 ii = 0;
1202 while (rr > r[ii] / r[j - 1])
1203 {
1204 ii++;
1205 }
1206 myfree(r);
1207 if (ii < to)
1208 return ii;
1209 else
1210 return ++ii;
1211 }
1212
1213 long
migration_from(long to,proposal_fmt * proposal)1214 migration_from (long to, proposal_fmt * proposal)
1215 {
1216 long ii = 0;
1217 MYREAL *r = proposal->world->migproblist[to];
1218 MYREAL rr = UNIF_RANDUM ();
1219 while (rr > r[ii] && ii < proposal->world->numpop-1)
1220 {
1221 ii++;
1222 }
1223 #ifdef TREEDEBUG
1224 if(r[ii] == 0.0)
1225 {
1226 warning("DEBUG: migproblist does not contain a migration probability\n");
1227 }
1228 #endif
1229 if (ii < to)
1230 return ii;
1231 else
1232 return ++ii;
1233 }
1234
1235 void
chooseTarget(proposal_fmt * proposal,timelist_fmt * timevector,node ** bordernodes,long * bordernum)1236 chooseTarget (proposal_fmt * proposal, timelist_fmt * timevector,
1237 node ** bordernodes, long *bordernum)
1238 {
1239 long actualpop = -99;
1240 node *rb = crawlback (proposal->root->next);
1241 *bordernum = 0;
1242 proposal->target = NULL;
1243 proposal->realtarget = NULL;
1244 if (proposal->migr_table_counter == 0)
1245 actualpop = proposal->origin->pop;
1246 else
1247 actualpop = proposal->migr_table[proposal->migr_table_counter - 1].from;
1248 if (rb->tyme < proposal->time)
1249 {
1250 error ("Wrong Time for action in chooseTarget()\n");
1251 }
1252 findbordernodes (rb, proposal, actualpop, &bordernodes, &proposal->listsize, bordernum,
1253 &(*timevector).tl, (*timevector).T);
1254 if (*bordernum > 0)
1255 {
1256 // found elegible linages to coalesce to
1257 proposal->target = bordernodes[RANDINT (0L, (*bordernum) - 1)];
1258 if (proposal->target != rb)
1259 {
1260 proposal->tsister = showsister (proposal->target);
1261 proposal->realtsister = crawlback (proposal->tsister)->back;
1262 }
1263 else
1264 proposal->tsister = NULL;
1265 proposal->realtarget = proposal->target;
1266 if (proposal->target->type == 'm')
1267 proposal->target = crawlback (showtop (proposal->target)->next);
1268 }
1269 else
1270 {
1271 // no lineage to coalesce to was found.
1272 proposal->target = NULL;
1273 proposal->tsister = NULL;
1274 proposal->realtsister = NULL;
1275 proposal->realtarget = NULL;
1276 }
1277 }
1278
1279 void
findbordernodes(node * theNode,proposal_fmt * proposal,long pop,node *** bordernodes,long * allocsize,long * bordernum,vtlist ** tyme,long gte)1280 findbordernodes (node * theNode, proposal_fmt * proposal, long pop,
1281 node *** bordernodes, long *allocsize, long *bordernum, vtlist ** tyme,
1282 long gte)
1283 {
1284 node *tmp, *back;
1285 // we search on the old tree and reaching the oback node that was excised from
1286 // from the residual tree we jump over it and go up the sister branch and back is
1287 // going further down than the oback node
1288 if (theNode == proposal->oback)
1289 {
1290 tmp = showtop (crawlback (proposal->osister)->back);
1291 back = showtop (proposal->oback->back);
1292 }
1293 else
1294 {
1295 tmp = showtop (theNode);
1296 back = showtop (theNode->back);
1297 }
1298 if (pop == tmp->pop && pop == back->actualpop && tmp->tyme < proposal->time
1299 && back->tyme > proposal->time)
1300 {
1301 // lineage end points bracket the the new proposed time,
1302 // therefore this branch needs to be included
1303 (*bordernodes)[(*bordernum)++] = tmp;
1304 return;
1305 }
1306 else
1307 {
1308 // the lineage does not fit the populations (it may fit the time, so)
1309 // the if here checks this
1310 if (back->tyme < proposal->time)
1311 {
1312 // this makes sure that we looked at all lineages that
1313 // may be electable to receive a coalescent have been looked
1314 // at, if we reach here the time of the proposal is older than
1315 // all available lineages on this branch and we safely can stop
1316 // following this branch.
1317 return;
1318 }
1319 // if we are still working on interior nodes and the proposal time is still
1320 // younger than then lineage bracketing nodes we need to continue our search
1321 // up in the tree following right and left branches, but need to stop when
1322 // we encounter a tip node (this is important with dated tips because these can
1323 // be interspersed in the timelist.
1324 if (tmp->type != 't')
1325 {
1326 if (tmp->next->back != NULL)
1327 findbordernodes (tmp->next->back, proposal, pop, bordernodes, allocsize,
1328 bordernum, tyme, gte);
1329 if (tmp->type != 'm' && tmp->next->next->back != NULL)
1330 findbordernodes (tmp->next->next->back, proposal, pop,
1331 bordernodes, allocsize, bordernum, tyme, gte);
1332 }
1333 }
1334 }
1335
1336 /*
1337 boolean
1338 same_pop(node * up, MYREAL tyme, long pop)
1339 {
1340 node *oldnn = showtop(up->back);
1341 node *nn = up;
1342 while (nn->tyme < tyme) {
1343 oldnn = nn;
1344 nn = showtop(nn->back);
1345 }
1346 if (oldnn->pop == pop && nn->actualpop == pop)
1347 return TRUE;
1348 else
1349 return FALSE;
1350 }
1351 */
1352
1353
1354 /* -----------------------------------------------------------------------
1355 simulates two lineages at once, if we are extending below the last node */
1356 int
pre_population(proposal_fmt * proposal,vtlist * ltime,long gte,long * slider)1357 pre_population (proposal_fmt * proposal, vtlist * ltime, long gte,
1358 long *slider)
1359 {
1360 boolean coalesced = FALSE;
1361 boolean choice = FALSE;
1362 long pop1 = -99, pop2 = -98;
1363 // long msta1=0, msto1=0;
1364 // long msta2=0, msto2=0;
1365 MYREAL ux;
1366 MYREAL rate = proposal->world->options->mu_rates[proposal->world->locus];
1367 MYREAL invrate = 1./rate;
1368 MYREAL age1, denom, rr, r0, r1, horizon, mm, mm2;
1369 // MYREAL denom_invrate;
1370 if (gte > 0)
1371 proposal->realtarget = ltime[gte - 1].eventnode;
1372 else
1373 proposal->realtarget = ltime[0].eventnode->next->back; //?????gte
1374 if (proposal->realtarget == proposal->oback)
1375 {
1376 proposal->realtarget = crawlback (proposal->osister)->back;
1377 }
1378 if (proposal->realtarget->type == 'm')
1379 {
1380 proposal->target = crawlback (proposal->realtarget->next);
1381 if (proposal->target == proposal->oback)
1382 {
1383 proposal->target = proposal->osister;
1384 }
1385 }
1386 else
1387 {
1388 proposal->target = proposal->realtarget;
1389 }
1390 proposal->tsister = NULL;
1391 pop2 = proposal->realtarget->pop;
1392 pop1 =
1393 proposal->migr_table_counter >
1394 0 ? proposal->migr_table[proposal->migr_table_counter -
1395 1].from : proposal->origin->pop;
1396 age1 =
1397 MAX (proposal->realtarget->tyme,
1398 proposal->migr_table_counter >
1399 0 ? proposal->migr_table[proposal->migr_table_counter -
1400 1].time : proposal->origin->tyme);
1401 horizon = MAX (proposal->oback->tyme, age1);
1402 while (age1 < horizon)
1403 {
1404 mm = proposal->mig0list[pop1] * invrate;
1405 //mm = 0.0;
1406 //msta1 = mstart(pop1,proposal->numpop);
1407 //msto1 = mend(pop1,proposal->numpop);
1408 // for (i = msta1; i < msto1; i++)
1409 //{
1410 //mm += proposal->param0[i];
1411 //}
1412 if (pop1 == pop2)
1413 {
1414 denom = mm + (2. / (proposal->param0[pop1] * rate));
1415 // denom_invrate = denom * invrate;
1416 ux = UNIF_RANDUM();
1417 // proposal->time = age1 - LOG(ux) / denom_invrate;
1418 proposal->time = age1 - LOG(ux) / denom;
1419 age1 = proposal->time;
1420 if (age1 < horizon)
1421 {
1422 rr = UNIF_RANDUM ();
1423 r0 = (2. / (proposal->param0[pop1] * rate)) / denom;
1424 if (rr < r0)
1425 {
1426 return 1;
1427 }
1428 }
1429 }
1430 else
1431 {
1432 denom = mm;
1433 proposal->time = age1 - LOG (UNIF_RANDUM ()) / denom;
1434 age1 = proposal->time;
1435 }
1436 if (age1 < horizon)
1437 {
1438 if (!migrate (proposal, proposal->origin))
1439 // &proposal->old_migr_table_counter, MIGRATION_AIR))
1440 {
1441 return 0;
1442 }
1443 pop1 =
1444 proposal->migr_table_counter >
1445 0 ? proposal->migr_table[proposal->migr_table_counter -
1446 1].from : proposal->origin->pop;
1447 }
1448 }
1449 age1 = horizon;
1450 while (!coalesced)
1451 {
1452 //mm = mm2 = 0;
1453 //msta1 = mstart(pop1,proposal->numpop);
1454 //msto1 = mend(pop1,proposal->numpop);
1455 //msta2 = mstart(pop2,proposal->numpop);
1456 //msto2 = mend(pop2,proposal->numpop);
1457 //for (i = msta1; i < msto1; i++)
1458 //{
1459 //mm += proposal->param0[i];
1460 //}
1461
1462 //limits treesize
1463 if(proposal->time > 1000000)
1464 return 0;
1465
1466 mm = proposal->mig0list[pop1] * invrate;
1467 mm2 = proposal->mig0list[pop2] * invrate;
1468
1469 if (pop1 == pop2)
1470 {
1471 denom = 2. * mm + (2. / proposal->param0[pop1] * rate);
1472 // denom_invrate = denom * invrate;
1473 proposal->time = age1 - LOG (UNIF_RANDUM ()) / denom;
1474 age1 = proposal->time;
1475 rr = UNIF_RANDUM ();
1476 r0 = ((2. / proposal->param0[pop1] * rate) / denom);
1477 r1 = r0 + mm / denom;
1478 if (rr < r0)
1479 {
1480 return 1;
1481 }
1482 else
1483 {
1484 if (rr < r1)
1485 {
1486 choice = TRUE;
1487 }
1488 else
1489 {
1490 choice = FALSE;
1491 }
1492 }
1493 }
1494 else
1495 { /*pop1 not equal pop2 */
1496 //for (i = msta2; i < msta2; i++)
1497 //{
1498 //mm2 += proposal->param0[i];
1499 //}
1500 denom = mm + mm2;
1501 // denom_invrate = denom * invrate;
1502 proposal->time = age1 - LOG (UNIF_RANDUM ()) / denom;
1503 age1 = proposal->time;
1504 if (RANDUM () < (mm / denom))
1505 {
1506 choice = TRUE;
1507 }
1508 else
1509 {
1510 choice = FALSE;
1511 }
1512 }
1513 if (choice)
1514 {
1515 if (!migrate (proposal, proposal->origin))
1516 // &proposal->old_migr_table_counter, MIGRATION_AIR))
1517 {
1518 return 0; /* migration limit reached */
1519 }
1520 pop1 =
1521 proposal->migr_table_counter >
1522 0 ? proposal->migr_table[proposal->migr_table_counter -
1523 1].from : proposal->origin->pop;
1524 }
1525 else
1526 {
1527 if (!migrateb (proposal, proposal->realtarget))
1528 // &proposal->old_migr_table_counter2,
1529 // MIGRATION_IN_TREE))
1530 {
1531 return 0; /* migration limit reached */
1532 }
1533 pop2 =
1534 proposal->migr_table_counter2 >
1535 0 ? proposal->migr_table2[proposal->migr_table_counter2 -
1536 1].from : proposal->realtarget->pop;
1537 }
1538 }
1539 error ("Reached the end of function without coalescing");
1540 return -1; /*makes the compiler happy */
1541 }
1542
1543 #ifdef TESTING2
1544 ///
1545 /// freeing the proposal structure, this is the replacement of the real free_proposal method
1546 /// and does not allocate, but reuses old memory
1547 //void free_proposal(proposal_fmt *proposal)/
1548 //{
1549 // do nothing
1550 //}
reset_simple_proposal_variables(proposal_fmt ** proposal)1551 void reset_simple_proposal_variables(proposal_fmt **proposal)
1552 {
1553 (*proposal)->mig_removed = FALSE;
1554 (*proposal)->rr = 0.0;
1555 (*proposal)->origin = NULL;
1556 (*proposal)->target = NULL;
1557 (*proposal)->realtarget = NULL;
1558 (*proposal)->tsister = NULL;
1559 (*proposal)->realtsister = NULL;
1560 (*proposal)->osister = NULL;
1561 (*proposal)->realosister = NULL;
1562 (*proposal)->ocousin = NULL;
1563 (*proposal)->realocousin = NULL;
1564 (*proposal)->oback = NULL;
1565 (*proposal)->realoback = NULL;
1566 //
1567 (*proposal)->connect = NULL;
1568 (*proposal)->likelihood = 0.0;
1569 (*proposal)->time = 0.0;
1570 (*proposal)->v = 0.0;
1571 (*proposal)->vs = 0.0;
1572 //
1573 #ifdef UEP
1574 (*proposal)->ueplikelihood = 0.0;
1575 #endif
1576 (*proposal)->migr_table_counter=0;
1577 (*proposal)->migr_table_counter2=0;
1578 (*proposal)->timeslice = 0;
1579 //
1580 (*proposal)->treelen = 0.0;
1581 #ifdef BEAGLE
1582 (*proposal)->parentid = 0;
1583 (*proposal)->leftid = 0;
1584 (*proposal)->rightid = 0;
1585 #endif
1586 }
1587
1588 void
reset_proposal(proposal_fmt ** proposal,world_fmt * world)1589 reset_proposal (proposal_fmt ** proposal, world_fmt *world)
1590 {
1591 const long listsize = 2 * (world->sumtips + 2);
1592 const long newsize = 4 * listsize;
1593 const long mal = (*proposal)->world->data->maxalleles[(*proposal)->world->locus];
1594 (*proposal)->likelihood = -HUGE;
1595 // pointers and values to outside structures
1596 (*proposal)->world = world;
1597 (*proposal)->datatype = world->options->datatype;
1598 (*proposal)->sumtips = world->sumtips;
1599 (*proposal)->numpop = world->numpop;
1600 (*proposal)->endsite = world->data->seq[0]->endsite;
1601 (*proposal)->fracchange = world->data->seq[0]->fracchange;
1602 (*proposal)->param0 = world->param0;
1603 (*proposal)->root = world->root;
1604 (*proposal)->migration_model = world->options->migration_model;
1605 // precalculated values
1606 (*proposal)->mig0list = world->mig0list;
1607 (*proposal)->design0list = world->design0list;
1608 (*proposal)->listsize = listsize;
1609 // ..line_t are also reset with this
1610 long i;
1611 for(i=0;i<newsize;i++)
1612 (*proposal)->nodedata[i]=NULL;
1613 // memset((*proposal)->nodedata,0,sizeof(node *) * newsize);
1614 memset((*proposal)->mf,0, sizeof(MYREAL)*(2 * (*proposal)->endsite)); // (*proposal)->mt is also freed with this
1615 // resetting pointers
1616 reset_simple_proposal_variables(proposal);
1617
1618 if (strchr (SEQUENCETYPES, (*proposal)->datatype))
1619 {
1620 zero_xseq(&(*proposal)->xf,(*proposal)->world->data->seq[0]->endsite, (*proposal)->world->options->rcategs);
1621 zero_xseq(&(*proposal)->xt,(*proposal)->world->data->seq[0]->endsite, (*proposal)->world->options->rcategs);
1622 }
1623 else
1624 {
1625 memset((*proposal)->xf.a, 0, mal);
1626 memset((*proposal)->xt.a, 0, mal);
1627 }
1628 memset((*proposal)->migr_table, 0, sizeof(migr_table_fmt) * (*proposal)->old_migr_table_counter);
1629 memset((*proposal)->migr_table2, 0, sizeof(migr_table_fmt) * (*proposal)->old_migr_table_counter2);
1630 #ifdef UEP
1631
1632 if ((*proposal)->world->options->uep)
1633 {
1634 memset((*proposal)->ueplike, 0, sizeof(MYREAL) * world->numpop * world->data->uepsites);
1635 for (j = 1; j < world->data->uepsites; ++j)
1636 (*(*proposal))->ueplike[j] = (*(*proposal))->ueplike[0] + j * world->numpop;
1637
1638 memset((*proposal)->uf.s, 0, world->data->uepsites * sizeof(pair));
1639 memset((*proposal)->ut.s, 0, world->data->uepsites * sizeof(pair));
1640 memset((*proposal)->umf, 0, world->data->uepsites * sizeof(MYREAL));
1641 memset((*proposal)->umt, 0, world->data->uepsites * sizeof(MYREAL));
1642 }
1643 #endif
1644 }
1645 #endif
1646 ///
1647 /// freeing the proposal structure, this will be replaced by a function that resets permanently
1648 /// allocated structure [reset_proposal()]
1649 void
free_masterproposal(proposal_fmt * proposal)1650 free_masterproposal (proposal_fmt * proposal)
1651 {
1652 //static long count=0;
1653 myfree(proposal->aboveorigin);
1654 // myfree(proposal->bordernodes);
1655 // myfree(proposal->line_f);
1656 // myfree(proposal->line_t);
1657 // printf("%li ",count++);fflush(stdout);
1658 myfree(proposal->nodedata); // ..line_t are also freed with this
1659 myfree(proposal->mf); // proposal->mt is also freed with this
1660 if (strchr (SEQUENCETYPES, proposal->datatype))
1661 {
1662 myfree(proposal->xf.s[0]);
1663 myfree(proposal->xt.s[0]);
1664 myfree(proposal->xf.s);
1665 myfree(proposal->xt.s);
1666 }
1667 else
1668 {
1669 myfree(proposal->xf.a);
1670 myfree(proposal->xt.a);
1671 }
1672 myfree(proposal->migr_table);
1673 myfree(proposal->migr_table2);
1674 #ifdef UEP
1675
1676 if (proposal->world->options->uep)
1677 {
1678 myfree(proposal->ueplike[0]);
1679 myfree(proposal->ueplike);
1680 myfree(proposal->uf.s);
1681 myfree(proposal->ut.s);
1682 myfree(proposal->umf);
1683 myfree(proposal->umt);
1684
1685 }
1686 #endif
1687 myfree(proposal);
1688 }
1689 //#endif /*TESTING*/
1690
1691
1692 void
free_timevector(timelist_fmt * timevector)1693 free_timevector (timelist_fmt * timevector)
1694 {
1695 // long i;
1696 // for (i = 0; i < timevector->allocT; i++)
1697 // {
1698 // myfree(timevector->tl[i].lineages);
1699 // }
1700 myfree(timevector->lineages);
1701 myfree(timevector->tl);
1702 myfree(timevector);
1703 }
1704
1705 ///
1706 /// using a global timevector to save malloc/free pair calls and
1707 /// so to save time
1708 void
free_timevector_new(timelist_fmt * timevector)1709 free_timevector_new (timelist_fmt * timevector)
1710 {
1711 // memset(timevector->lineages, 0, sizeof(long) * );
1712 // (timevector->tl);
1713 }
1714
1715 /*----------------------------------------------------------*
1716 * rejection/acceptance of the new tree according to the likelihood
1717 * and an acceptance ratio which is higher the better the
1718 * likelihood values are (-> Metropolis)
1719 */
1720 boolean
acceptlike(world_fmt * world,proposal_fmt * proposal,long g,timelist_fmt * tyme)1721 acceptlike (world_fmt * world, proposal_fmt * proposal, long g,
1722 timelist_fmt * tyme)
1723 {
1724 //proposal changes when chain is stuck, to accept some new trees without
1725 //considering data to remove out of sticky area [hopefully]
1726 static long not_accepted = 0;
1727
1728 const long limit = MIGRATION_LIMIT * world->numpop;
1729
1730 long rm = 0L;
1731 long rmc = rmigrcount (proposal);
1732
1733 MYREAL rr;
1734 MYREAL expo;
1735
1736 #ifdef UEP
1737 node *first;
1738 #endif
1739
1740 rm = proposal->migr_table_counter + proposal->migr_table_counter2
1741 + world->migration_counts - rmc;
1742
1743 if (rm > limit)
1744 {
1745 // warning might disturb parallel runs and might be not really useful at all
1746 // if (report || g != oldg)
1747 // {
1748 // warning ("migration limit (%li) exceeded: %li\n",
1749 // MIGRATION_LIMIT * world->numpop, rm);
1750 // warning
1751 // ("results may be underestimating migration rates for this chain\n");
1752 // report = FALSE;
1753 // oldg = g;
1754 // }
1755 #ifdef BEAGLE
1756 printf("Reset to: %f\n",world->likelihood[g]);
1757 //set_beagle_dirty(proposal->origin,proposal->target,showtop(world->root->next->back));
1758 //calcLnL(world, world->beagle->instance);
1759 #endif
1760 return FALSE;
1761 }
1762 #ifdef UEP
1763 if (world->options->uep)
1764 {
1765 first = first_uep2 (proposal->world->root->next->back,
1766 proposal->world->root->next->back,
1767 proposal->world->data->uepsites);
1768 proposal->firstuep = first_uep (first, proposal->world->root,
1769 proposal->world->data->uepsites);
1770 proposal->ueplikelihood = pseudo_ueplikelihood (world, proposal);
1771 proposal->likelihood = pseudotreelikelihood (world, proposal);
1772 }
1773 else
1774 {
1775 proposal->likelihood = pseudotreelikelihood (world, proposal);
1776 }
1777 #else
1778 // printf("DEBUG: world->likellihood[%li]=%f heat=%f\n",g,world->likelihood[g],world->heat);
1779 proposal->likelihood = pseudotreelikelihood (world, proposal);
1780 //printf("DEBUG: proposal->likellihood=%f heat=%f\n",proposal->likelihood, world->heat);
1781 #endif
1782
1783 //if(proposal->likelihood <= -HUGE && world->likelihood[g] <= -HUGE)
1784 // {
1785 // if(!world->in_burnin)
1786 // {
1787 //long burnin = world->options->burn_in;
1788 //world->options->burn_in *= 10;
1789 //fprintf(stdout,"likelihood is -inf, insert a burn-in phase\n");
1790 // burnin_chain(world);
1791 //world->options->burn_in = burnin;
1792 // }
1793 // return TRUE;
1794 // }
1795 //@@@@@@@@@@@@@@@@@@@@@ <<<<<<
1796 // if(world->cold)
1797 // {
1798 // printf("%f %f %f\n",world->likelihood[g], proposal->likelihood, -world->likelihood[g]+ proposal->likelihood);
1799 // }
1800 if (world->likelihood[g] < proposal->likelihood)
1801 {
1802 not_accepted = 0;
1803 return TRUE;
1804 }
1805 if (!world->options->heating)
1806 {
1807 expo = proposal->likelihood - world->likelihood[g];
1808 rr = LOG (RANDUM ());
1809 if (rr < expo)
1810 {
1811 not_accepted = 0;
1812 return TRUE;
1813 }
1814 }
1815 else
1816 {
1817 expo = (proposal->likelihood - world->likelihood[g]) * world->heat;
1818 rr = LOG (RANDUM ());
1819 if (rr < expo)
1820 {
1821 not_accepted = 0;
1822 return TRUE;
1823 }
1824 }
1825 if(world->options->prioralone)
1826 {
1827 return TRUE;
1828 }
1829 return FALSE;
1830 }
1831
1832 long
rmigrcount(proposal_fmt * proposal)1833 rmigrcount (proposal_fmt * proposal)
1834 {
1835 node *p;
1836 long count = 0;
1837 for (p = proposal->origin; p != proposal->oback; p = showtop (p->back))
1838 {
1839 if (p->type == 'm')
1840 count++;
1841 }
1842 return count;
1843 }
1844
1845
1846 ///
1847 /// generates the time interval to the next event (migration or coalescence)
1848 /// and also decides what event happens.
1849 /// This funcion is modified by the rate of the mutation rate that can be
1850 /// estimated in a Bayesian context, and be fixed in a ML context
1851 /// this rate is only influencing the time and not the eventtype
1852 /// because the rate cancels out of the ratio that chooses the event.
1853 MYREAL
eventtime(proposal_fmt * proposal,long pop,vtlist * tentry,char * event)1854 eventtime (proposal_fmt * proposal, long pop, vtlist * tentry, char *event)
1855 {
1856 // static boolean mig_force = TRUE;
1857 MYREAL interval;
1858 MYREAL lines;
1859 MYREAL denom;
1860 MYREAL invdenom;
1861 MYREAL rate = proposal->world->options->mu_rates[proposal->world->locus];
1862 MYREAL invrate = 1./ rate;
1863 MYREAL mm = proposal->mig0list[pop] * invrate ;
1864 //long design0 = proposal->design0list[pop];
1865 //could be used to calculate p(gn|go): MYREAL nu;
1866 lines = 2.0 * tentry->lineages[pop];
1867 denom = mm + (lines * (1.0 / (proposal->param0[pop]*rate)));
1868 invdenom = 1.0 / denom;
1869 interval = (-(LOG (/*nu = */UNIF_RANDUM ())) * invdenom) ;
1870 //could be used to calculate p(gn|go): proposal->nu += nu; //adds the probabilities (=random numbers);
1871 //[overcautious]
1872 //if (interval < 0.0)
1873 // error("interval time is negative");
1874 if (lines > 0.0)
1875 {
1876 if ((UNIF_RANDUM ()) < (mm * invdenom))
1877 {
1878 *event = 'm';
1879 return interval;
1880 }
1881 else
1882 {
1883 *event = 'c';
1884 return interval;
1885 }
1886 }
1887 else
1888 {
1889 // printf("mcmc1.c 1653 pop = %li, lines = %li\n", pop, tentry->lineages[pop]);
1890 *event = 'm';
1891 return interval;
1892 }
1893 }
1894
1895 /*--------------------------------------------------------*
1896 * showsister()
1897 * find the sisternode, by going down the branch and up on
1898 * the other side again, neglecting the migration nodes.
1899 */
1900 node *
showsister(node * theNode)1901 showsister (node * theNode)
1902 {
1903 node *tmp = crawlback (theNode);
1904
1905 if (tmp->next->top)
1906 {
1907 return crawlback (tmp->next->next);
1908 }
1909 else
1910 {
1911 if (tmp->next->next->top)
1912 {
1913 return crawlback (tmp->next);
1914 }
1915 else
1916 {
1917 error ("error in treestructure, cannot find sisternode\n");
1918 }
1919 }
1920 return NULL;
1921 }
1922
1923 void
count_migrations(node * p,long * count)1924 count_migrations (node * p, long *count)
1925 {
1926 if (p->type != 't')
1927 {
1928 if (p->type == 'm')
1929 {
1930 *count += 1;
1931 count_migrations (p->next->back, count);
1932 }
1933 else
1934 {
1935 count_migrations (p->next->back, count);
1936 count_migrations (p->next->next->back, count);
1937 }
1938 }
1939 }
1940
1941 MYREAL
prob_tree(world_fmt * world,timelist_fmt * tyme)1942 prob_tree (world_fmt * world, timelist_fmt * tyme)
1943 {
1944 long j, pop;
1945 MYREAL mm, cc, ss = 0;
1946 long tlfrom;
1947 boolean usem = world->options->usem;
1948 MYREAL pk;
1949 const MYREAL mu_rate = world->options->mu_rates[world->locus];
1950 //const MYREAL lmu_rate = world->options->lmu_rates[world->locus];
1951 tyme->tl[0].interval = tyme->tl[0].age;
1952 for (j = 1; j < tyme->T; j++)
1953 {
1954 tyme->tl[j].interval = tyme->tl[j].age - tyme->tl[j - 1].age;
1955 }
1956 for (j = 0; j < tyme->T - 1; j++)
1957 {
1958 mm = cc = 0.0;
1959 for (pop = 0; pop < world->numpop; pop++)
1960 {
1961 mm += tyme->tl[j].lineages[pop] * world->mig0list[pop];
1962 cc += tyme->tl[j].lineages[pop] * (tyme->tl[j].lineages[pop] -
1963 1) / (world->param0[pop] * mu_rate);
1964 }
1965 ss += -(tyme->tl[j].interval) * (mm + cc);
1966 tlfrom = tyme->tl[j].from;
1967 if (tlfrom == tyme->tl[j].to)
1968 ss += LOG2 - LOG (world->param0[tlfrom]*mu_rate);
1969 else
1970 { //@@
1971 pk = usem ? world->param0[tlfrom]/mu_rate : world->param0[tlfrom] /(mu_rate * world->param0[pop]);
1972 ss += LOG (pk);
1973 }
1974 }
1975 return ss;
1976 }
1977