1 /*------------------------------------------------------
2 Maximum likelihood estimation
3 of migration rate and effectice population size
4 using a Metropolis-Hastings Monte Carlo algorithm
5 -------------------------------------------------------
6 M C M C 2 R O U T I N E S
7
8 Tree changing routines
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-2004 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: mcmc2.c 2158 2013-04-29 01:56:20Z beerli $
37 -------------------------------------------------------*/
38 /* \file mcmc2.c
39 Contains low-level routines to manipulate the treechanger
40 */
41
42 #include "migration.h"
43 #include "sighandler.h"
44 #include "tree.h"
45 #include "uep.h"
46 #include "tools.h"
47 #include "mcmc.h"
48 #include "migrate_mpi.h"
49 #ifdef BEAGLE
50 #include "calculator.h"
51 #endif
52 #ifdef DMALLOC_FUNC_CHECK
53 #include <dmalloc.h>
54 #endif
55
56 #define comment(a) FPRINTF(stderr,a)
57
58
59 /* prototypes ------------------------------------------- */
60 /* changes the tree to the proposed tree by rebuilding parts
61 of the tree and also inserts new migration events. */
62 void coalesce1p (proposal_fmt * proposal);
63 /* and its pretend version */
64 void pretendcoalesce1p (proposal_fmt * proposal);
65
66 /* private functions */
67 /* family of tree manipulation routines: dependent on the position of
68 target and origin in the tree different routines are used.
69 oback: there is no change in the tree topology.
70 ocousin: target is the cousin-node, the branch is inserted in its
71 former sister branch.
72 rbcoa: target is the bottommost node, the ripped branch will be
73 bring the new bottommost node.
74 stancoa: all other cases. */
75 void coat_oback (proposal_fmt * proposal);
76 void coat_ocousin (proposal_fmt * proposal);
77 void coat_obbcoa (proposal_fmt * proposal);
78 void coat_rbcoa (proposal_fmt * proposal);
79 void coat_stancoa (proposal_fmt * proposal);
80 /* pretend versions of the above */
81 void target_oback (proposal_fmt * proposal);
82 void target_obbcoa (proposal_fmt * proposal);
83 void target_ocousin (proposal_fmt * proposal);
84 void target_rbcoa (proposal_fmt * proposal);
85 void target_stancoa (proposal_fmt * proposal);
86 /* subroutines in target_stancoa */
87 void t_s_upper (proposal_fmt * proposal, node * connect, node * obb);
88 node *t_s_obbncon (proposal_fmt * proposal, node * obb);
89 void t_s_tncon (proposal_fmt * proposal, node * obb);
90 void t_s_tcon (proposal_fmt * proposal, node * uppern);
91 long erase_migr_nodes (world_fmt *world, node * up);
92 node *findcrossing (node ** ptrl1, node ** ptrl2);
93 void connectnodes (world_fmt *world, node * mother, node * brother, node * sister);
94 void gotoroot (node * origin, node ** ptrlist);
95 void adjust (node * theNode, MYREAL tyme, long level);
96 void localevaluate (node * mother);
97 void copy_x (proposal_fmt * proposal, xarray_fmt xx1, xarray_fmt xx2);
98 void fix_root_pop (world_fmt *world, node * p);
99 void pseudoevaluate (proposal_fmt * proposal, xarray_fmt x, MYREAL *lx,
100 node * mother, node * newdaughter, MYREAL v);
101 node *crawl_down (node * theNode, MYREAL tyme);
102
103
104 /*=========================================================*/
105
106 void
coalesce1p(proposal_fmt * proposal)107 coalesce1p (proposal_fmt * proposal)
108 {
109 #ifdef MESS
110 int skip=1;
111 #endif
112 //static long count=0;
113 node *obb = NULL;
114 erase_migr_nodes (proposal->world, proposal->origin);
115 if (proposal->migr_table_counter2 > 0) //|| proposal->mig_removed)
116 erase_migr_nodes (proposal->world, proposal->realtarget);
117 if (proposal->target == proposal->ocousin)
118 {
119 coat_ocousin (proposal);
120 if (proposal->oback->tyme > showtop (proposal->oback->back)->tyme
121 && crawlback (proposal->oback)->type != 'r')
122 error ("problem with time encountered"); /* was >= ** */
123 }
124 else
125 {
126 if (((proposal->target == proposal->oback)
127 || (proposal->target == proposal->osister)
128 || (proposal->target == proposal->origin)))
129 {
130 //xcode obb = proposal->oback->back;
131 coat_oback (proposal);
132 if (proposal->oback->tyme > showtop (proposal->oback->back)->tyme)
133 { /* was >= ** */
134 if (proposal->oback->back->type == 'r')
135 {
136 adjust (proposal->root, proposal->oback->tyme + 10000, 0L);
137 comment ("Root time adjusted: was to small");
138 }
139 else
140 error ("problem with time encountered");
141 }
142 }
143 else
144 {
145 if (crawlback (proposal->target)->type == 'r')
146 {
147 coat_rbcoa (proposal);
148 if (proposal->oback->tyme >
149 showtop (proposal->oback->back)->tyme)
150 { /* was >= ** */
151 if (proposal->oback->back->type == 'r')
152 {
153 adjust (proposal->root, proposal->oback->tyme + 10000,
154 0L);
155 comment ("Root time adjusted: was to small");
156 }
157 else
158 error ("problem with time encountered");
159 }
160 }
161 else
162 {
163 obb = showtop (crawlback (proposal->oback));
164 if ((obb != NULL) && (obb == proposal->target))
165 {
166 coat_obbcoa (proposal);
167 if (proposal->oback->tyme > showtop (proposal->oback->back)->tyme) /* was >= ** */
168 error ("problem with time encountered");
169 }
170 else
171 {
172 coat_stancoa (proposal);
173 if (proposal->oback->tyme > showtop (proposal->oback->back)->tyme) /* was >= ** */
174 error ("problem with time encountered");
175 }
176 }
177 }
178 }
179 set_pop (proposal->oback, proposal->origin->actualpop,
180 proposal->origin->actualpop);
181 // free_migr_node(proposal->origin, crawlback(proposal->origin));
182 // if(++count % 10000 == 0)
183 // {
184 // printf("%li heat %f> %li %li %li %li\n", count, proposal->world->heat, proposal->migr_table_counter, proposal->migr_table_counter2, proposal->world->node_collection_count, proposal->world->node_collection_allocated);
185 //
186 // }
187
188 if (proposal->migr_table_counter > 0)
189 {
190 set_pop (proposal->oback,
191 proposal->migr_table[proposal->migr_table_counter - 1].from,
192 proposal->migr_table[proposal->migr_table_counter - 1].from);
193 insert_migr_node (proposal->world, proposal->origin,
194 crawlback (proposal->origin), proposal->migr_table,
195 &(proposal->migr_table_counter));
196 }
197 proposal->migr_table_counter = 0;
198 if (proposal->migr_table_counter2 > 0)
199 {
200 insert_migr_node (proposal->world, proposal->realtarget,
201 crawlback (proposal->target), proposal->migr_table2,
202 &(proposal->migr_table_counter2));
203 proposal->migr_table_counter2 = 0;
204 }
205 /* if (proposal->oback == crawlback(proposal->root->next)) */
206 fix_root_pop (proposal->world, crawlback (proposal->root->next));
207 // if(proposal->world->in_last_chain)
208 traverse_tagnew(crawlback (proposal->root->next),proposal->origin);
209 #ifdef MESS
210 //traverse_checker(proposal->root);
211 //if(skip==0)
212 // {
213 FPRINTF (stdout,
214 "\n[& Comment: Locus %li, best log likelihood = %f (to file),heat=%f]\n",
215 proposal->world->locus + 1, proposal->world->likelihood[proposal->world->G], proposal->world->heat);
216
217 debugtreeout (stdout, crawlback (proposal->world->root->next),
218 crawlback (proposal->world->root->next), 0);
219 // }
220 #endif
221 }
222
223
224
225 void
pretendcoalesce1p(proposal_fmt * proposal)226 pretendcoalesce1p (proposal_fmt * proposal)
227 {
228 if (proposal->target == proposal->ocousin)
229 {
230 target_ocousin (proposal);
231 }
232 else
233 {
234 if (((proposal->target == proposal->oback)
235 || (proposal->target == proposal->osister)
236 || (proposal->target == proposal->origin)))
237 {
238 target_oback (proposal);
239 }
240 else
241 {
242 if (crawlback (proposal->target)->type == 'r')
243 {
244 target_rbcoa (proposal);
245 }
246 else
247 {
248 if ((showtop (crawlback (proposal->oback)) != NULL)
249 && (showtop (crawlback (proposal->oback)) ==
250 proposal->target))
251 {
252 target_obbcoa (proposal);
253 }
254 else
255 {
256 target_stancoa (proposal);
257 }
258 }
259 }
260 }
261 }
262
263 /*================================================================*/
264
265 void
coat_oback(proposal_fmt * proposal)266 coat_oback (proposal_fmt * proposal)
267 {
268 node *obb, *obm;
269 #ifdef TREEDEBUG
270 printf("coat_oback\n");
271 #endif
272 obb = showtop (crawlback (proposal->oback));
273 connectnodes (proposal->world,showtop (proposal->oback->back), proposal->osister, NULL);
274 adjust_time (proposal->oback, proposal->time);
275 obm = showtop (crawl_down (proposal->osister, proposal->time)->back);
276 connectnodes (proposal->world, proposal->oback, proposal->origin, proposal->osister);
277 proposal->oback->back = NULL;
278 if (obm->type == 'm')
279 {
280 connectnodes (proposal->world,obm, proposal->oback, NULL);
281 }
282
283 if (obb->type == 'r')
284 {
285 connectnodes (proposal->world,obb, proposal->oback, NULL);
286 erase_migr_nodes(proposal->world,proposal->oback);
287 }
288 else
289 {
290 connectnodes (proposal->world,obb, proposal->oback, proposal->ocousin);
291 }
292 adjust (obb, obb->tyme, 2L);
293 set_dirty (proposal->oback);
294 localevaluate (proposal->oback);
295 }
296
297 void
coat_ocousin(proposal_fmt * proposal)298 coat_ocousin (proposal_fmt * proposal)
299 {
300 node *obb, *rcback;
301 #ifdef TREEDEBUG
302 printf("coat_ocousin\n");
303 #endif
304
305 obb = showtop (crawlback (proposal->oback));
306 connectnodes (proposal->world,showtop (proposal->oback->back), proposal->osister, NULL);
307 rcback = showtop (crawl_down (proposal->ocousin, proposal->time)->back);
308 adjust_time (proposal->oback, proposal->time);
309 proposal->oback->back = NULL;
310 connectnodes (proposal->world,rcback, proposal->oback, NULL);
311 connectnodes (proposal->world,proposal->oback, proposal->ocousin, proposal->origin);
312 connectnodes (proposal->world,obb, proposal->oback, proposal->osister);
313 adjust (obb, obb->tyme, 2L);
314 set_dirty (proposal->oback);
315 set_dirty (obb);
316 if (crawlback (obb)->type != 'r')
317 localevaluate (showtop (crawlback (obb)));
318 }
319
320 void
coat_rbcoa(proposal_fmt * proposal)321 coat_rbcoa (proposal_fmt * proposal)
322 {
323 node *obb, *root;
324 #ifdef TREEDEBUG
325 printf("coat_rbcoa\n");
326 #endif
327
328 if (proposal->ocousin != NULL)
329 {
330 obb = showtop (crawlback (proposal->oback));
331 root = proposal->root;
332 connectnodes (proposal->world,showtop (proposal->oback->back), proposal->osister, NULL);
333 connectnodes (proposal->world,obb, proposal->ocousin, proposal->osister);
334 adjust (obb, obb->tyme, 2L);
335 adjust_time (proposal->oback, proposal->time);
336 proposal->oback->back = NULL;
337 connectnodes (proposal->world,proposal->oback, proposal->origin, proposal->target);
338 connectnodes (proposal->world,root, proposal->oback, NULL);
339 adjust (proposal->oback, proposal->time, 2L);
340 set_dirty (obb);
341 set_dirty (proposal->oback);
342 localevaluate (showtop (crawlback (obb)));
343 }
344 else
345 {
346 error ("error in coat_rbcoa\n");
347 }
348 }
349
350 void
coat_obbcoa(proposal_fmt * proposal)351 coat_obbcoa (proposal_fmt * proposal)
352 {
353 node /**obb,*/ * proposalack;
354 #ifdef TREEDEBUG
355 printf("coat_obbcoa\n");
356 #endif
357
358 if (proposal->ocousin != NULL)
359 {
360 /* obb = showtop(crawlback(proposal->oback)); */
361 connectnodes (proposal->world,showtop (proposal->oback->back), proposal->osister, NULL);
362 // proposalack = showtop (crawl_down (proposal->target, proposal->time)->back);
363 proposalack = showtop (crawlback (proposal->target));
364 adjust_time (proposal->oback, proposal->time);
365 connectnodes (proposal->world,proposal->target, proposal->ocousin, proposal->osister);
366 adjust (proposal->target, proposal->target->tyme, 1L);
367 proposal->oback->back = NULL;
368 adjust_time (proposal->oback, proposal->time);
369 if (showtop (proposal->realtarget->back)->type == 'm')
370 {
371 connectnodes (proposal->world,showtop (proposal->realtarget->back), proposal->oback,
372 NULL);
373 }
374 connectnodes (proposal->world,proposal->oback, proposal->origin, proposal->target);
375 connectnodes (proposal->world,proposalack, proposal->tsister, proposal->oback);
376 adjust (proposalack, proposalack->tyme, 2L);
377 set_dirty (proposal->target);
378 set_dirty (proposal->oback);
379 localevaluate (proposalack);
380 }
381 else
382 {
383 error ("error in coat_obbcoa\n");
384 }
385 }
386
387 void
coat_stancoa(proposal_fmt * proposal)388 coat_stancoa (proposal_fmt * proposal)
389 {
390 node *obb, *proposalack;
391 #ifdef TREEDEBUG
392 printf("coat_stancoa\n");
393 #endif
394
395 if (proposal->ocousin != NULL)
396 {
397 obb = showtop (crawlback (proposal->oback));
398 //proposalack = showtop (crawl_down (proposal->target, proposal->time)->back);
399 proposalack = showtop (crawlback (proposal->target));
400 connectnodes (proposal->world,showtop (proposal->oback->back), proposal->osister, NULL);
401 proposal->oback->back = NULL;
402 connectnodes (proposal->world,obb, proposal->osister, proposal->ocousin);
403 adjust_time (proposal->oback, proposal->time);
404 if (showtop (proposal->realtarget->back)->type == 'm')
405 {
406 connectnodes (proposal->world,showtop (proposal->realtarget->back), proposal->oback,
407 NULL);
408 }
409 connectnodes (proposal->world,proposal->oback, proposal->origin, proposal->target);
410 connectnodes (proposal->world,proposalack, proposal->oback, proposal->tsister);
411 adjust (proposalack, proposalack->tyme, 3L);
412 adjust (obb, obb->tyme, 3L);
413 set_dirty (obb);
414 set_dirty (proposal->oback);
415 localevaluate (proposalack);
416 localevaluate (obb);
417 }
418 else
419 {
420 adjust_time (proposal->oback, proposal->time);
421 obb = showtop (crawlback (proposal->oback));
422 connectnodes (proposal->world,showtop (proposal->oback->back), proposal->osister, NULL);
423 //proposalack = showtop (crawl_down (proposal->target, proposal->time)->back);
424 proposalack = showtop (crawlback (proposal->target));
425 adjust_time (proposal->oback, proposal->time);
426 proposal->oback->back = NULL;
427 if (showtop (proposal->realtarget->back)->type == 'm')
428 connectnodes (proposal->world,showtop (proposal->realtarget->back), proposal->oback,
429 NULL);
430 connectnodes (proposal->world,proposal->oback, proposal->origin, proposal->target);
431 connectnodes (proposal->world,proposalack, proposal->oback, proposal->tsister);
432 /* either proposalack is osister or a descendent of her
433 other case are disallowed by time constraints */
434 connectnodes (proposal->world,obb, proposal->osister, NULL);
435 adjust (proposalack, proposalack->tyme, 3L);
436 adjust (obb, obb->tyme, 3L);
437 set_dirty (proposal->oback);
438 localevaluate (proposalack);
439
440 }
441 }
442
443
444 void
target_oback(proposal_fmt * proposal)445 target_oback (proposal_fmt * proposal)
446 {
447 node *obb = showtop (crawlback (proposal->oback));
448 #ifdef TREEDEBUG
449 printf("target_oback\n");
450 #endif
451 copy_x (proposal, proposal->xf, proposal->origin->x);
452 #ifdef UEP
453
454 if(proposal->world->options->uep)
455 copy_uepx(proposal,proposal->uf, proposal->origin->ux);
456 #endif
457
458 memcpy (proposal->mf, proposal->origin->scale,
459 sizeof (MYREAL) * proposal->endsite);
460 proposal->v = (proposal->time - proposal->origin->tyme);
461 proposal->vs = (proposal->time - proposal->osister->tyme);
462 #ifdef BEAGLE
463 long id0, id1, id2;
464 long bid1, bid2;
465 double v1, v2;
466 id0 = proposal->oback->id;
467 id1 = proposal->origin->id;
468 bid1 = proposal->origin->bid;
469 v1 = proposal->v;
470 id2 = proposal->osister->id;
471 bid2 = proposal->osister->bid;
472 v2 = proposal->vs;
473 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
474 evaluate_beagle_instances_proposal (proposal,obb, proposal->oback, id0, proposal->oback->bid,
475 fabs (proposal->time - obb->tyme));
476 #else
477 pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
478 proposal->osister->x, proposal->osister->scale, proposal->vs);
479 pseudoevaluate (proposal, proposal->xf, proposal->mf, obb, proposal->oback,
480 fabs (proposal->time - obb->tyme));
481 #endif
482 }
483
484 void
target_ocousin(proposal_fmt * proposal)485 target_ocousin (proposal_fmt * proposal)
486 {
487 node *obb;
488 #ifdef TREEDEBUG
489 printf("target_ocousin\n");
490 #endif
491
492 obb = showtop (crawlback (proposal->oback));
493 copy_x (proposal, proposal->xf, proposal->origin->x);
494 #ifdef UEP
495
496 if(proposal->world->options->uep)
497 copy_uepx(proposal,proposal->uf, proposal->origin->ux);
498 #endif
499
500 memcpy (proposal->mf, proposal->origin->scale,
501 sizeof (MYREAL) * proposal->endsite);
502 copy_x (proposal, proposal->xt, proposal->ocousin->x);
503 #ifdef UEP
504
505 if(proposal->world->options->uep)
506 copy_uepx(proposal,proposal->ut, proposal->ocousin->ux);
507 #endif
508
509 memcpy (proposal->mt, proposal->ocousin->scale,
510 sizeof (MYREAL) * proposal->endsite);
511 proposal->v = (proposal->time - proposal->origin->tyme);
512 proposal->vs = fabs (proposal->ocousin->tyme - proposal->time);
513 #ifdef BEAGLE
514 long id0, id1, id2;
515 long bid1, bid2;
516 double v1, v2;
517 id0 = proposal->oback->id;
518 id1 = proposal->origin->id;
519 bid1 = proposal->origin->bid;
520 v1 = proposal->v;
521 id2 = proposal->ocousin->id;
522 bid2 = proposal->ocousin->bid;
523 v2 = proposal->vs;
524 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
525 #else
526 pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
527 proposal->xt, proposal->mt, proposal->vs);
528 #endif
529 proposal->v = fabs (proposal->time - obb->tyme);
530 proposal->vs = obb->tyme - proposal->osister->tyme;
531 #ifdef BEAGLE
532 id0 = obb->id;
533 id1 = new_id(proposal->oback->id,proposal->world->sumtips);
534 bid1 = proposal->oback->bid;
535 v1 = proposal->v;
536 id2 = proposal->osister->id;
537 bid2 = proposal->osister->bid;
538 v2 = proposal->vs;
539 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
540 evaluate_beagle_instances_proposal (proposal,
541 showtop (crawlback (obb)), obb, obb->id, obb->bid,
542 (showtop (crawlback (obb))->tyme - obb->tyme));
543
544 #else
545 pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
546 proposal->osister->x, proposal->osister->scale, proposal->vs);
547 pseudoevaluate (proposal, proposal->xf, proposal->mf,
548 showtop (crawlback (obb)), obb,
549 (showtop (crawlback (obb))->tyme - obb->tyme));
550 #endif
551 }
552
553 void
target_rbcoa(proposal_fmt * proposal)554 target_rbcoa (proposal_fmt * proposal)
555 {
556 node *obb;
557 #ifdef TREEDEBUG
558 printf("target_rbcoa\n");
559 #endif
560 if (proposal->ocousin != NULL)
561 {
562 obb = showtop (crawlback (proposal->oback));
563 proposal->vs = obb->tyme - proposal->osister->tyme;
564 proposal->v = obb->tyme - proposal->ocousin->tyme;
565 copy_x (proposal, proposal->xt, proposal->osister->x);
566 #ifdef UEP
567
568 if(proposal->world->options->uep)
569 copy_uepx(proposal,proposal->ut, proposal->osister->ux);
570 #endif
571
572 memcpy (proposal->mt, proposal->osister->scale,
573 sizeof (MYREAL) * proposal->endsite);
574 #ifdef BEAGLE
575 long id0, id1, id2;
576 long bid1, bid2;
577 double v1, v2;
578 id0 = obb->id;
579 id1 = proposal->osister->id;
580 bid1 = proposal->osister->bid;
581 v1 = proposal->vs;
582 id2 = proposal->ocousin->id;
583 bid2 = proposal->ocousin->bid;
584 v2 = proposal->v;
585 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
586 evaluate_beagle_instances_proposal(proposal,
587 showtop (crawlback (obb)), obb, obb->id, obb->bid,
588 showtop (crawlback (obb))->tyme - obb->tyme);
589 #else
590 pseudonuview (proposal, proposal->xt, proposal->mt, proposal->vs,
591 proposal->ocousin->x, proposal->ocousin->scale,
592 proposal->v);
593 pseudoevaluate (proposal, proposal->xt, proposal->mt,
594 showtop (crawlback (obb)), obb,
595 showtop (crawlback (obb))->tyme - obb->tyme);
596 #endif
597 proposal->v = proposal->time - proposal->origin->tyme;
598 proposal->vs = proposal->time - proposal->target->tyme;
599 copy_x (proposal, proposal->xf, proposal->origin->x);
600 #ifdef UEP
601
602 if(proposal->world->options->uep)
603 copy_uepx(proposal,proposal->uf, proposal->origin->ux);
604 #endif
605
606 memcpy (proposal->mf, proposal->origin->scale,
607 sizeof (MYREAL) * proposal->endsite);
608 #ifdef BEAGLE
609 id0 = proposal->oback->id;
610 id1 = proposal->origin->id;
611 bid1 = proposal->origin->bid;
612 v1 = proposal->v;
613 id2 = new_id(proposal->target->id, proposal->world->sumtips);//target == root->next->back
614 bid2 = proposal->target->bid;
615 v2 = proposal->vs;
616 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
617 #else
618 pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
619 proposal->xt, proposal->mt, proposal->vs);
620 #endif
621 }
622 else
623 {
624 error ("error in target_rbcoa\n");
625 }
626 }
627
628 void
target_obbcoa(proposal_fmt * proposal)629 target_obbcoa (proposal_fmt * proposal)
630 {
631 node *obb, *obbb;
632 #ifdef TREEDEBUG
633 printf("target_obbcoa\n");
634 #endif
635
636 if (proposal->ocousin != NULL)
637 {
638 obb = showtop (crawlback (proposal->oback));
639 proposal->vs = obb->tyme - proposal->osister->tyme;
640 proposal->v = obb->tyme - proposal->ocousin->tyme;
641 copy_x (proposal, proposal->xt, proposal->osister->x);
642 #ifdef UEP
643
644 if(proposal->world->options->uep)
645 copy_uepx(proposal,proposal->ut, proposal->osister->ux);
646 #endif
647
648 memcpy (proposal->mt, proposal->osister->scale,
649 sizeof (MYREAL) * proposal->endsite);
650 #ifdef BEAGLE
651 long id0, id1, id2;
652 long bid1, bid2;
653 double v1, v2;
654 id0 = obb->id;
655 id1 = proposal->osister->id;
656 bid1 = proposal->osister->bid;
657 v1 = proposal->vs;
658 id2 = proposal->ocousin->id;
659 bid2 = proposal->ocousin->bid;
660 v2 = proposal->v;
661 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
662 #else
663 pseudonuview (proposal, proposal->xt, proposal->mt, proposal->vs,
664 proposal->ocousin->x, proposal->ocousin->scale,
665 proposal->v);
666 #endif
667 proposal->v = proposal->time - proposal->origin->tyme;
668 proposal->vs = fabs (obb->tyme - proposal->time);
669 copy_x (proposal, proposal->xf, proposal->origin->x);
670 #ifdef UEP
671
672 if(proposal->world->options->uep)
673 copy_uepx(proposal,proposal->uf, proposal->origin->ux);
674 #endif
675
676 memcpy (proposal->mf, proposal->origin->scale,
677 sizeof (MYREAL) * proposal->endsite);
678 #ifdef BEAGLE
679 id0 = proposal->oback->id;
680 id1 = proposal->origin->id;
681 bid1 = proposal->origin->bid;
682 v1 = proposal->v;
683 id2 = new_id(obb->id, proposal->world->sumtips);
684 bid2 = obb->bid;
685 v2 = proposal->vs;
686 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
687 proposal->v = fabs (showtop (crawlback (obb))->tyme - proposal->time);
688 obbb = showtop (crawlback (obb));
689 /*check this carefully: fakes obb and uses id of oback which is inserted between obb and obbb*/
690 /*this should guarantee the proper evaluation but feeding the right ids into beagle*/
691 evaluate_beagle_instances_proposal (proposal,obbb, obb, proposal->oback->id, proposal->oback->bid,
692 proposal->v);
693 #else
694 pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
695 proposal->xt, proposal->mt, proposal->vs);
696 proposal->v = fabs (showtop (crawlback (obb))->tyme - proposal->time);
697 obbb = showtop (crawlback (obb));
698 pseudoevaluate (proposal, proposal->xf, proposal->mf, obbb, /* is this correct??????*/ obb,
699 proposal->v);
700 #endif
701 }
702 else
703 {
704 error ("error in target_obbcoa\n");
705 }
706 }
707
708 void
target_stancoa(proposal_fmt * proposal)709 target_stancoa (proposal_fmt * proposal)
710 {
711 node *obb, *nn, *oldnn, *d1, *d2, *proposalack, *uppern = NULL; /* o R */
712 node **double_line;
713 #ifdef TREEDEBUG
714 printf("target_stancoa\n");
715 #endif
716 double_line = (node **) mycalloc (2, sizeof (node *));
717 if (proposal->ocousin != NULL)
718 {
719 obb = showtop (crawlback (proposal->oback));
720 gotoroot (proposal->target, proposal->line_t);
721 double_line[0] = proposal->osister;
722 /*findcrossing needs a last NULL element */
723 if (findcrossing (double_line, proposal->line_t) != NULL)
724 {
725 t_s_upper (proposal, proposal->osister, obb);
726 }
727 else
728 {
729 double_line[0] = proposal->ocousin;
730 /*findcrossing needs a last NULL element */
731 if (findcrossing (double_line, proposal->line_t) != NULL)
732 {
733 t_s_upper (proposal, proposal->ocousin, obb);
734 }
735 else
736 {
737 gotoroot (obb, proposal->line_f);
738 proposal->connect =
739 findcrossing (proposal->line_f, proposal->line_t);
740 proposal->vs = obb->tyme - proposal->osister->tyme;
741 proposal->v = obb->tyme - proposal->ocousin->tyme;
742 copy_x (proposal, proposal->xt, proposal->osister->x);
743 #ifdef UEP
744
745 if(proposal->world->options->uep)
746 copy_uepx(proposal,proposal->ut, proposal->osister->ux);
747 #endif
748
749 memcpy (proposal->mt, proposal->osister->scale,
750 sizeof (MYREAL) * proposal->endsite);
751 if (obb!=proposal->connect) {
752 uppern = t_s_obbncon (proposal, obb);
753 }
754 else
755 {
756 uppern = proposal->connect;
757 }
758 proposal->v = proposal->time - proposal->origin->tyme;
759 copy_x (proposal, proposal->xf, proposal->origin->x);
760 #ifdef UEP
761
762 if(proposal->world->options->uep)
763 copy_uepx(proposal,proposal->uf, proposal->origin->ux);
764 #endif
765
766 memcpy (proposal->mf, proposal->origin->scale,
767 sizeof (MYREAL) * proposal->endsite);
768 if (proposal->target != proposal->connect)
769 {
770 // t_s_tncon (proposal, obb);
771 t_s_tncon (proposal, uppern);
772 }
773 else
774 {
775 t_s_tcon (proposal, uppern);
776 }
777 }
778 }
779 }
780 else
781 {
782 proposal->v = proposal->time - proposal->origin->tyme;
783 proposal->vs = proposal->time - proposal->target->tyme;
784 copy_x (proposal, proposal->xf, proposal->origin->x);
785 #ifdef UEP
786
787 if(proposal->world->options->uep)
788 copy_uepx(proposal,proposal->uf, proposal->origin->ux);
789 #endif
790
791 memcpy (proposal->mf, proposal->origin->scale,
792 sizeof (MYREAL) * proposal->endsite);
793 copy_x (proposal, proposal->xt, proposal->target->x);
794 #ifdef UEP
795
796 if(proposal->world->options->uep)
797 copy_uepx(proposal,proposal->ut, proposal->target->ux);
798 #endif
799
800 memcpy (proposal->mt, proposal->target->scale,
801 sizeof (MYREAL) * proposal->endsite);
802 //CHECK HERE memcpy(proposal->mt,proposal->osister->scale, sizeof(MYREAL)*proposal->endsite);
803 #ifdef BEAGLE
804 long id0, id1, id2;
805 long bid1, bid2;
806 double v1, v2;
807 id0 = proposal->oback->id;
808 id1 = proposal->origin->id;
809 bid1 = proposal->origin->bid;
810 v1 = proposal->v;
811 id2 = proposal->target->id;
812 bid2 = proposal->target->bid;
813 v2 = proposal->vs;
814 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
815 #else
816 pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
817 proposal->xt, proposal->mt, proposal->vs);
818 #endif
819 proposal->v =
820 fabs (((proposalack =
821 showtop (crawlback (proposal->target)))->tyme -
822 proposal->time));
823 if (proposalack != proposal->oback)
824 {
825 children (proposalack, &d1, &d2);
826 if (d1 == proposal->target)
827 {
828 d1 = d2;
829 d2 = proposal->target;
830 }
831 #ifdef BEAGLE
832 id0 = proposalack->id;
833 id1 = new_id(proposal->oback->id,proposal->world->sumtips);
834 bid1 = proposal->oback->bid;
835 v1 = proposal->v;
836 id2 = d1->id;
837 bid2 = d1->bid;
838 v2 = d1->v;
839 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
840 #else
841 pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
842 d1->x, d1->scale, d1->v);
843 #endif
844 oldnn = proposalack;
845 nn = showtop (crawlback (oldnn));
846 while (nn != proposal->oback)
847 {
848 children (nn, &d1, &d2);
849 if (d1 == oldnn)
850 {
851 d1 = d2;
852 d2 = oldnn;
853 }
854 #ifdef BEAGLE
855 id0 = nn->id;
856 id1 = new_id(oldnn->id,proposal->world->sumtips);
857 bid1 = oldnn->bid;
858 v1 = oldnn->v;
859 id2 = d1->id;
860 bid2 = d1->bid;
861 v2 = d1->v;
862 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
863 #else
864 pseudonuview (proposal, proposal->xf, proposal->mf, oldnn->v,
865 d1->x, d1->scale, d1->v);
866 #endif
867 oldnn = nn;
868 nn = showtop (crawlback (nn));
869 }
870 }
871 }
872 myfree(double_line);
873 }
874
875 void
t_s_upper(proposal_fmt * proposal,node * connect,node * obb)876 t_s_upper (proposal_fmt * proposal, node * connect, node * obb)
877 {
878 node *nn, *oldnn, *d1, *d2;
879 proposal->v = proposal->time - proposal->origin->tyme;
880 proposal->vs = proposal->time - proposal->target->tyme;
881 copy_x (proposal, proposal->xf, proposal->origin->x);
882 #ifdef UEP
883
884 if(proposal->world->options->uep)
885 copy_uepx(proposal,proposal->uf, proposal->origin->ux);
886 #endif
887
888 memcpy (proposal->mf, proposal->origin->scale,
889 sizeof (MYREAL) * proposal->endsite);
890 #ifdef BEAGLE
891 long id0, id1, id2;
892 long bid1, bid2;
893 double v1, v2;
894 id0 = proposal->oback->id;
895 id1 = proposal->origin->id;
896 bid1 = proposal->origin->bid;
897 v1 = proposal->v;
898 id2 = proposal->target->id;
899 bid2 = proposal->target->bid;
900 v2 = proposal->vs;
901 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
902 #else
903 pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
904 proposal->target->x, proposal->target->scale, proposal->vs);
905 #endif
906 nn = showtop (crawlback (proposal->target));
907 oldnn = proposal->target;
908 proposal->vs = fabs (nn->tyme - proposal->time);
909 #ifdef BEAGLE
910 children (nn, &d1, &d2);
911 if (d1 == oldnn)
912 {
913 d1 = d2;
914 d2 = oldnn;
915 }
916 id0 = nn->id;
917 id1 = new_id(proposal->oback->id, proposal->world->sumtips);
918 bid1 = proposal->oback->bid;
919 v1 = proposal->vs;
920 id2 = d1->id;
921 bid2 = d1->bid;
922 v2 = d1->v;
923 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
924 proposal->vs = nn->v;
925 if (!(nn != connect && nn->type!='r'))
926 {
927 nn = showtop (crawlback (nn));
928 }
929 else
930 {
931 #endif
932 while (nn != connect && nn->type!='r')
933 {
934 children (nn, &d1, &d2);
935 if (d1 == oldnn)
936 {
937 d1 = d2;
938 d2 = oldnn;
939 }
940 #ifdef BEAGLE
941 id0 = nn->id;
942 id1 = oldnn->id;
943 bid1 = oldnn->bid;
944 v1 = proposal->vs;
945 id2 = d1->id;
946 bid2 = d1->bid;
947 v2 = d1->v;
948 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
949 #else
950 pseudonuview (proposal, proposal->xf, proposal->mf, proposal->vs, d1->x,
951 d1->scale, d1->v);
952 #endif
953 proposal->vs = nn->v;
954 oldnn = nn;
955 nn = showtop (crawlback (nn));
956 }
957 children (nn, &d1, &d2);
958 if (d1 == oldnn)
959 {
960 d1 = d2;
961 d2 = oldnn;
962 }
963 #ifdef BEAGLE
964 id0 = nn->id;
965 id1 = new_id(oldnn->id, proposal->world->sumtips);
966 bid1 = oldnn->bid;
967 v1 = proposal->vs;
968 id2 = d1->id;
969 bid2 = d1->bid;
970 v2 = d1->v;
971 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
972 #else
973 pseudonuview (proposal, proposal->xf, proposal->mf, proposal->vs, d1->x,
974 d1->scale, d1->v);
975 #endif
976 proposal->v = obb->tyme - proposal->osister->tyme;
977 proposal->vs = obb->tyme - proposal->ocousin->tyme;
978 if (connect == proposal->ocousin)
979 {
980 #ifdef BEAGLE
981 id0 = obb->id;
982 id1 = proposal->osister->id;
983 bid1 = proposal->osister->bid;
984 v1 = proposal->v;
985 id2 = proposal->ocousin->id;
986 bid2 = proposal->ocousin->bid;
987 v2 = proposal->vs;
988 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
989 #else
990 pseudonuview (proposal, proposal->xf, proposal->mf, proposal->vs,
991 proposal->osister->x, proposal->osister->scale,
992 proposal->v);
993 #endif
994 }
995 else
996 {
997 #ifdef BEAGLE
998 id0 = obb->id;
999 id1 = proposal->osister->id;
1000 bid1 = proposal->osister->bid;
1001 v1 = proposal->v;
1002 id2 = proposal->ocousin->id;
1003 bid2 = proposal->ocousin->bid;
1004 v2 = proposal->vs;
1005 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
1006 #else
1007 pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
1008 proposal->ocousin->x, proposal->ocousin->scale,
1009 proposal->vs);
1010 #endif
1011 }
1012 #ifdef BEAGLE
1013 }
1014 evaluate_beagle_instances_proposal (proposal, showtop (crawlback (obb)), obb, obb->id, obb->bid, obb->v);
1015 #else
1016 pseudoevaluate(proposal,proposal->xf,proposal->mf,showtop (crawlback (obb)), obb, obb->v );
1017 #endif
1018 }
1019
1020 node *
t_s_obbncon(proposal_fmt * proposal,node * obb)1021 t_s_obbncon (proposal_fmt * proposal, node * obb)
1022 {
1023
1024 node *nn, *oldnn, *d1, *d2;
1025 /* FPRINTF(stdout,"t_s_obbncon\n"); */
1026 #ifdef BEAGLE
1027 long id0, id1, id2;
1028 long bid1, bid2;
1029 double v1, v2;
1030 // printf("????????\n");
1031 id0 = obb->id;
1032 id1 = proposal->osister->id;
1033 bid1 = proposal->osister->bid;
1034 v1 = proposal->vs;
1035 id2 = proposal->ocousin->id;
1036 bid2 = proposal->ocousin->bid;
1037 v2 = proposal->v;
1038 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
1039 #else
1040 pseudonuview (proposal, proposal->xt, proposal->mt, proposal->vs,
1041 proposal->ocousin->x, proposal->ocousin->scale, proposal->v);
1042 #endif
1043 nn = showtop (crawlback (obb));
1044 oldnn = obb;
1045 proposal->vs = fabs (nn->tyme - obb->tyme);
1046 while (nn != proposal->connect)
1047 {
1048 children (nn, &d1, &d2);
1049 if (d1 == oldnn)
1050 {
1051 d1 = d2;
1052 d2 = oldnn;
1053 }
1054 #ifdef BEAGLE
1055 id0 = nn->id; // will be offset in prepare_beagle_instances_proposal()
1056 id1 = new_id(oldnn->id, proposal->world->sumtips);
1057 bid1 = oldnn->bid;
1058 v1 = proposal->vs;
1059 id2 = d1->id;
1060 bid2 = d1->bid;
1061 v2 = d1->v;
1062 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
1063 #else
1064 pseudonuview (proposal, proposal->xt, proposal->mt, proposal->vs, d1->x,
1065 d1->scale, d1->v);
1066 #endif
1067 proposal->vs = nn->v;
1068 oldnn = nn;
1069 nn = showtop (crawlback (nn));
1070 }
1071 return oldnn;
1072 }
1073
1074 void
t_s_tncon(proposal_fmt * proposal,node * obb)1075 t_s_tncon (proposal_fmt * proposal, node * obb)
1076 {
1077 node *nn, *oldnn, *d1, *d2;
1078 #ifdef BEAGLE
1079 long id0, id1, id2;
1080 long bid1, bid2;
1081 double v1, v2;
1082 id0 = proposal->oback->id;//id0 will point to the new id+nodep_boundary
1083 id1 = proposal->origin->id;
1084 bid1 = proposal->origin->bid;
1085 v1 = proposal->v;
1086 id2 = proposal->target->id;
1087 bid2 = proposal->target->bid;
1088 v2 = proposal->time - proposal->target->tyme;
1089 // connect oback, origin, target
1090 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
1091 #else
1092 pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
1093 proposal->target->x, proposal->target->scale,
1094 proposal->time - proposal->target->tyme);
1095 #endif
1096 nn = showtop (crawlback (proposal->target));
1097 oldnn = proposal->target;
1098 // calculates the time between the newly inserted oback and the back of target
1099 proposal->v = fabs (nn->tyme - proposal->time);
1100 #ifdef BEAGLE
1101 children (nn, &d1, &d2);
1102 if (d1 == oldnn)
1103 {
1104 d1 = d2;
1105 }
1106 id0 = nn->id;
1107 id1 = new_id(proposal->oback->id,proposal->world->sumtips);
1108 bid1 = proposal->oback->bid;
1109 v1 = proposal->v;
1110 id2 = d1->id;
1111 bid2 = d1->bid;
1112 v2 = d1->v;
1113 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
1114 proposal->v = nn->v;
1115 oldnn = nn;
1116 if(nn != proposal->connect)
1117 {
1118 nn = showtop (crawlback (nn));
1119 #endif
1120 while (nn != proposal->connect)
1121 {
1122 children (nn, &d1, &d2);
1123 if (d1 == oldnn)
1124 {
1125 d1 = d2;
1126 d2 = oldnn;
1127 }
1128 #ifdef BEAGLE
1129 id0 = nn->id;
1130 id1 = new_id(oldnn->id,proposal->world->sumtips);
1131 bid1 = oldnn->bid;
1132 v1 = oldnn->v;
1133 id2 = d1->id;
1134 bid2 = d1->bid;
1135 v2 = d1->v;
1136 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
1137 #else
1138 pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v, d1->x, d1->scale, d1->v);
1139 #endif
1140 proposal->v = nn->v;
1141 oldnn = nn;
1142 nn = showtop (crawlback (nn));
1143 }
1144 if (obb != proposal->connect)
1145 {
1146 #ifdef BEAGLE
1147 id0 = nn->id;
1148 id1 = new_id(oldnn->id, proposal->world->sumtips);
1149 bid1 = oldnn->bid;
1150 v1 = proposal->v;
1151 id2 = new_id(obb->id, proposal->world->sumtips);
1152 bid2 = obb->bid;
1153 v2 = obb->v;
1154 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
1155 #else
1156 pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
1157 proposal->xt, proposal->mt, proposal->vs);
1158 #endif
1159 }
1160 #ifdef BEAGLE
1161 }
1162 #endif
1163 proposal->v = proposal->connect->v;
1164 #ifdef BEAGLE
1165 evaluate_beagle_instances_proposal (proposal,
1166 showtop (crawlback (proposal->connect)), proposal->connect,
1167 proposal->connect->id, proposal->connect->bid,
1168 proposal->v);
1169 #else
1170 pseudoevaluate (proposal, proposal->xf, proposal->mf,
1171 showtop (crawlback (proposal->connect)), proposal->connect,
1172 proposal->v);
1173 #endif
1174 }
1175
1176 void
t_s_tcon(proposal_fmt * proposal,node * uppern)1177 t_s_tcon (proposal_fmt * proposal, node * uppern)
1178 {
1179 node *obb, *d1, *d2;
1180 children (proposal->target, &d1, &d2);
1181 if (d1 == uppern)
1182 {
1183 d1 = d2;
1184 d2 = uppern;
1185 }
1186 proposal->vs = fabs (proposal->target->tyme - showtop (uppern)->tyme);
1187 #ifdef BEAGLE
1188 long id0, id1, id2;
1189 long bid1, bid2;
1190 double v1, v2;
1191 id0 = proposal->target->id;
1192 id1 = new_id(uppern->id, proposal->world->sumtips);
1193 bid1 = uppern->bid;
1194 v1 = proposal->vs;
1195 id2 = d1->id;
1196 bid2 = d1->bid;
1197 v2 = d1->v;
1198 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
1199 #else
1200 pseudonuview (proposal, proposal->xt, proposal->mt, proposal->vs, d1->x, d1->scale, d1->v); /* eval target */
1201 #endif
1202 proposal->vs = proposal->time - proposal->target->tyme;
1203 #ifdef BEAGLE
1204 id0 = proposal->oback->id;
1205 id1 = proposal->origin->id;
1206 bid1 = proposal->origin->bid;
1207 v1 = proposal->v;
1208 id2 = proposal->target->id;
1209 bid2 = proposal->target->bid;
1210 v2 = proposal->vs;
1211 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
1212 #else
1213 pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v, proposal->xt, proposal->mt, proposal->vs); /*eval oback */
1214 #endif
1215 obb = showtop (crawlback (proposal->target));
1216 proposal->v = fabs (obb->tyme - proposal->time);
1217 if (obb->type == 'r')
1218 return;
1219 else
1220 {
1221 children (obb, &d1, &d2);
1222 if (d1 == proposal->target)
1223 {
1224 d1 = d2;
1225 d2 = proposal->target;
1226 }
1227 #ifdef BEAGLE
1228 id0 = obb->id;
1229 id1 = new_id(proposal->oback->id, proposal->world->sumtips);
1230 bid1 = proposal->oback->bid;
1231 v1 = proposal->v;
1232 id2 = d1->id;
1233 bid2 = d1->bid;
1234 v2 = d1->v;
1235 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
1236 evaluate_beagle_instances_proposal (proposal,
1237 showtop (crawlback (obb)), obb, obb->id, obb->bid,
1238 showtop (crawlback (obb))->tyme - obb->tyme);
1239 #else
1240 pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v, d1->x,
1241 d1->scale, d1->v);
1242 pseudoevaluate (proposal, proposal->xf, proposal->mf,
1243 showtop (crawlback (obb)), obb, obb->v);
1244 #endif
1245 }
1246 }
1247
1248
1249 /*-----------------------------------------------------------------
1250 erase all migration nodes between proposal->origin and proposal->oback
1251 */
1252 long
erase_migr_nodes(world_fmt * world,node * up)1253 erase_migr_nodes (world_fmt *world, node * up)
1254 {
1255 long deleted = 0;
1256 node *theNode, *down, *oldNode;
1257 down = crawlback (up);
1258 theNode = showtop(up->back);
1259 // printf("-------------\nup=%li\n",up->id);
1260 while (theNode->type == 'm')
1261 {
1262
1263 oldNode = theNode;
1264 theNode = showtop(theNode->back);
1265 free_mignodelet(oldNode,world);
1266 deleted++;
1267 //printf("-");
1268 }
1269 //printf("down=%li deleted=%li\n-------------\n",down->id, 2*deleted);
1270 down->back = up;
1271 up->back = down;
1272 //if(world->heat < 0.001)
1273 // printf("%i> heat=%f deleted=%li\n",myID,1./world->heat,deleted);
1274 return deleted;
1275 }
1276
1277 void
erase_migr_nodes2(world_fmt * world,node * up)1278 erase_migr_nodes2 (world_fmt *world, node * up)
1279 {
1280 int skip=0;
1281 long deleted = 0;
1282 node *theNode, *oldNode;
1283 //node *down
1284 if(skip!=0)
1285 return;
1286 //xcode down = crawlback (up);
1287 theNode = showtop(up->back);
1288 // printf("-------------\nup=%li\n",up->id);
1289 while (theNode->type == 'm')
1290 {
1291 oldNode = theNode;
1292 theNode = showtop(theNode->back);
1293 free_mignodelet(oldNode,world);
1294 deleted++;
1295 //printf("-");
1296 }
1297 }
1298 /*-------------------------------------------------------
1299 Connects and adjusts three nodes, the first is the
1300 mother node and the next two are child nodes, if one
1301 of the child nodes is NULL it just connects mother with
1302 the not-NULL child.
1303 */
1304 void
connectnodes(world_fmt * world,node * mother,node * brother,node * sister)1305 connectnodes (world_fmt *world, node * mother, node * brother, node * sister)
1306 {
1307 node *tmp;
1308 if(mother->top!=1)
1309 error("mother is not on top");
1310 if(brother !=0 && brother->top!=1)
1311 error("brother is not on top");
1312 if(sister!=NULL && sister->top!=1)
1313 error("sister is not on top");
1314 if ((sister != NULL) && (brother != NULL))
1315 {
1316 if ((mother == brother) || (mother == sister)
1317 || (sister == brother))
1318 {
1319 error ("connectnodes() conflict");
1320 }
1321 tmp = crawl_down (brother, mother->tyme);
1322 #ifdef MESS
1323 if(tmp->back != NULL && tmp->back->type == 'm' && showtop(tmp->back)->tyme > mother->tyme)
1324 {
1325 printf("new code 1 \n");
1326 erase_migr_nodes2(world,showtop(tmp));
1327 }
1328 #endif
1329 mother->next->back = tmp;
1330 tmp->back = mother->next;
1331
1332 tmp = crawl_down (sister, mother->tyme);
1333 #ifdef MESS
1334 if(tmp->back != NULL && tmp->back->type == 'm' && showtop(tmp->back)->tyme > mother->tyme)
1335 {
1336 printf("new code 2\n");
1337 erase_migr_nodes2(world,showtop(tmp));
1338 }
1339 #endif
1340
1341 mother->next->next->back = tmp;
1342 tmp->back = mother->next->next;
1343 }
1344 else
1345 {
1346 if (sister == NULL)
1347 {
1348 tmp = crawl_down (brother, mother->tyme);
1349 #ifdef MESS
1350 if(tmp->back != NULL && tmp->back->type == 'm' && showtop(tmp->back)->tyme > mother->tyme)
1351 {
1352 printf("new code 3\n");
1353 erase_migr_nodes2(world,showtop(tmp));
1354 }
1355 #endif
1356
1357 mother->next->back = tmp;
1358 tmp->back = mother->next;
1359 }
1360 else
1361 {
1362
1363 if (brother == NULL)
1364 {
1365 if (mother->type == 'm')
1366 {
1367 tmp = crawl_down (sister, mother->tyme);
1368 #ifdef MESS
1369 if(tmp->back != NULL && tmp->back->type == 'm' && showtop(tmp->back)->tyme > mother->tyme)
1370 {
1371 printf("new code 4 \n");
1372 erase_migr_nodes2(world,showtop(tmp));
1373 }
1374 #endif
1375 mother->next->back = tmp;
1376 tmp->back = mother->next;
1377 }
1378 else
1379 {
1380 tmp = crawl_down (sister, mother->tyme);
1381 #ifdef MESS
1382 if(tmp->back != NULL && tmp->back->type == 'm' && showtop(tmp->back)->tyme > mother->tyme)
1383 {
1384 printf("new code 5\n");
1385 erase_migr_nodes2(world,showtop(tmp));
1386 }
1387 #endif
1388
1389 mother->next->next->back = tmp;
1390 tmp->back = mother->next->next;
1391 }
1392 }
1393 else
1394 {
1395 error ("Single, lonely rootnode detected in connectenodes()\n");
1396 }
1397 }
1398 }
1399 }
1400
1401 void
gotoroot(node * origin,node ** ptrlist)1402 gotoroot (node * origin, node ** ptrlist)
1403 {
1404 node *theNode;
1405 long i = 0;
1406
1407 for (theNode = origin;
1408 (theNode != NULL) && (crawlback (theNode)->type != 'r');
1409 theNode = showtop (crawlback (theNode)))
1410 {
1411 ptrlist[i++] = theNode;
1412 }
1413 ptrlist[i] = theNode; /* adds root->back to the list */
1414 ptrlist[i+1] = NULL; /* guarantees an end for the array)*/
1415 }
1416
1417 void
adjust(node * theNode,MYREAL tyme,long level)1418 adjust (node * theNode, MYREAL tyme, long level)
1419 {
1420 if (level < 0 || theNode == NULL)
1421 return;
1422
1423 theNode->tyme = tyme;
1424 if (theNode->type == 'r')
1425 {
1426 theNode->v = 1;
1427 theNode->length = MYREAL_MAX;
1428 if (theNode->next->back != NULL)
1429 adjust (crawlback (theNode->next), crawlback (theNode->next)->tyme,
1430 level);
1431 if (theNode->next->next->back != NULL)
1432 adjust (crawlback (theNode->next->next),
1433 crawlback (theNode->next->next)->tyme, level);
1434 }
1435 else if (theNode->type == 't')
1436 {
1437 //#ifndef TESTINGDATE
1438 //theNode->tyme = 0.0;
1439 //#endif
1440 theNode->length = lengthof (theNode);
1441 ltov (theNode);
1442 return;
1443 }
1444 else if ((theNode->type != 't'))
1445 {
1446 if (theNode->type != 'm')
1447 {
1448 theNode->length = lengthof (theNode);
1449 ltov (theNode);
1450 if (theNode->next->back != NULL)
1451 {
1452 theNode->next->tyme = tyme;
1453 theNode->next->v = theNode->v;
1454 theNode->next->length = theNode->length;
1455 adjust (crawlback (theNode->next),
1456 crawlback (theNode->next)->tyme, level - 1);
1457 }
1458 if (theNode->next->next->back != NULL)
1459 {
1460 theNode->next->next->tyme = tyme;
1461 theNode->next->next->v = theNode->v;
1462 theNode->next->next->length = theNode->length;
1463 adjust (crawlback (theNode->next->next),
1464 crawlback (theNode->next->next)->tyme, level - 1);
1465 }
1466 }
1467 else
1468 adjust (crawlback (theNode->next), crawlback (theNode->next)->tyme,
1469 level);
1470 }
1471 }
1472
1473
1474 /* calculates x-array down the tree assuming that only one line
1475 is affected by the change of the sub-likelihoods above that line
1476 BUT does NOT calculate the tree-likelihood as evaluate() does.
1477 THIS CHANGES THE TREE
1478 */
1479 void
localevaluate(node * mother)1480 localevaluate (node * mother)
1481 {
1482 node *nn = NULL;
1483
1484 if (mother->type != 'r')
1485 {
1486 set_dirty (mother);
1487 for (nn = mother; crawlback (nn)->type != 'r';
1488 nn = showtop (crawlback (nn)))
1489 {
1490 set_dirty (nn);
1491 }
1492 }
1493 }
1494
1495
1496 ///
1497 /// copy the sequence data
1498 void
copy_x(proposal_fmt * proposal,xarray_fmt xx1,xarray_fmt xx2)1499 copy_x (proposal_fmt * proposal, xarray_fmt xx1, xarray_fmt xx2)
1500 {
1501 long locus = proposal->world->locus;
1502 long endsite = proposal->endsite;
1503 long rcategs = proposal->world->options->rcategs;
1504 long *maxalleles = proposal->world->data->maxalleles;
1505
1506 #ifdef ALTIVEC
1507 const size_t sitelikesize = sizeof (FloatVec) * rcategs;
1508 #else
1509 const size_t sitelikesize = sizeof (sitelike) * rcategs;
1510 #endif
1511 long i;
1512 // long j;
1513
1514 switch (proposal->datatype)
1515 {
1516 case 'a':
1517 case 'b':
1518 case 'm':
1519 memcpy (xx1.a, xx2.a, sizeof (MYREAL) * maxalleles[locus]);
1520 break;
1521 case 's':
1522 case 'n':
1523 case 'h':
1524 case 'u':
1525 case 'f':
1526 for (i = 0; i < endsite; i++)
1527 {
1528 //for (j = 0; j < rcategs; j++)
1529 //{
1530 memcpy (xx1.s[i], xx2.s[i], sitelikesize);
1531 //}
1532 }
1533 break;
1534 }
1535 }
1536
1537 void
fix_root_pop(world_fmt * world,node * p)1538 fix_root_pop (world_fmt *world, node * p)
1539 {
1540 if (crawlback (p) != p->back)
1541 {
1542 erase_migr_nodes (world, p);
1543 }
1544 if (p->back->actualpop != p->actualpop)
1545 {
1546 p->back->pop = p->back->next->pop = p->back->next->next->pop = p->pop;
1547 p->back->actualpop = p->back->next->actualpop = p->actualpop;
1548 p->back->next->next->actualpop = p->actualpop;
1549 }
1550 }
1551
1552
1553 /* transfers an x-array down the tree assuming that only one line
1554 is affected by the change of the sub-likelihoods above that line
1555 BUT does NOT calculate the tree-likelihood as evaluate() does.
1556 DOES NOT CHANGE THE TREE
1557 */
1558 void
pseudoevaluate(proposal_fmt * proposal,xarray_fmt x,MYREAL * lx,node * mother,node * newdaughter,MYREAL v)1559 pseudoevaluate (proposal_fmt * proposal, xarray_fmt x, MYREAL *lx,
1560 node * mother, node * newdaughter, MYREAL v)
1561 {
1562 node *nn = NULL, *d1 = NULL, *d2 = NULL, *oldnn = NULL;
1563
1564 if (mother->type != 'r')
1565 {
1566 children (mother, &d1, &d2);
1567 if (d1 == newdaughter)
1568 {
1569 d1 = d2;
1570 d2 = newdaughter;
1571 }
1572 #ifdef BEAGLE
1573 long id0, id1, id2;
1574 long bid1, bid2;
1575 double v1, v2;
1576 id0 = mother->id;
1577 id1 = new_id(proposal->oback->id,proposal->world->sumtips);
1578 bid1 = proposal->oback->bid;
1579 v1 = v;
1580 id2 = d1->id;
1581 bid2 = d1->bid;
1582 v2 = d1->v;
1583 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
1584 #else
1585 pseudonuview (proposal, x, lx, v, d1->x, d1->scale, d1->v);
1586 #endif
1587 oldnn = mother;
1588 nn = showtop (crawlback (mother));
1589 while (nn->type != 'r')
1590 {
1591 children (nn, &d1, &d2);
1592 if (d1 == oldnn)
1593 {
1594 d1 = d2;
1595 d2 = oldnn;
1596 }
1597 #ifdef BEAGLE
1598 id0 = nn->id;
1599 id1 = new_id(d2->id,proposal->world->sumtips);
1600 bid1 = d2->bid;
1601 v1 = d2->v;
1602 id2 = d1->id;
1603 bid2 = d1->bid;
1604 v2 = d1->v;
1605 prepare_beagle_instances_proposal(proposal, id0, id1, bid1, v1, id2, bid2, v2, proposal->world->beagle);
1606 #else
1607 pseudonuview (proposal, x, lx, oldnn->v, d1->x, d1->scale, d1->v);
1608 #endif
1609 oldnn = nn;
1610 nn = showtop (crawlback (nn));
1611 }
1612 }
1613 }
1614
1615 node *
findcrossing(node ** ptrl1,node ** ptrl2)1616 findcrossing (node ** ptrl1, node ** ptrl2)
1617 {
1618 long i = 0, j = 0;
1619
1620 /* assumes that there is an NULL element at the end */
1621 for (i = 0; ptrl1[i] != NULL; j = 0, i++)
1622 {
1623 while ((ptrl2[j] != NULL) && (ptrl1[i] != ptrl2[j]))
1624 j++;
1625 if (ptrl2[j] != NULL)
1626 {
1627 break;
1628 }
1629 }
1630 return ptrl1[i];
1631 }
1632
1633 node *
crawl_down(node * theNode,MYREAL tyme)1634 crawl_down (node * theNode, MYREAL tyme)
1635 {
1636 node *tN = showtop(theNode);
1637 node *otmp, *tmp = tN->back;
1638
1639 otmp = tN;
1640 if (tmp == NULL)
1641 return otmp;
1642 while ((tmp->type == 'm') && showtop (tmp)->tyme < tyme)
1643 {
1644 #ifdef TREEDEBUG
1645 printf("Crawldown: start=%x<%x> (%li<%li>) %li %f %f\n",theNode, tN, theNode->id, tN->id, tmp->id, showtop (tmp)->tyme, tyme);
1646 #endif
1647 otmp = tmp->next;
1648 tmp = tmp->next->back;
1649 if (tmp == NULL)
1650 return otmp;
1651 }
1652 return otmp;
1653 }
1654