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